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

今回扱うデータとしては、イルミナ社の 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 でインポート処理をすることで逆多重化（demultiplexed）され、ひとつファイル（アーティファクト）となります。ここで得られたファイル名は demux.qza とします。

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

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

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

demux.qzv は視覚化ツールから見ることができます。ウェブ版ビューワーから見るのが手軽です。

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 を超えるよう指定します。それぞれの値は先ほどのクオリティスコア分布をみて決定します。両端が同程度にクオリティスコアを維持している長さを採用すると、フィーチャーが得やすいでしょう。

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

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

$ 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 オプションをつけてもう一度実行し、出力されたエラー内容を見て、解決を図っていきます。

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）の視覚化データを作成します。たとえば次のようにコマンドを実行します。

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

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

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

$ 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 コマンドを利用します。

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

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

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

系統樹の生成

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

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

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

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

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

メタデータの作成

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

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

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

$ 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 が生成してあれば、それを見て決めてもいいでしょう。

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

$ 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 に記述した情報から、試料グループの有意性解析を行うことができます。

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

$ 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

$ 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

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

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

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

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

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

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

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

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

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

$ 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 があるおかげで、プログラマーでなくてもだいぶ扱いやすくなった印象です。プログラマーとしても複雑さから離れることができるので、ウェットな方面の知識が得やすくなりました。