細菌と代謝物の共起確率推定ツール mmvec の仕様

# 概要

本稿は、microbiome データと metabolome データを用いて細菌と代謝物の**共起確率（co-occurrence probabilities）**を推定する解析ツール [**mmvec**](https://github.com/biocore/mmvec) の詳細な仕様をまとめたものである。**mmvec** の理論背景などについては、以下の記事に大変わかりやすく解説されている。

https://qiita.com/khigashi02/items/48a216839b13af6c0826

# 原著論文

[Morton, et al. "Learning representations of microbe–metabolite interactions." Nature Methods 16, 1306–1314 (2019)](https://www.nature.com/articles/s41592-019-0616-3)

# 環境

* macOS Big Sur （v11.2.3） 

# mmvec のインストール

[公式ドキュメント](https://github.com/biocore/mmvec) にあるように、mmvec は standalone 版と Qiime2 版が用意されている。standalone 版は `` pip `` でインストールすることもできるが、バージョンの依存関係などからインストールできないことも多い。今回はパッケージのバージョンを管理するために ``conda`` で構築した仮想環境上にインストールする。

## 仮想環境の構築

mmvec は、菌叢解析プラットフォーム [Qiime2](https://qiita.com/keisuke-ota/items/6399b2f2f7459cd9e418) の [version 2020.6](https://docs.qiime2.org/2020.6/install/native/#install-qiime-2-within-a-conda-environment) をインストールする際に用意する仮想環境上にインストールすることができる。この仮想環境では、standalone 版と Qiime2 版の両方をインストールすることができる。[Qiime2 公式ドキュメント](https://docs.qiime2.org/2020.6/install/native/#install-qiime-2-within-a-conda-environment) に従い、まずは仮想環境を作る。

```terminal:ターミナル
## 仮想環境に使う yaml ファイルをダウンロードする
wget https://data.qiime2.org/distro/core/qiime2-2020.6-py36-osx-conda.yml 
## yaml ファイルをもとに仮想環境を構築する
conda env create -n qiime2-2020.6 --file qiime2-2020.6-py36-osx-conda.yml
## 使い終わった yaml ファイルを消去する
rm qiime2-2020.6-py36-osx-conda.yml
```

仮想環境 ``qiime2-2020.6`` に入る。この仮想環境では [Qiime2](https://qiita.com/keisuke-ota/items/6399b2f2f7459cd9e418) を用いた解析ができる。

```terminal:ターミナル
conda activate qiime2-2020.6
```

仮想環境から抜けるには ``conda deactivate`` と入力する。

<details><summary> **仮想環境の容量を節約して standalone 版のみをインストールしたい場合** </summary><div>

必要最低限のパッケージのみが入っている ``python=3.6`` の仮想環境を用意してもよい。

```terminal:ターミナル
conda create -n mmvec python=3.6
conda activate mmvec
```

</div></details>

## TensorFlow のインストール

mmvec は **TensorFlow** を用いて細菌と代謝物の共存確率を推定する。ただし最新版の TensorFlow には対応しておらず、[**TensorFlow v1.14**](https://github.com/tensorflow/tensorflow/releases/tag/v1.14.0) を指定してインストールする必要がある。

```terminal:ターミナル
conda install -c conda-forge tensorflow=1.14
```

## NumPy のダウングレード

**TensorFlow v1.14** から最新の NumPy をインポートすると ``FutureWarning`` が発生する。これは NumPy のバージョンを 1.16.4 にダウングレードすると回避できる。警告が表示されても解析可能であるため、そのまま先に進んでもよい。

```terminal:ターミナル
conda install -c conda-forge numpy=1.16.4
```

<details><summary> **必要最低限のパッケージで仮想環境を構築した場合** </summary><div>

``conda`` ではなく ``pip3`` でダウングレードする。 

```terminal:ターミナル
pip3 install numpy==1.16.4
```

</div></details>

## mmvec のインストール

以上でパッケージ依存関係が整備されたので、mmvec をインストールする。

```terminal:ターミナル
conda install mmvec -c conda-forge
```

これで standalone 版と Qiime2 版の　mmvec が同時にインストールされる。**どちらも ``paired-omics`` コマンドで解析することができる。**それぞれの ``--help`` を確認しておこう。

<details><summary> **standalone 版の ``paired-omics``** </summary><div>

解析パラメータのデフォルト値は、[公式 script](https://github.com/biocore/mmvec/blob/master/scripts/mmvec) から調べることができる。

```terminal:ターミナル
$ mmvec paired-omics --help
Usage: mmvec paired-omics [OPTIONS]

Options:
  --microbe-file TEXT             Input microbial abundances
  --metabolite-file TEXT          Input metabolite abundances
  --metadata-file TEXT            Input sample metadata file
  --training-column TEXT          Column in the sample metadata specifying
                                  which samples are for training and testing.

  --num-testing-examples INTEGER  Number of samples to randomly select for
                                  testing

  --min-feature-count INTEGER     Minimum number of samples a microbe needs to
                                  be observed in order to not filter out

  --epochs INTEGER                Number of epochs to train
  --batch-size INTEGER            Size of mini-batch
  --latent-dim INTEGER            Dimensionality of shared latent space. This
                                  is analogous to the number of PC axes.

  --input-prior FLOAT             Width of normal prior for input embedding.
                                  Smaller values will regularize parameters
                                  towards zero. Values must be greater than 0.

  --output-prior FLOAT            Width of normal prior for input embedding.
                                  Smaller values will regularize parameters
                                  towards zero. Values must be greater than 0.

  --arm-the-gpu                   Enables GPU support
  --learning-rate FLOAT           Gradient descent decay rate.
  --beta1 FLOAT                   Gradient decay rate for first Adam momentum
                                  estimates

  --beta2 FLOAT                   Gradient decay rate for second Adam momentum
                                  estimates

  --clipnorm FLOAT                Gradient clipping size.
  --checkpoint-interval INTEGER   Number of seconds before a storing a
                                  summary.

  --summary-interval INTEGER      Number of seconds before a storing a
                                  summary.

  --summary-dir TEXT              Summary directory to save cross validation
                                  results.

  --embeddings-file TEXT          Path to save the embeddings learned from the
                                  model. If this is not specified, then this
                                  will be saved under `--summary-dir`.

  --ranks-file TEXT               Path to save the ranks learned from the
                                  model. If this is not specified, then this
                                  will be saved under `--summary-dir`.

  --ordination-file TEXT          Path to save the ordination learned from the
                                  model. If this is not specified, then this
                                  will be saved under `--summary-dir`.

  --help                          Show this message and exit.
```

</div></details>


<details><summary> **Qiime2 版の ``paired-omics``** </summary><div>

Qiime2 の入力は ``Inputs:``（入力成果物名）、``Parameters:``（解析時のパラメータ）、``Outputs:``（出力成果物名）に区別される。この中で ``[required]`` は必須事項であり、それ以外のパラメータは省略すると ``[default:]`` の値で解析される。

```terminal:ターミナル
$ qiime mmvec paired-omics --help
Usage: qiime mmvec paired-omics [OPTIONS]

  Performs bi-loglinear multinomial regression and calculates the
  conditional probability ranks of metabolite co-occurence given the microbe
  presence.

Inputs:
  --i-microbes ARTIFACT FeatureTable[Frequency]
                          Input table of microbial counts.          [required]
  --i-metabolites ARTIFACT FeatureTable[Frequency]
                          Input table of metabolite intensities.    [required]
Parameters:
  --m-metadata-file METADATA...
    (multiple arguments   Sample metadata table with covariates of interest.
     will be merged)                                                [optional]
  --p-training-column TEXT
                          The metadata column specifying which samples are
                          for training/testing. Entries must be marked `Train`
                          for training examples and `Test` for testing
                          examples.                                 [optional]
  --p-num-testing-examples INTEGER
                          The number of random examples to select if
                          `training-column` isn't specified.      [default: 5]
  --p-min-feature-count INTEGER
                                                                 [default: 10]
  --p-epochs INTEGER      The total number of iterations over the entire
                          dataset.                              [default: 100]
  --p-batch-size INTEGER  The number of samples to be evaluated per training
                          iteration.                             [default: 50]
  --p-latent-dim INTEGER                                          [default: 3]
  --p-input-prior NUMBER  Width of normal prior for the microbial
                          coefficients. Smaller values will regularize
                          parameters towards zero. Values must be greater than
                          0.                                      [default: 1]
  --p-output-prior NUMBER Width of normal prior for the metabolite
                          coefficients. Smaller values will regularize
                          parameters towards zero. Values must be greater than
                          0.                                      [default: 1]
  --p-learning-rate NUMBER
                          Gradient descent decay rate.        [default: 1e-05]
  --p-summary-interval INTEGER
                                                                 [default: 60]
Outputs:
  --o-conditionals ARTIFACT FeatureData[Conditional]
                          Mean-centered Conditional log-probabilities.
                                                                    [required]
  --o-conditional-biplot ARTIFACT PCoAResults % Properties('biplot')
                          Biplot of microbe-metabolite vectors.     [required]
Miscellaneous:
  --output-dir PATH       Output unspecified results to a directory
  --verbose / --quiet     Display verbose output to stdout and/or stderr
                          during execution of this action. Or silence output
                          if execution is successful (silence is golden).
  --examples              Show usage examples and exit.
  --citations             Show citations and exit.
  --help                  Show this message and exit.

```

</div></details>

**↑↑エポック数（epochs）**や**学習率（learning-rate）**といった機械学習において非常に重要なパラメータのデフォルト値が standalone 版と Qiime2 版で異なることに注意しよう。

# mmvec paired-omics コマンドの入力について

上述のとおり、standalone 版も Qiime2 版も ``paired-omics`` コマンドで mmvec を実行する。主な入力事項について以下にまとめた。再現性を得るためにはどれも重要なパラメータである。

| 入力情報（＊は必須）| standalone 版| Qiime2 版| 注釈 |
| :--- | :--- | :--- | :--- |
| Microbiome データファイル（＊）| --microbe-file 【tsv か biom ファイル】| --i-microbes 【qza ファイル】| qza は Qiime2 専用のデータ形式。biom データからインポート可能 |
| Metabolome データファイル（＊）| --metabolite-file 【tsv か biom ファイル】| --i-metabolites 【qza ファイル】| サンプル名を Microbiome データと揃えること|
| メタデータファイル| --metadata-file 【tsv ファイル】| --m-metadata-file 【tsv ファイル】| サンプルを学習データとテストデータに分ける|
| メタデータのカラム名| --training-column 【テキスト】|--p-training-column 【テキスト】| メタデータファイルがなければ必要ない|
| テストサンプルの数 | --num-testing-examples 【整数】|--p-num-testing-examples 【整数】| メタデータがない場合、学習データとテストデータはランダムに分けられる。原著論文ではサンプル数の 1 割弱をテストデータに用いている。|
| 細菌の最低 read 数| --min-feature-counts 【整数】| --p-min-feature-count 【整数】|Microbiome データは sparse なことが多い。ここで指定した read 数より少ない細菌を除外する|
| エポック数**（重要）**| --epochs 【整数】| --p-epochs 【整数】|学習データ １ サンプルあたりの学習回数。計算時間と比例するパラメータ。デフォルトでは小さめに設定されているため学習不足になりやすい。一方で大きく設定すると過学習が起こりうる。|
| バッチサイズ | --batch-size 【整数】| --p-batch-size 【整数】|学習データをいくつかに分けて学習させる際に、ひとまとまりにされる学習データの数|
| 潜在ベクトルの次元数| --latent-dim 【整数】| --p-latent-dim 【整数】|潜在ベクトルは主成分分析（PCA）の主成分のように扱える。可視化に有効|
| 1 つ目のアダムモーメント| --beta1 【1未満の小数】| なし |機械学習のハイパーパラメータ。モデルの精度などに影響|
| 2 つ目のアダムモーメント| --beta2 【1以下の小数】| なし |機械学習のハイパーパラメータ。モデルの精度などに影響|
| Microbiome の事前分布を定めるパラメータ| --input-prior 【0 以上の実数】| --p-input-prior  【0 以上の実数】|mmvec では共役事前分布を用いた [map 推定](https://qiita.com/jyori112/items/80b21422b15753a1c5a4)で共存確率を求める。ここでは micobiome の共役事前分布の分散を指定する|
| Metabolome の事前分布を定めるパラメータ| --output-prior 【0 以上の実数】| --p-output-prior  【0 以上の実数】|ここでは metabolome の共役事前分布の分散を指定する|
|GPU による高速化| --arm-the-gpu  【なし】|なし| GPU で高速に学習させる。Mac 非対応。|
|学習率 **（重要）**| --learning-rate 【0 以上の実数】| --p-learning-rate  【0 以上の実数】|大きいほど収束速度が早まる。しかし大きすぎると発散してしまい収束しない|
|出力フォルダ（＊）| --summary-dir 【テキスト】| --output-dir 【テキスト】|同じ名前のフォルダがあると自動で上書きされるため注意|
|学習経過を何秒おきにメモするか **（重要）**|--summary-interval 【整数】|なし|**TensorBoard** で学習経過を確認するために必要。この値が大きいと、学習経過がメモされず学習不足や過学習を確認できない。|

# standalone 版について

## Microbiome データと Metabolome データの形式

Microbiome データと Metabolome データは次のような tsv で用意する。両ファイルでサンプル名を揃えること。

```:Microbiome.tsv
	サンプルa	サンプルb	サンプルc
細菌A	5	10	12
細菌B	3	4	7
細菌C	0	1	2
```

```:Metabolome.tsv
	サンプルa	サンプルb	サンプルc
代謝物A	5	10	12
代謝物B	3	4	7
代謝物C	0	1	2
```

## メタデータの形式

学習データ（Train）とテストデータ（Test）を以下の tsv ファイルで指定する。

```:Meta.tsv
サンプル名	Testing
サンプルa	Test
サンプルb	Train
サンプルc	Train
サンプルd	Train
```

## 学習経過の確認

mmvec を実行したディレクトリで ``tensorboard --logdir .`` を実行すると TensorBoard で学習経過を確認できる。TensorBoard で確認可能な学習経過は精度評価指標 **rmse** と損失関数 **logloss** である。mmvec 実行中は ``--summary-interval`` で指定した秒数ごとに学習経過が計算されるので、この値を小さく設定するほどグラフは滑らかになる。**rmse** や **logloss** が下がり続けている間は学習途中といえる。**rmse** や **logloss** が収束していない場合はエポック数を増やす必要がある。

![スクリーンショット 2021-07-21 12.45.18.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/612858/4e91c3e1-8a5e-b27c-404c-1092eb6ce67b.png)

　上図の場合、**rmse** や **logloss** はほぼ収束しているためエポック数は適切である可能性が高い。

# Qiime2 版について
## ``qza`` ファイルの作成

Qiime2 版の mmvec を使用したい場合、 biom データからインポートした ``qza`` ファイルが必要である。手元に tsv ファイルしかない場合、まずは ``biom convert`` をつかって biom ファイルを得る。``biom convert`` は仮想環境構築時に自動でインストールされる。biom ファイルの形式については[ biom 公式ドキュメント](https://biom-format.org/documentation/format_versions/biom-2.0.html)も参考にするとよい。

```terminal:ターミナル
biom convert -i microbiome.tsv -o microbiome.biom --table-type="OTU table" --to-hdf5
biom convert -i metabolome.tsv -o metabolome.biom --table-type="Metabolite table" --to-hdf5
```

biom ファイルから ``qza`` ファイルへのインポートは以下のように行う。[Qiime2 公式ドキュメント](https://docs.qiime2.org/2020.6/tutorials/importing/#feature-table-data)も参考にするとよい。

```terminal:ターミナル
qiime tools import \
        --input-path microbiome.biom \
        --output-path microbiome.qza \
        --type FeatureTable[Frequency]

qiime tools import \
        --input-path metabolome.biom \
        --output-path metabolome.qza \
        --type FeatureTable[Frequency]
```

# standalone 版と Qiime2 版の比較

**特別な理由がない限り standalone 版を勧める。**
（特別な理由の例： 菌叢データを qza ファイルで既に持っている。Qiime2 の操作に慣れている。）

| | standalone 版| Qiime2 版| 
| :--- | :--- | :--- | 
| メリット|tsv で入力可能。乱数が固定されておりパラメータが同じであれば同じ解析結果が得られる|heatmap や emperior といった可視化機能が揃っている|
| デメリット|可視化などの付属機能はない|``qza`` ファイルのインポートが煩雑。一部パラメータを使えない。GPU 非対応。乱数が固定されていないため毎回解析結果が異なる|