初心者の菌叢解析 Qiime2で解析(10) 多様性解析 ~α多様性~

前回までの解析で菌叢の組成を可視化するところまでいきました。今回は多様性解析の方法について記載します。

1．多様性について

菌叢の多様性はその環境を考察する上で非常に重要です。腸内細菌叢では宿主の疾患や薬の有効性に多様性が関与していることが報告されています。

多様性にはα多様性とβ多様性がありますが、今回はα多様性について解析していきます。多様性については以下のサイトがわかりやすく示されていました。

α多様性とβ多様性 - Qiita

2．レアファクションカーブの作成

まずはレアファクションカーブを作成していきます。解析用ディレクトリ（Qiime2_test）に移動し、Qiime2をアクティベートします。

cd /Users/ユーザー名/Desktop/Qiime2_test 
conda activate qiime2-2021.2

以下のコマンドを実行し系統樹のアーティファクトを作成します。
作成には時間がかかりますので、ご注意ください。また、先日まではここの記載を失念していまして、2021年11月25日に追記いたしました。
無事解析が終了すると「aligned-rep-seqs.qza」「masked-aligned-rep-seqs.qza」「unrooted-tree.qza」「rooted-tree.qza」の4つのファイルが生成されます。

qiime phylogeny align-to-tree-mafft-fasttree \
--i-sequences rep-seqs.qza \
--o-alignment aligned-rep-seqs.qza \
--o-masked-alignment masked-aligned-rep-seqs.qza \
--o-tree unrooted-tree.qza \
--o-rooted-tree rooted-tree.qza

以下のコマンドを実行し、Qiime2_testのフォルダの中に「alpha-diversity」というフォルダを作成し、その中に「alpha-rarefaction」というフォルダを作成します。

mkdir alpha-diversity
mkdir alpha-diversity/alpha-rarefaction

以下のコマンドでレアファクションカーブを生成します。
今回は「observed features」「shannon」「faith pd」という3つの指標で解析します。

レアファクションカーブではサンプリングの回数（sample depth）が多い方がレアな菌を抽出できます。ただ、サンプリングの回数は得られたRead数等に依存する為、その都度考える必要があります。
今回は10000回のサンプリングでレアファクションカーブを作成します。

qiime diversity alpha-rarefaction  \
--i-table table.qza  \
--i-phylogeny rooted-tree.qza  \
--p-metrics observed_features  \
--p-metrics shannon  \
--p-metrics faith_pd  \
--p-max-depth 10000  \
--m-metadata-file sample-metadata.txt  \
--o-visualization alpha-diversity/alpha-rarefaction/sample_depth-10000.qzv

「--p-max-depth 10000」の値を変更すればsample depthを変更できます。
うまくいけば「sample_depth-10000.qzv」ファイルが「alpha-rarefaction」フォルダに生成されていると思います。

Qiime2 Viewで確認すると以下の様になりました。

赤枠の「Metric」の部分をクリックすれば他の指標を見ることが出来ます。また、左のボタンからCSV形式で保存も可能です。

今回の結果からは大阪湾の方がα多様性が高いことが読み取れます。さらに、マウスの糞便データではsample depthが1000程度で横ばいになっており、1000程度のサンプリング回数で、マウスのα多様性（faith_pd）は十分解析可能ということになります。もちろん指標が変われば横ばいになるsample depthも変わりますので、ご注意ください。

3．α多様性指数の算出(2022年3月4日変更）

レアファクションカーブを参考にし、α多様性指数の算出を行っていきます。
今回はSample depthを10,000で解析してみます。
Shannon indexはSample depthが1000程度で頭打ちになりましたので、
このような場合は1000程度のSample depthでの解析で問題ありません。
解析したい指数に応じてレアファクションカーブで確認し、算出を行うことが良いと思います。

3ｰ1．多様性指数の算出

それではまずはディレクトリを作成していきます。
「alpha-diversity」の中に、「Sample_depth_10000」を作成します。
次にその中に、「alpha-diversity-index」と「alpha-group-significance」を作成していきます。

mkdir alpha-diversity/Sample_depth_10000
mkdir alpha-diversity/Sample_depth_10000/alpha-diversity-index alpha-diversity/Sample_depth_10000/alpha-group-significance

次にSample depthが10000の時のα多様性指数の算出を行っていきます。
以下のコマンドを実行してください。

qiime diversity core-metrics-phylogenetic \
--i-phylogeny rooted-tree.qza \
--i-table table.qza \
--p-sampling-depth 10000 \
--m-metadata-file sample-metadata.txt \
--output-dir alpha-diversity/Sample_depth_10000/core-metrics-results

こちらのコマンドは次回の投稿で解析するβ多様性解析でも使用します。
「--p-sampling-depth 10000」の部分がSample depthの設定になります。
ここの値を変化させることで、そのSample depthに合わせた解析が可能になります。
生成される「core-metrics-results」では多くのファイルが生成されますが、今回はα多様性に絞って解析を進めていきます。

3-2. 多様性指数の可視化

ここでは以下の4つのα多様性について解析が可能になります。
「evenness」「faith_pd」「observed_features」「shannon」

次にそれぞれの指数についてqzvファイルを作成し、可視化していきます。
以下のコマンドを実行します。

qiime metadata tabulate \
--m-input-file alpha-diversity/Sample_depth_10000/core-metrics-results/evenness_vector.qza \
--o-visualization alpha-diversity/Sample_depth_10000/alpha-diversity-index/evenness_index

qiime metadata tabulate \
--m-input-file alpha-diversity/Sample_depth_10000/core-metrics-results/faith_pd_vector.qza \
--o-visualization alpha-diversity/Sample_depth_10000/alpha-diversity-index/faith_pd_index

qiime metadata tabulate \
--m-input-file alpha-diversity/Sample_depth_10000/core-metrics-results/observed_features_vector.qza \
--o-visualization alpha-diversity/Sample_depth_10000/alpha-diversity-index/observed_features_index

qiime metadata tabulate \
--m-input-file alpha-diversity/Sample_depth_10000/core-metrics-results/shannon_vector.qza \
--o-visualization alpha-diversity/Sample_depth_10000/alpha-diversity-index/shannon_index

qzvファイルが「alpha-diversity-index」ディレクトリの中に生成されたと思います。
「faith_pd」の値をQiime2 viewで確認してみます。
同時に生成した他の指数につても同様に確認できます。

Sample depth 10000における「faith pd」指数の結果

大阪湾の菌叢の方が「faith_pd」の値が高いことが分かります。
実際に有意な変化なのか統計的に調べてみます。

3ｰ3．多様性指数の統計処理

それでは統計的に検定を行ってみます。
以下のコマンドを実行してください。

qiime diversity alpha-group-significance \
--i-alpha-diversity alpha-diversity/Sample_depth_10000/core-metrics-results/evenness_vector.qza \
--m-metadata-file sample-metadata.txt \
--o-visualization alpha-diversity/Sample_depth_10000/alpha-group-significance/evenness_significance.qzv

qiime diversity alpha-group-significance \
--i-alpha-diversity alpha-diversity/Sample_depth_10000/core-metrics-results/faith_pd_vector.qza \
--m-metadata-file sample-metadata.txt \
--o-visualization alpha-diversity/Sample_depth_10000/alpha-group-significance/faith_pd_significance.qzv

qiime diversity alpha-group-significance \
--i-alpha-diversity alpha-diversity/Sample_depth_10000/core-metrics-results/observed_features_vector.qza \
--m-metadata-file sample-metadata.txt \
--o-visualization alpha-diversity/Sample_depth_10000/alpha-group-significance/observed_features_significance.qzv

qiime diversity alpha-group-significance \
--i-alpha-diversity alpha-diversity/Sample_depth_10000/core-metrics-results/shannon_vector.qza \
--m-metadata-file sample-metadata.txt \
--o-visualization alpha-diversity/Sample_depth_10000/alpha-group-significance/shannon_significance.qzv

「alpha-group-significance」の中に各指数の統計解析結果が出力されます。
生成された「faith_pd_significance.qzv」をQiime2 Viewで確認してみます。

Sample depth 10000における「faith pd」指数の統計結果

統計処理はKruskal-Wallis検定で行っています。
P値を確認すると0.05以下ということで、有意差ありと判定できそうです。
多重比較を行う場合は下にスクロースすると「Kruskal-Wallis (pairwise)」がありますので、そこのQ値を確認するのが良いと思います。

マウス腸内細菌叢と大阪湾ではそもそもが環境が違うので、α多様性を直接比較することはあまり意味がないとも思いますが、大阪湾の方がfaith-pdが有意に高いことが分かります。

データをダウンロードも出来ますので、確認してみてください。

読んでいただきありがとうございました。
次回はβ多様性について書きます。