QIIME 2 でアンプリコンシーケンスを解析

アンプリコンシーケンスのデータ解析部分をバイオインフォマティクスプラットフォームの [QIIME 2](https://qiime2.org) を使ってやっていくにあたり、どうコマンドラインを打ち込んでいくかという備忘録です。

今回扱うデータとしては、イルミナ社の MiSeq による 16S rRNA 菌叢解析のプロトコールに基づいたアンプリコンシーケンス結果になります。

解析処理をすべて行うと次の視覚化データファイルが得られます。

* `demux.qzv` - シーケンスとクオリティの概要
* `table.qzv` - フィーチャー（OTU）
* `rep-seqs.qzv` - 代表配列
* `jaccard_emperor.qzv` - Jaccard 指数（類似度距離）
* `bray_curtis_emperor.qzv` - Bray-Curtis 指数（類似度距離）
* `unweighted_unifrac_emperor.qzv` - Unweighted UniFrac （質的距離）
* `weighted_unifrac_emperor.qzv` - Weighted UniFrac （量的距離）
* `faith-pd-group-significance.qzv` - 試料グループ間の系統学的多様性
* `evenness-group-significance.qzv` - 試料グループ間の均等度
* `unweighted-unifrac-body-site-significance.qzv` - 試料グループ間の類似度検定
* `taxonomy.qzv` - フィーチャーの系統分類
* `taxa-bar-plots.qzv` - 試料の系統構成比率

## 入力ファイルとその概要の生成

シーケンスはペアエンドリードシーケンスになりまして、その結果データファイルは複数あるのですが、 [QIIME 2 でインポート処理](https://qiita.com/HiguchiRyusuke/items/90ca05c60ab8210c2d7a#入力データの作成)をすることで逆多重化（demultiplexed）され、ひとつファイル（アーティファクト）となります。ここで得られたファイル名は `demux.qza` とします。

このデータファイルは、QIIME 2 におけるセマンティックタイプが `SampleData[PairedEndSequencesWithQuality]` となっています。これは公式チュートリアルで扱っているアーティファクトのセマンティックタイプとは異なるもので、使用するコマンドも一部で異なっています。

入力ファイルを得たら解析の始まりです。まずはデータファイルの概要を「視覚化」します。次のようなコマンドになります。

```console
qiime demux summarize \
  --i-data demux.qza \
  --o-visualization demux.qzv
```

`demux.qzv` は視覚化ツールから見ることができます。[ウェブ版ビューワー](https://view.qiime2.org/)から見るのが手軽です。

`demux.qzv` をビューワーで開くと、「Interactive Quality Plot」タブからクオリティスコア分布のグラフが見られるはずです。縦軸にクオリティスコア、横軸に配列長でプロットされています。ペアエンドリード形式では図が 2 つ表示されます。左がフォワードのリードの、右がリバースのリードのものになります。

これらの分布を元に次の品質管理（クオリティコントロール、ノイズ除去）工程のパラメーターを決定します。

## 品質管理（クオリティコントロール）処理

品質管理処理を行うことで、データセット中の精度の高いものだけを採用し、解析に適した「代表配列（フィーチャー）」を獲得することができます。QIIME 2 では品質管理処理の結果ファイルとして**フィーチャーテーブル**（`FeatureTable[Frequency]` アーティファクト）と、その**代表配列**（`FeatureData[Sequence]` アーティファクト）が得られます。「フィーチャー」 は QIIME 1 での **OTU** または **BIOM** に相当するそうです。

QIIME 2 の品質管理処理にはいくつかアルゴリズムが用意されています。いずれも QIIME 1 で標準となっている「97% OTU」の結果よりもいい品質があるそうなので、安心して導入してよさそうです。

### DADA2 によるノイズ除去

今回のペアエンドリードシーケンスに対してのノイズ除去は DADA2 による `dada2 denoise-paired` を利用します。実行には 2 つパラメーターを、フォワード側（`f`）とリバース側（`r`）でそれぞれ指定します。

ひとつめは「採用する最大配列長」となる `--p-trunc-len-(r|f) <length>` オプションです。`<length>` 塩基以降の配列を除去（truncate）します。

もうひとつは「先頭側のカット長」となる `--p-trim-left-(r|f) <position>` オプションです。各シーケンス結果の最初の `<position>` 塩基を削除（trim）します。

採用する最大配列長には結合したペアエンドリードが十分重なる長さを指定することが求められます。V3-V4 領域を含むシーケンスを例とすると、足し合わせて 460bp を超えるよう指定します。それぞれの値は先ほどのクオリティスコア分布をみて決定します。両端が同程度にクオリティスコアを維持している長さを採用すると、フィーチャーが得やすいでしょう。

プライマー配列は取り除く必要があるので、先頭側のカット長はプライマー配列長を通常指定します。

たとえば次のようなコマンドになります。

```console
$ qiime dada2 denoise-paired \
    --i-demultiplexed-seqs demux.qza \
    --o-table table.qza \
    --o-representative-sequences rep-seqs.qza \
    --p-trim-left-f 17 \
    --p-trim-left-r 21 \
    --p-trunc-len-f 280 \
    --p-trunc-len-r 220
```

これにより、フィーチャーテーブル `table.qza` とその代表配列の詳細情報が `rep-seqs.qza` へと出力されます。

#### トラブルシューティング

`dada2 denoise-paired` 時に次のようなエラーで止まる場合があります。

```
Plugin error from dada2:

  An error was encountered while running DADA2 in R (return code 1), please inspect stdout and stderr to learn more.

Debug info has been saved to /tmp/qiime2-q2cli-err-jfu9_pbl.log
```

表示されているパスにあるログをみるか、`--verbose` オプションをつけてもう一度実行し、出力されたエラー内容を見て、解決を図っていきます。

```log
4) Remove chimeras (method = consensus)
Error in isBimeraDenovoTable(unqs[[i]], ..., verbose = verbose) : 
  Input must be a valid sequence table.
Calls: removeBimeraDenovo -> isBimeraDenovoTable
```

もし上記の内容のエラーとなった場合、原因にはペアエンドリードが重ならなかったことがあげられます。この場合にはパラメーターの `--p-trunc-len-(r|f)` を、足し合わせた長さが増幅領域を超える長さまで十分に伸ばしてあげることで解決されるかもしれません。（参考 - https://forum.qiime2.org/t/qiime-dada2-denoise-paired-end-sequences/1062/9 ）

### フィーチャーテーブルの概要を視覚化

得られたフィーチャー情報の視覚化データを作成することで、品質管理処理の結果を確認できます。2 つのファイルそれぞれにコマンドを実行します。

`feature-table summarize` コマンドではフィーチャーテーブル（`table.qza`）の視覚化データを作成します。たとえば次のようにコマンドを実行します。

```console
$ qiime feature-table summarize \
  --i-table table.qza \
  --o-visualization table.qzv
```

`table.qzv` をビューワーで見ると、品質管理によってどれくらいの配列（フィーチャー）をサンプルから得られたのかが確認できます。あまり配列が得られていないようでしたら、品質管理のパラメーターを変えたやり直しを行うとよいでしょう。

`feature-table tabulate-seqs` コマンドでは代表配列（`rep-seqs.qza`）の視覚化データを作成します。たとえば次のようにコマンドを実行します。

```console
$ qiime feature-table tabulate-seqs \
  --i-data rep-seqs.qza \
  --o-visualization rep-seqs.qzv
```

`rep-seqs.qzv` ではフィーチャーの代表配列を見ることができます。2 つのファイルは「フィーチャー ID」でマッピングされており、各代表配列には BLAST へのリンクがつけられています。


## 多様性解析

得られたフィーチャー情報によって系統学的多様性解析が行えます。次の手順で行います。

* 多重配列整列
* 系統樹の生成
* メタデータの作成（無ければ）
* 多様性指標の生成と視覚化

一連の解析処理によってα多様性とβ多様性の指標データが得られます。多様性の説明をすると、α多様性とは試料から取り出された代表配列の多様性であり、β多様性とは試料間の代表配列集合の不一致性であることになります。

### 多重配列整列（マルチプルアライメント）

多重配列整列は MAFFT を利用しますが、こちらも QIIME 2 から実行できます。`qiime alignment mafft` コマンドを利用します。

```console
$ qiime alignment mafft \
  --i-sequences rep-seqs.qza \
  --o-alignment aligned-rep-seqs.qza
```

続いてこの結果に対してマスクを行います。マスクによって系統学的にノイズと思われる高ギャップ部分を抽出（除去）してくれるそうです。

```console
$ qiime alignment mask \
  --i-alignment aligned-rep-seqs.qza \
  --o-masked-alignment masked-aligned-rep-seqs.qza
```

### 系統樹の生成

多様性解析といえば系統樹です。マスク済み多重配列整列を使って系統樹を作成します。アルゴリズムとしては FastTree を使います。次のようにコマンドを実行します。

```console
$ qiime phylogeny fasttree \
  --i-alignment masked-aligned-rep-seqs.qza \
  --o-tree unrooted-tree.qza
```

これで得られるのは無根系統樹です。有根系統樹を得るにはさらに「根」を決定する処理を行います。

```console
$ qiime phylogeny midpoint-root \
  --i-tree unrooted-tree.qza \
  --o-rooted-tree rooted-tree.qza
```

以上で多様性解析のもとになるデータの生成は完了です。

### メタデータの作成

試料間の分析を行うにあたって、メタデータ（試料についての情報）ファイルを作っておく必要があります。メタデータが詳細であるほど、視覚化データの観覧時に使用できるオプションが増えることになります。

メタデータファイルはタブ区切りテキストで作ります。最初の列に「試料 ID」を入れることだけが必須です。あとは自由に列を作れるので、グループごとに解析をしたいときはその属性名で列を作るとよいでしょう。ここでは「`Group`」という列を作ったものとします。[QIIME 1 のときと同じ記述方法](http://qiime.org/documentation/file_formats.html#metadata-mapping-files)でも大丈夫だそうです。「試料 ID」は `table.qzv` の「Sample ID」と一致しているかを確認するとよいでしょう。

メタデータを視覚化データにしてビューワーで見ることもできます。ここでは `sample-metadata.tsv` という名前で作ったものとします。視覚化データにしてビューワーから確認することができます。次のようなコマンドを実行します。

```console
$ qiime metadata tabulate \
  --m-input-file sample-metadata.tsv \
  --o-visualization tabulated-sample-metadata.qzv
```

### 系統的指標の生成と視覚化

ここまでで生成してきたファイルたちに統計学的処理を与えて、多様性解析の主要な指標と視覚化データを作成します。`diversity core-metrics-phylogenetic`
コマンドで行いますが、ここでは `--p-sampling-depth` という重要なパラメーター指定があります。

多様性の計算は各試料（場所・個体）内から得られた代表配列の数をもとにします。代表配列の数は元々の試料のリード数（濃度）に依存している可能性があり、リード数が少なければ得られるフィーチャーも少なくなってしまうということです。

ですのでサンプリング濃度をリード数の低い試料にあわせる「希薄化（レアファクション）」計算を行うことで、その偏りを補正しようという考えだそうです。「希薄化」については http://pedsurgery.wp.xdomain.jp/?p=213 で詳しく解説されていまして、ここの説明にも参考にさせていただきました。

さて `--p-sampling-depth` の求め方ですが、ビューワーで `table.qzv` の「Interactive Sample Detail」から行うとよいそうです。すべての試料が収まる最大の値が適切です。希薄化しすぎると元々の結果から離れてしまうことになりますので、サンプルを切り捨てて解析することも場合によっては適切と考えられます。ほかの判断方法として、今回は省略しましたが希薄化度合いのデータ `alpha-rarefaction.qzv` が生成してあれば、それを見て決めてもいいでしょう。

たとえば次のようなコマンドで実行します。

```console
$ qiime diversity core-metrics-phylogenetic \
  --i-phylogeny rooted-tree.qza \
  --i-table table.qza \
  --m-metadata-file sample-metadata.tsv \
  --p-sampling-depth 841 \
  --output-dir core-metrics-results
```

UniFrac 距離をビューワーで 3 次元で確認できる視覚化データなどが生成されます。

#### 試料をグループ化した多様性解析

`sample-metadata.tsv` に記述した情報から、試料グループの有意性解析を行うことができます。

簡単にコマンド例だけ紹介します。α多様性については次の通りです。

```console
$ qiime diversity alpha-group-significance \
  --i-alpha-diversity core-metrics-results/faith_pd_vector.qza \
  --m-metadata-file sample-metadata.tsv \
  --o-visualization core-metrics-results/faith-pd-group-significance.qzv
```

```console
$ qiime diversity alpha-group-significance \
  --i-alpha-diversity core-metrics-results/evenness_vector.qza \
  --m-metadata-file sample-metadata.tsv \
  --o-visualization core-metrics-results/evenness-group-significance.qzv
```

β多様性については次の通りです。

```console
$ qiime diversity beta-group-significance \
  --i-distance-matrix core-metrics-results/unweighted_unifrac_distance_matrix.qza \
  --m-metadata-file sample-metadata.tsv \
  --m-metadata-category Group \
  --o-visualization core-metrics-results/unweighted-unifrac-body-site-significance.qzv \
  --p-pairwise
```

 
## 系統分類の推定

今回は 16S rRNA の菌叢解析ですが、どの系統で構成されているかを図表化することができます。参照するマーカー遺伝子データベースはあらかじめ増幅した領域について鍛えておくことでよりよい結果が得られるそうです。

マーカー遺伝子データベースは、[オフィシャルで用意されている](https://docs.qiime2.org/2017.12/data-resources/)ものをダウンロードして利用できます。「鍛えられた」データーベースは 515F/806R 領域のものしか用意されていないため、今回は Greengenes の汎用的な全長のもの `gg-13-8-99-nb-classifier.qza` を利用します。

系統分類処理は次のように行います。

```console
$ qiime feature-classifier classify-sklearn \
  --i-classifier gg-13-8-99-nb-classifier.qza \
  --i-reads rep-seqs.qza \
  --o-classification taxonomy.qza
```

これを視覚化データにするにはさらに次のコマンドを実行します。

```console
$ qiime metadata tabulate \
  --m-input-file taxonomy.qza \
  --o-visualization taxonomy.qzv
```

ここで得られた生物分類について、精度が気になるときは Feature ID をもとに `table.qzv` につけられた BLAST リンクをたどってみるとよいでしょう。求めている分類精度によっては独自のデータベースの BLAST を利用したくなるかもしれません。ローカル BLAST による分類処理も QIIME 2 で用意されています。

最後に、得られた系統分類データを、はじめの方のフィーチャーテーブルと照らし合わせることで、試料ごとの生物系統の構成の視覚化データを作成できます。

```console
$ qiime taxa barplot \
  --i-table table.qza \
  --i-taxonomy taxonomy.qza \
  --m-metadata-file sample-metadata.tsv \
  --o-visualization taxa-bar-plots.qzv
```


## おわりに

バイオインフォマティクスとしてアンプリコンシーケンスを扱おうとすると、コマンドラインインターフェイスや Python （Ruby または Perl）スクリプトを頑張って書いていかなければならなかったものです。いまでは統合されたツールとして QIIME 2 があるおかげで、プログラマーでなくてもだいぶ扱いやすくなった印象です。プログラマーとしても複雑さから離れることができるので、ウェットな方面の知識が得やすくなりました。