読者です 読者をやめる 読者になる 読者になる

ピエール瀧になりたい

備忘録として

メタゲノムデータ解析その1

統計の勉強がてらTara Oceansのメタゲノムデータ解析をしています。

Tara OceansとはTara号という小さな研究船で世界中の海を旅してその生物や物理化学を調べるプロジェクトで、
微生物群集に関してのメタゲノムデータは以下のサイトで割といい感じに公開されています。

Tara Oceans

せっかく公開データを使った解析をやるのだからブログに進捗をちょいちょい上げていこうと思います。
僕の理解が違っているところもあると思うのでコメントもらえたほうが勉強になるだろうし。

というわけで初回はまずTara Oceansの微生物群集に関する以下の論文に触れつつ解析の目標などを明確にしていきたいと思います。

Sunagawa, Shinichi, et al. "Structure and function of the global ocean microbiome." Science 348.6237 (2015): 1261359.

f:id:kmooog:20170414072852j:plain

Tara Oceansでは世界中の海から243の微生物サンプルを取っています。そのうち139が原核生物由来なので今回はそれを使います。
それぞれのサイトで表層 (SRF)、 深部クロロフィル極大 (DCM)、 中深海水層 (MESO)でサンプルを取得します(3点で取れていないサイトもあります)。

いろいろ他にも面白い結果を出しているのですがレビューが目的では無いので関係する部分だけ触れます。

f:id:kmooog:20170414074116p:plain

なかなかかっこいいビジュアライズですね。
右の図は各環境パラメータ間のスピアマン相関係数を算出しています。
一例として、リン酸濃度と硝酸濃度は強い正の相関を示し、温度と酸素濃度は強い逆相関を示しています。

そして右に表される相関係数の行列と、微生物群集のデータとの間での偏マンテル検定の結果が左側のエッジです。

偏マンテル検定は生態学でよく使われる統計手法のようですが初見でした。
以下に詳しく書いてあって参考になりました。

Mantel's Test

目的としては環境パラメータから種の分布を説明するときに使います。

環境パラメータは大体相互に影響し合っている(例えば光強度が強ければ光合成が活発になり酸素濃度が増加するだろう)ので一つの変数で因果関係を説明することが困難です。また、時系列データの場合などには自己相関があるので、通常の相関分析が適応できないそうです。

まず、自己相関がある場合、通常の線形回帰を行うことができません。回帰関係を


y_i = \alpha + \beta x_i + \epsilon_i

としたときに、誤差項である\epsilon_i\epsilon_i-1に依存してしまうからです。

具体的な例としては、例えば光強度と植物の量の関係を見たいときに、たまたまある点で調査する直前に山火事があって植物が全部燃えたとします。これは\epsilon_iに関してはただの誤差(または外れ値)として現れますが、同地点での次の調査や地理的に近い地点での調査である\epsilon_i+1の値にも影響を与えてしまいます。

通常の線形回帰分析は誤差が独立に生じることを前提とするので、これではまずいのです。
相関係数の計算に関しても同様の仮定をおいているのでおかしくなります。

ではどうするか。マンテル検定ではこの問題を解決するために、
二つの距離行列を解析対象にします。
具体的には、「生物の遺伝的距離と空間的な距離は関係あるのか?」というようなことを解析対象とします。

統計量としては

\displaystyle  z = \sum_{i = 1}^N  \sum_{j = 1}^N x_{ij} y_{ij}

を用います。実際に用いるのはこれを正規化したものですが理解のためには簡単な方でいいでしょう。
ここでNはサンプリング地点の数で、xyはパラメータです(種の相関行列と環境パラメーターの相関行列)。

地点1と地点2、地点1と地点3.....におけるある生物aの遺伝的距離と空間的な距離の積を取って足したものですね。
距離行列の相関が高ければzは大きくなり、距離行列が逆相関すればzは負の値をとるような感じです。


今回の論文では、地点間の微生物群集の類似度、環境パラメータの類似度をユークリッド距離で算出し、
さらに空間分布を考慮した手法を用いたらしいです(偏マンテル検定)。
偏マンテル検定は「偏相関のように特定の変数の影響を排除したマンテル検定」だそうですが結構理論的にもめてる分野らしいのであまり深く突っ込まないことにしました。

空間分布は単純に球面状での距離であるHaversine distancesを用いており、海流の影響とか全く考えてないのでそれでいいのか感がちょっとあります。自分が解析するとしたら、そこはむしろ考慮しない気がします(微生物群集は高等生物に比べ地理的隔離の影響がはるかに小さいことが予想されるためです)。
やるのだとしたらまずHaversine distancesで表される空間分布と生物層に相関があるのかからだと思うのですが、
やってないとは思えないので事前の解析で色々やった結果この手法に落ち着いたか、レビューに言われたかみたいな感じでしょうか。

まぁそういうことを言い出すと地点間の微生物群集の類似度を単純なユークリッド距離で計算しているのはどうなのかとかも気になるところですが、割愛します。

無駄に長くなりましたが肝心な点は、微生物群集は温度、酸素濃度と非常に強い相関を持っており、
温度と酸素濃度が微生物群集を決める鍵だということです。
これ自体はずっと昔から言われてることですが、データ量が多くてビジュアライズがかっこいいので正義です。


f:id:kmooog:20170414105719p:plain


そして本題はこの図です。
この図では、環境変数である温度と酸素濃度を微生物群集から予測しています。
手法としてはelastic net modelを用いています。

最近ではメタゲノムにもよく使われる統計手法のようですね。
回帰問題を解く際の誤差関数の正則化項の手法を工夫し、
次元圧縮を行えるようにしたもののようです。

単一の観測変数から予測する場合は場合は以下のサイトがわかりやすかったです。
breakbee.hatenablog.jp

しかし、メタゲノムを用いて予測する場合は各OTUの存在量からなる多数の観測変数からの予測になるので、その場合どうするのでしょう。
各OTUに関してelastic net modelを用いて回帰曲線を描いて、それらを重み付けして合計するのか(非線形重回帰分析)、
それとも各OTUは全て1次に相関するものと仮定して計算するのか(重回帰分析)、
などと考えていたのですが"Compositional data (see above) were normalized to ranks"と書いてあり、
順序尺度に変換してありました。これで1次関数での相関は仮定できますが、大分情報量をロスしてるのが気になるところです。
Cross-validationをしているので計算量的にこのあたりが限界なのかもしれません。

図の説明に入ります。

elastic net modelを使って学習すると、かなり精度よくメタゲノムデータから水温を推定できるというのが右の図です。これは結構かっこいいです。
External validationとして、異なるプロジェクトであるGrobal Ocean Samplingのデータを使っています。

そして左の図は、表層データで学習させて作ったモデルを深部クロロフィル極大に適応しても結構精度よく水温が推定できます。
酸素濃度は無理です。だから温度が一番大事です。という図です。

面白いし納得できるデータだと思います。

面白かったので自分でもデータをいじってみたくなったので、

・手法改良で温度や酸素濃度の予測精度をもうちょっと上げられないか?
・他の環境変数を遺伝子データからなんとか予測できないか?
・個々の遺伝子は環境パラメータとどのような相関関係を持つのか?

あたりのことを今後ちまちまやっていこうと思います。

前置きのつもりだったんですがダラダラ書いてたら長くなってしまいました。。。