メタゲノムアドベントカレンダー、最終日の記事。ぎりぎりの登録&投稿ですみません...
マイクロバイオームとメタボロームの相関を検出する手法、microbe–metabolite vectors、略してmmvecの使い方の紹介。
(Morton, et al. "Learning representations of microbe–metabolite interactions." Nature Methods 16, 1306–1314 (2019))
word2vecと語感をそろえてきてる（？）ことからわかるように、マイクロバイオームとメタボロームのペアデータセットから、それらの関係性を表現した潜在ベクトルを学習する手法。
問題設定
マイクロバイオーム解析において、相関はすごく難しい。相関は因果を意味しない、とかって話ではなく、そもそもそれ以前の問題として、マイクロバイオームデータを別のメタデータやオミックスデータと相関分析して出てきた相関関係は信用ならない場合が多い。
根っこの問題は、マイクロバイオームが本質的に組成データ(Compositional data)である（として扱うことが多い）ため。
系統組成、遺伝子機能組成など、相対存在量として記述されたデータで単純な相関係数を計算すると、偽陽性が非常に多く検出されやすいことが知られている。
たとえば、単純な例として、ある環境メタデータの増大に伴うマイクロバイオームの変動を時系列解析したとする。ある系統が環境メタデータの増大に応答してすばやく増殖し、他のほとんどの種はゆっくりと増殖したとする。このデータを相対存在量に変換してピアソンやスピアマン相関係数を計算すると、ある系統は環境メタデータと正に相関して増殖し、他のほとんどの系統は環境メタデータと負に相関して減少していくように見えてしまう（組成データにおける和が一定の制約＝100%のパイを奪い合ってしまうため）。その結果大量の偽陽性が検出されることになる。
この問題はこれまでも、それこそピアソンの時代から繰り返し繰り返し議論されてきて、いくつかの対処法が提案されている。SparCC, Proportionality, SPIEC-EASIなど。
しかし、SparCCなど存在量のlog ratio変換に依存した方法は、マイクロバイオームデータの中だけで相関分析する場合など、単一のデータセットの中で計算するぶんには問題ないのだけど、トータルのカウントが大きく異なる複数のデータセットのあいだで計算する場合は、log ratio変換の前提が成り立たないため妥当な計算とはならない。Meta16SとFungal ITSのあいだ、くらいのスケールの違いならまだマシだけど、Meta16Sとメタボロームのあいだだと、スケールがまるで異なるため深刻な問題になる。
マイクロバイオームとメタボロームのあいだの相互作用は研究の題材としてすごく重要なんだけど、「相関」が正しく計算できない。じゃあどうしたらいいか、ということに取り組んだのがこの論文。
※この手法が既存手法と比較して本当に検出精度が高いかどうか、まだいろいろと議論があるっぽいので注意。
反論（適切に変換すればピアソンなどの線形な手法のほうがmmvecよりずっと精度が高い、という主張）
反論への反論（テストしてるデータセットがダメ、実データはすごくスパースで、そういうデータでテストするとやっぱりmmvecが勝つ、と反論）
mmvecのアルゴリズム
「相関」を計算することはやめて、微生物（$\mu$）が存在したときの代謝物質（$\nu$）の条件付き確率$P(\nu|\mu)$を推定しよう、というコンセプト。
上の図はmmvecのリポジトリより。
この図がちょっとミスリーディングで、オートエンコーダとかVAEっぽく見えるけど、実際はシンプルな確率モデルをtensorflowでMAP推定してるだけ。なんかやたらと、ニューラルネットワークで云々みたいなこと主張してるけど、これがニューラルネットワークモデルなのかどうかよくわからない（NNの定義がよくわからないのでなんともいえないけど、これがNNならロジスティック回帰もNNだと思う）
生成モデル
微生物$\mu \in \{1, ...N\}$、代謝物質$\nu \in \{1, ...M\}$のそれぞれについて、以下のように正規分布からp次元の微生物ベクトル$u_\mu$, 代謝物質ベクトル$v_\nu$を生成する。（低次元pはユーザが設定する）
$$u_\mu \sim \mathcal{N}(0, \sigma_\mu I)$$
$$v_\nu \sim \mathcal{N}(0, \sigma_\nu I)$$
マイクロバイオームデータからサンプリングされた配列の系統アノテーション$\mu$について、代謝物質$v$の条件付き確率を、
$$P(v|\mu) = \frac {exp \left( v_\nu \cdot u_\mu + v_{\nu0} + u_{\mu0} \right)} {\sum_j exp \left( v_j \cdot u_\mu + v_{j0} + u_{\mu0} \right)}$$
とする。サンプル$k$の代謝物質の存在量$y_k$（総数$n$）はこの多項分布$q$からのサンプリングで、
$$q = [ P(v_1|\mu), ...P(v_M|\mu)]$$
$$y_k \sim Multinomial(n, q)$$
とする。以上のシンプルなモデル。
これは、それぞれの微生物$\mu$ごとの代謝物質$\nu$の条件付き確率$P(\nu|\mu)$のテーブル（をCenter logratio変換したもの）を、微生物ベクトルを並べた行列$U$と、代謝物質ベクトルを並べた行列$V$に行列分解していることに相当する。
推論は、tensorflowのadam optimizerを使って、以下の対数尤度を最大化している。サンプル$k$の配列集合を$x_k$とすると、
$$L = L_y + L_U + L_V$$
$$L_y = \sum_k \sum_{r \in x_k} ln Multinomial(y_k | q)$$
$$L_U = \sum_\mu \sum_{p=1}^{p} ln \mathcal{N}(U_{\mu,p} | 0, \sigma_\mu)$$
$$L_V = \sum_\nu \sum_{p=1}^{p} ln \mathcal{N}(V_{\nu,p} | 0, \sigma_\nu)$$
TensorBoardで計算グラフを見てみると↓のような感じ。
マイクロバイオームデータ中のすべての配列一本一本について対数尤度を計算して、パラメータをアップデートしていく。
論文メソッド中の記述によると、こうやって学習された$u$ベクトル, $v$ベクトルはどちらもp次元の割合ベクトルに直接変換できる（Aitchison simplexからのcenter logratio変換とみなせる）ため、トピックモデルにおけるトピックの組成と同じような解釈ができるらしい。
土壌微生物データに適用してみる
インストールは単純に、pip install mmvecなど。
学習に結構時間がかかるので、GPUを使って計算する場合はtensorflow-gpu を入れておく。
QIIME2にプラグインとして統合する方法もあるみたいだけど、QIIME使ったことないのでよくわからない。GPUが使えなかったり制限があるみたいなのでスタンドアロン版使った方がいいと思う。
題材として扱うのは、mmvec論文のFig.3でも扱われている土壌微生物群集のメタゲノム・メタボロームデータ。
Swenson, et al. "Linking soil biology and chemistry in biological soil crust using isolate exometabolomics." Nature Communications 9, 19 (2018)
mmvec論文のFig.3aの再現を試みる。
まずはデータの読み込み。
すごくありがたいことに、土壌微生物論文はSupplementary Data 1, Supplementary Data 4で19サンプルの系統組成、代謝物質組成をエクセルファイルとして配布してくれているので、そのままロードして少しだけ整形してtsvのテーブルファイルとして保存する。
まずはマイクロバイオームデータ。
import pandas as pd
microbe = pd.read_excel('./data/41467_2017_2356_MOESM7_ESM.xlsx', \
skiprows=0, header=1, index_col=0)
microbe = microbe.drop(['phylum', 'class', 'order', 'family', 'genus', 'species'], axis=1)
microbe = microbe.drop('9hr_late', axis=1) # No corresponding sample in metabolites
microbe = microbe * 1000.
microbe.to_csv('./data/microbe.tsv', sep='\t')
print(microbe.shape)
print(microbe.iloc[:3, :3])
(468, 19)
3min_early 3min_earlymid 3min_latemid
rplO ID
rplo 1 (Cyanobacteria) 4121.0 4489.0 2032.0
rplo 10 (Firmicutes) 0.0 0.0 0.0
rplo 100 (Proteobacteria) 132.0 56.0 16.0
19サンプル、検出された系統は468。
次にメタボロームデータ。
metabolite = pd.read_excel('./data/41467_2017_2356_MOESM4_ESM.xlsx', \
skiprows=0, header=1, index_col=0, \
usecols='A:AG')
metabolite = metabolite.iloc[:, 13:]
# 微妙にサンプルのラベルが違うので、マイクロバイオームデータに合わせる
metabolite = metabolite.rename(columns=lambda l: l.split('_')[0] +'_'+ l.split('_')[2])
metabolite = metabolite.loc[:, microbe.columns]
metabolite.to_csv('./data/metabolite.tsv', sep='\t')
print(metabolite.shape)
print(metabolite.iloc[:3, :3])
(85, 19)
3min_early 3min_earlymid 3min_latemid
Name
(2,3-dihydroxy-3-methylbutanoate) 1.00 4.609149e+02 1.379355e+02
(2,5-diaminohexanoate) 2782242.25 1.883228e+06 2.531814e+06
(3-hydroxypyridine) 2764132.25 3.014544e+06 3.817224e+06
19サンプル、検出している代謝物質は85種類。
この２つのファイルを入力として、mmvec paired-omicsコマンドで学習を実行する。
以下、オプション：
--microbe-file, --metabolite-file: 系統組成データ、メタボロームデータ。リポジトリのREADMEではbiom-formatのファイルしか入力できないっぽく書かれてるけど、わざわざbiomに変換しなくても普通のタブ区切りテキストファイルを指定して実行できる。
--summary-dir: 結果を出力するディレクトリ
--epochs: エポック数
--batch-size: ミニバッチのサイズ
--latent-dim: 潜在ベクトルの次元
--input-prior: 事前分布のハイパーパラメータ$\sigma_\mu$
--output-prior: 事前分布のハイパーパラメータ$\sigma_\nu$
--arm-the-gpu: GPU使う場合はオンにする
--learning-rate, --beta1, --beta2: 学習率、Adamのパラメータ
mmvec paired-omics \
--microbe-file ./data/microbe.tsv \
--metabolite-file ./data/metabolite.tsv \
--summary-dir ./summary
エポックをどの程度にすればいいかは、データセットによって大きく異なるので一概には言えない。
さいわい、RMSEや対数尤度を記録してくれるので、実行した場所でtensorboard --logdir .すればTensorBoard上で収束の確認はできる。今回の場合↓
ls summary
latent_dim_3_input_prior_1.00_output_prior_1.00_beta1_0.90_beta2_0.95
latent_dim_3_input_prior_1.00_output_prior_1.00_beta1_0.90_beta2_0.95_embedding.txt
latent_dim_3_input_prior_1.00_output_prior_1.00_beta1_0.90_beta2_0.95_ordination.txt
latent_dim_3_input_prior_1.00_output_prior_1.00_beta1_0.90_beta2_0.95_ranks.txt
学習が終わると3種類のファイルを出力する。
それぞれのファイルがなんなのか、ろくに説明がないのでわかりにくいが、ソースコードを読む限り以下のような結果が書き出されているらしい。
** XXX_embedding.txt
- それぞれの微生物、代謝物ごとの潜在ベクトル（u, v）
- XXX_ranks.txt
- 行列U, 行列Vを掛け算して中心化したマトリックス。条件付き確率$P(\nu|\mu)$をCenter logratio変換した値と一致。
- XXX_ordination.txt
- 上のranksファイルのマトリックスを特異値分解（u, s, v）して計算した、微生物の埋め込み（dot(u,s)）と代謝物の埋め込み（v）。論文中のbiplotはこの埋め込みで描かれている。
たぶん一番使い勝手がいいのはXXX_ranks.txtファイルで、普通の相関分析と同じように扱える。ただし、負の値は負の相関を意味するわけではなく、「全代謝物質の平均的な条件付き確率よりも下回っている」ことを意味している点に注意。
rankfile = './summary/latent_dim_3_input_prior_1.00_output_prior_1.00_beta1_0.90_beta2_0.95_ranks.txt'
ranks = pd.read_csv(rankfile, sep='\t', index_col=0)
print(ranks.shape)
print(ranks.iloc[:3, :3])
(85, 191)
rplo 1 (Cyanobacteria) rplo 100 (Proteobacteria) rplo 101 (Proteobacteria)
featureid
(2,3-dihydroxy-3-methylbutanoate) -5.076497 -0.229192 -3.435717
(2,5-diaminohexanoate) -1.053542 -1.119804 -0.508942
(3-hydroxypyridine) 0.055645 0.001804 0.797657
ここでは（mmvec論文同様）、Microcoleus vaginatus ("rplo 1 (Cyanobacteria)") の結果に焦点をあてる。
もとの土壌微生物論文では、Microcoleus vaginatusの分離株の実験で、以下の代謝物を産生・放出することを確認していた。
しかし、メタゲノムのSpearman相関分析の結果、この中の一部については正の相関係数が検出されたが、半分は負の相関となってしまい、分離株と環境中で矛盾した結果となってしまった（論文のFig.3）。
true_metabolites_text_offsets = {
'(3-methyladenine)': (-0., -0.15),
'7-methyladenine': (0.7, -0),
'4-guanidinobutanoate': (-1.9, -0.2),
'uracil': (0.6, -0.1),
'xanthine': (0, 0.1),
'hypoxanthine': (0.2, -0.2),
'(N6-acetyl-lysine)': (-1.7, 0.1),
'cytosine': (0, 0.1),
'N-acetylornithine': (0.2, 0.1),
'succinate': (-0.5, 0.2),
'adenosine': (-1.2, 0.1),
'guanine': (0.4, -0.2),
'adenine': (-0.7, -0.2)}
これらの代謝物に関して、mmvecによる条件付き確率の推定ではどのような値になっているか調べてみる。
target = 'rplo 1 (Cyanobacteria)' # Microcoleus vaginatus
Microcoleus = pd.DataFrame(ranks.loc[:, target].values, \
index=ranks.index, columns=['log CondProb'])
%matplotlib inline
import matplotlib.pyplot as plt
import seaborn as sns
df = Microcoleus.sort_values(by='log CondProb', ascending=False)
colors = ['red' if m in true_metabolites_text_offsets.keys() else 'gray' for m in df.index]
fig, ax = plt.subplots(figsize=(16, 6))
g = sns.barplot(x=df.index, y='log CondProb', \
palette=colors, data=df, ax=ax)
g.set_xticklabels(g.get_xticklabels(), rotation=90)
ax.set_ylabel('log conditional probabilities')
sns.despine()
赤色がMicrocoleus vaginatusが放出するはずの代謝物。
mmvecによる推定値では、元論文で負のSpearman相関となってしまった代謝物含めて、すべてプラスの値=平均以上の確率値としてちゃんと検出できていることがわかった。
それぞれの代謝物に関してSpearman相関も計算して、結果を比較してみる。
import scipy.stats
Microcoleus['Spearman rho'] = np.zeros(len(Microcoleus))[:, np.newaxis]
Microcoleus['P values'] = np.zeros(len(Microcoleus))[:, np.newaxis]
for m in Microcoleus.index:
rho, p = scipy.stats.spearmanr(microbe.loc[target, :].values, \
metabolite.loc[m, :].values)
Microcoleus.loc[m, 'Spearman rho'] = rho
Microcoleus.loc[m, 'P values'] = p
fig, ax = plt.subplots(figsize=(12, 12))
ax.scatter(Microcoleus['log CondProb'], Microcoleus['Spearman rho'], \
s=-10*np.log(Microcoleus['P values'].values))
for m in true_metabolites_text_offsets.keys():
x, y, _ = Microcoleus.loc[m, :].values
dx, dy = true_metabolites_text_offsets[m]
fc = 'red' if y < 0 else 'black'
ax.annotate(m, xy=(x, y), xytext=(x+dx, y+dy), fontsize=16, \
arrowprops=dict(facecolor=fc, shrink=0.01))
ax.axhline(y=0, color='gray', linestyle=':')
ax.axvline(x=0, color='gray', linestyle=':')
ax.set_xlabel('log conditional probabilities', fontsize=18)
ax.set_ylabel('Spearman correlations', fontsize=18)
plt.show()
円の大きさはSpearman相関のp-valueを示している。矢印がついているのが、Microcoleus vaginatusが放出するはず（正に相関するはず）の代謝物。赤色矢印が、分離株の結果とSpearman相関で矛盾となってしまっている代謝物。
結果、やはりSpearman相関ではいくつかの代謝物に関して負の相関係数が計算されてしまっているが、mmvecの計算ではいずれもプラスとして計算できていることがわかる。
「正解」がなんなのか定義するのが難しく、手法の有効性をどのように評価するのか、まだいろいろと議論がされているが、この結果を見る限りはそれなりに有望なアプローチっぽい。