こちらは創薬 Advent Calendar 2017（#souyakuAC2017）８日目の記事です。
創薬 Advent Calendar 2018（#souyakuAC2018）に前日譚である「テンソル分解を用いて遺伝子発現プロファイルからインシリコ創薬（Ⅱ）」を投稿しました。

インシリコ創薬の課題

いわゆるインシリコ創薬（コンピュータを用いた創薬）、バズワード的にはAI創薬とも最近は言われるが、は、徐々に実際の創薬につながるものも出てきたものの、実際に従来の実験ベースの創薬を淘汰したり、あるいは、せめて互角であったりする状態にあるとはなかなか言えない。

インシリコ創薬の２大潮流

最も多く使われているインシリコ創薬の潮流には構造ベースのものとリガンドベースのものがある。構造ベースのものは、標的となるタンパクの立体構造（これ自体が未知の場合には計算機で予測しなくてはならないが）に結合する低分子化合物をいわゆるドッキングシミュレーションで見つけるものである。この方法は非常に有望視されているがいろいろ問題点も多い。例えば

ドッキングシミュレーションに膨大な計算機資源が必要
タンパクの立体構造の情報が必要（未知の場合にはこれ自体が予測対象）
そもそも低分子化合物とタンパクの関係は量子力学の多体問題であり実効的に計算不能。近似をする必要があり、それが計算精度を下げる

など。一方、リガンドベースの場合は、タンパクの立体構造の情報は不要である。既知の低分子化合物があった場合、その化合物と「似た」化合物を探索することでインシリコ創薬を実行する。このため、上記のような問題は一切存在しない。しかし、別の問題も存在していて、それは「似た」化合物しかみつからない、ということである。あるタンパクに結合する低分子化合物は同じような構造をしているとは限らないし、そもそも、結合部位が違う場合には構造が近い必要さえない。さらに、もし、１つも有効な低分子化合物が知られていなかったら、そもそも、リガンドベースのインシリコ創薬は不可能である。構造ベースのインシリコ創薬は逆にこの様な問題点を持たないから、どちらかが優れているというよりは、むしろ、構造ベースとリガンドベースは相補的な関係にある。逆に言えば、これらの欠点のため、なかなか実際の創薬には届きにくい。

インシリコ創薬の新潮流

このような観点から第三の潮流が生まれた。それは遺伝子発現プロファイルからの創薬である。具体的には、実験動物や培養細胞になるべく多種類の薬剤を作用させて遺伝子発現プロファイルを記録しておく。そこで、「既知の」薬剤と同じ遺伝子発現プロファイルを持っている薬剤を探すことで創薬を行うのである。これは既知の薬剤が特定のたんぱくを標的とした薬剤であれば、そのたんぱくを標的とする別種の薬剤を探すことになるし、対象が疾患であれば、疾患に対する新たな創薬をすることになる。いずれにせよ、これらは「遺伝子発現プロファイル」を低分子化合物の特徴量であると考えて、低分子化合物を（効く薬剤と効かない薬剤に）分類している分類問題を統計学習あるいは機械学習で解く問題として扱っているとみなすことができる。この方法は低分子やタンパクの構造をまったく考えないため、第三の方法として有効であるが、しかし、既知の有効な薬剤が無い場合には実行しようがないというリガンドベースのインシリコ創薬と同じ問題を持っている。

教師なし学習によるインシリコ創薬

構造ベースの創薬はいわゆる教師無し学習であり、既知の有効な薬剤が無い場合でも実行可能であった。ここでは遺伝子発現プロファイルを使った教師なし学習によるインシリコ創薬の可能性について考えよう。

テンソル分解

かなり話が難しくなってしまうがここでテンソル分解といいうものを考える。テンソルとは添え字が３個以上ある行列のようなものである。たとえば、$i$番目の遺伝子の$j$番目の臓器における$k$番目の実験条件における発現量を$x_{ijk}\in \mathbb{R}^{N \times M \times K}$と書く、などという表記になる。テンソル分解ではこれを以下の様に展開する。
$$
x_{ijk} = \sum_{\ell_1 \ell_2 \ell_3} G(\ell_1,\ell_2,\ell_3) x_{\ell_1 i}x_{\ell_2 j} x_{\ell_3 k}
$$
ここで$ G(\ell_1,\ell_2,\ell_3) \in \mathbb{R}^{N \times M \times K}$はコアテンソル、$x_{\ell_1 i} \in \mathbb{R}^{N \times N}$,$x_{\ell_2 j} \in \mathbb{R}^{M \times M}$,$x_{\ell_3 k} \in \mathbb{R}^{K \times K}$,は特異値行列（直交行列であることを仮定する）。このあたりの詳細については関係データ学習の６．３タッカー分解などをみるとよいだろう。$x_{ijk}$と$G(\ell_1,\ell_2,\ell_3)$がまったく同じなのでこれではちっとも情報の縮約になっていないのだが、HOSVDというアルゴリズムを用いると絶対値の大きい少数個の$G$で打ち切っても精度よく近似できることが経験的に知られている。

テンソル分解を用いた遺伝子発現プロファイルの統合解析

$j$番目の化合物を作用させてから$k$番目の経過時間後の$i$番目の遺伝子の発現量を$x_{ijk}$とする。このようなデータが複数のラットの臓器に対して網羅的に蓄積されているデータセットがDrugMatrixである。これをなんらかの疾患の患者と健常者、計$S$人の遺伝子発現プロファイル($x_{is} \in \mathbb{R}^{ N \times S}$)と統合解析して、創薬候補化合物を推定することを考えよう。いろいろな方法が考えられるが、ここでは前述のテンソル分解を利用することを考えよう。そのために２つの行列の各成分の積から新たなテンソル$x_{ijks} = x_{ijk}x_{is}$を作成しよう。このテンソル$x_{ijks}\in \mathbb{R}^{N \times M\times K\times S}$をテンソル分解することで２つの遺伝子発現プロファイルを統合解析すると化合物で処理した時の遺伝子発現プロファイルと、疾患プロファイルの共通する遺伝子発現プロファイルが$x_{\ell_1 i}$として取り出されることになる。このテンソルあるいは行列の積から新しいテンソルを構成してテンソル分解を適用する一般的な枠組みについては文献に詳しく書いた。日本語での発表スライドとこのスライドの１ページにリンクしてある発表ビデオも参照してほしいが、脇道にそれ過ぎるのでここではテンソル分解の数理的な側面にはこれ以上詳しくは触れない。

テンソル分解を用いた教師なし学習による変数選択

我々が欲しいのは以下の様な遺伝子発現プロファイルを持つ遺伝子のリストである。

患者と健常者で発現に差がある
薬剤投与後の遺伝子発現プロファイルに時間依存的な変化がある。

これを実現するために$x_{ijks}$のテンソル分解
$$
x_{ijks} = \sum_{\ell_1 \ell_2 \ell_3 \ell_4} G(\ell_1,\ell_2,\ell_3,\ell_4) x_{\ell_1 i}x_{\ell_2 j} x_{\ell_3 k} x_{\ell_4 s}
$$
を考えよう。$x_{\ell_2 j}, x_{\ell_3 k}, x_{\ell_4 s}$はそれぞれ、化合物依存性、（投与後の）時間依存性、患者と健常者の差、を表現している。そこで、化合物特異的に、時間依存性及び患者と健常者の間の発現差を伴っている$\ell_2',\ell_3',\ell_4'$を探す。もちろん、教師無し学習の宿命でこのような成分が無いこともあり得る。以下ではそのような$\ell_2',\ell_3',\ell_4'$が見つかったとしよう。すると注目すべき(遺伝子の選択に使うべき）$x_{\ell_1' i}$は絶対値が大きな$G(\ell_1',\ell_2',\ell_3',\ell_4')$を伴う$\ell_1'$である。なぜなら、そのような$x_{\ell_1' i}$が元のテンソル$x_{ijks}$に対する寄与が最も大きいからである。

この様な$\ell_1'$が見つかったら、今度はこの中で特に（絶対）値が大きな$x_{\ell_1' i}$をもつ遺伝子$i$を見つける。これは以下の様な方法で行う。まず、帰無仮説として$x_{ijks}$が乱数であるという仮説を置き、これに反する$i$を探すことにする。この帰無仮説の元では$x_{\ell_1 i}$は中心極限定理のおかげで漸近的にガウス分布になることが期待される。そこで$x_{\ell_1 i}$がガウス分布に従っていると仮定し、その仮定では説明できないほど（絶対値が）大きな$x_{\ell_1 i}$をもっている$i$が$x_{ijks}$に有意に寄与しているとしよう。これを実行するために$x_{\ell_1 i}$に$\chi^2$分布（カイ二乗分布）を仮定して$P$値、$P_i$を付与する。
$$
P_i = P_{\chi^2} \left [> \sum_{\ell_1'} {\left ( \frac{x_{\ell_1' i}}{\sigma_{\ell_1'}}\right )^2} \right ]
$$
ただし$P_{\chi^2}[>x]$は引数が$x$より大きい場合の$\chi^2$分布による累積確率である。
付与された$P_i$は遺伝子が多数個あるために小さな$P$値が出やすくなってしまうという上述の問題を回避するため、Benjamini–Hochberg procedureで多重比較補正を行う。最後に補正された$P_i$が十分小さく（例えば０．０１以下）、ガウス分布には従っていないと思われる遺伝子を選択する。これがこの節の冒頭で述べた性質をもった遺伝子発現プロファイルが付与された遺伝子であると期待される。

標的タンパクの推定

これで薬剤によって発現が変化し、かつ、その変化が疾患によるものと一致していると期待できる遺伝子が特定できた。しかし、実際に化合物が結合しているのはタンパクである。一方、計測されているのはmRNAの発現量である。タンパクに化合物が結合することでそのタンパクをコードしている遺伝子のmRNAの量が変化しているとは思えないので、発現量が変化している遺伝子の中に標的タンパクは無いだろう。
化合物が実際に結合しているタンパクの他の遺伝子発現プロファイルに対する影響は、そのタンパクがコードされている遺伝子をノックアウトした場合に近いと期待される。そこで遺伝子を網羅的にノックアウトした場合の遺伝子発現プロファイルを参照することで標的遺伝子を推定する。これら一連の作業をイラスト化したのが以下になる。

計算例

この様な計算を実際におこなった（文献）。表２にその結果が書かれている。上から順に心臓病、心的外傷後ストレス障害（PTSD）、急性リンパ性白血病（ALL）、糖尿病（diabetes）、腎臓がん、肝硬変のそれぞれについて、ラットの心臓、脊髄、腎臓（糖尿病と腎臓がん）、肝臓に多種類の低分子化合物を作用させたときの遺伝子発現プロファイルと疾患における遺伝子発現プロファイルを比較して本手法で得られた低分子化合物と標的タンパクの組み合わせにどれくらい「当たり」（既知の低分子化合物―標的タンパクの関係）が含まれているかを$\chi^2$検定とフィッシャーの正確確率検定で調べた。ALL以外のすべてのケースで、予測された既知の低分子化合物―標的タンパク関係は既知のものと統計的に有意にオーバーラップしていた。この結果自体は既知の低分子化合物―標的タンパクの関係予測として考えた場合、必ずしも誇れるものではないが、完全な教師なし学習（疾患や標的タンパクに対する低分子化合物の既知の情報を全く使わない）であり、かつ、遺伝子発現プロファイル「だけ」から達成された成果であることを考えれば

教師データは不要だが膨大な計算機資源が必要な構造ベースのインシリコ創薬
膨大な計算資源は不要だが多くの教師データを必要とするリガンドベースや（従来の）遺伝子発現プロファイルベースのインシリコ創薬

とはまったく一線を画したインシリコ創薬であると言えるだろう。しかも、低分子化合物を作用した遺伝子発現プロファイルは任意の疾患遺伝子プロファイルを組み合わせることが可能なので低分子化合物の広範な標的疾患・標的タンパク予測に使えることが期待される。もっと言えば、低分子化合物を作用した遺伝子発現プロファイルの部分は核酸創薬でもペプチド創薬の様な中分子化合物でもなんでもよく、あるいは、なんらかの治療という意味で投薬である必要さえない。Aという疾患で有効だった治療法が別のBという疾患でも有効かどうかということを遺伝子発現プロファイルレベルで推定可能な方法だと言えるだろう。

結言

ここではテンソル分解を用いて複数の遺伝子発現プロファイルを統合解析し、教師なし学習でインシリコ創薬を行う試みを紹介した。複数の行列・テンソルの統合解析という意味ではこの枠組みはインシリコ創薬以外に、あるいは、ゲノム科学以外の広範なデータ解析に使用可能な手法である可能性もあるので、今後もいろいろな試みを続けていきたい。