メタゲノム解析
2019/9/17
DRY解析教本の「メタゲノム解析」をやってみる。
ikraの三島合宿にて。
参考：微生物群集解析ツール QIIME 2 を使う．Part1.（インストール，ファイルインポート編）（伊藤さん）
install
https://gyazo.com/0c964f7a7dbe67845fbefa0316563cc5
・qiime2（チャイム2）
Qiime2 （チャイム2）は、次世代シークエンサーによって得られた16S rRNA可変部位領域等の大量配列情報をもとにして微生物叢のプロファイリング解析を行う際に便利なパッケージングツール
condaで。
code:bash
$ conda install -c qiime2/label/r2019.1 qiime2
versionがいろいろあるが、本と揃える。
https://gyazo.com/6e980f9eba83e781f993d025ef79ba4a
＜エラー＞
code:bash
・qiime2のインストール方法
https://docs.qiime2.org/2019.7/install/
code:bash
Natively installing QIIME 2
Install Miniconda
Install QIIME 2 within a conda environment
Activate the conda environment
Installing QIIME 2 using Virtual Machines
Installing QIIME 2 using VirtualBox
Installing QIIME 2 using Amazon Web Services
Installing QIIME 2 using Docker
公式ではminicondaが推奨されている。
仮想環境を作ってinstall（ymlを使う）
https://gyazo.com/d5f1be03276da88783065ab7eef58882
＜エラー＞
code:bash
Collecting package metadata (repodata.json): done
Solving environment: failed
ResolvePackageNotFound:
- libgfortran=3.0.1
- llvm=4.0.1
- clangxx_osx-64=4.0.1
- libcxxabi=4.0.1
- gfortran_osx-64=4.8.5
- libcxx=4.0.1
- clang=4.0.1
- compiler-rt=4.0.1
- clang_osx-64=4.0.1
- ld64=274.2
- llvm-lto-tapi=4.0.1
- appnope=0.1.0
- cctools=895
- clangxx=4.0.1
→間違えて、Mac用でinstallしていた。。。
titanでは、linux用を選ぶ。
code:bash
$ wget https://data.qiime2.org/distro/core/qiime2-2019.7-py36-linux-conda.yml
$ conda env create -n qiime2-2019.7 --file qiime2-2019.7-py36-linux-conda.yml
# OPTIONAL CLEANUP
$ rm qiime2-2019.7-py36-linux-conda.yml
#activate
$ source activate qiime2-2019.7 #これは古い
$ conda activate qiime2-2019.7
$ qiime --help
無事動いた！
・rename
conda install rename
これはダメ。。。
（renameは入っていたが、別の標準のもの？）
↓
conda install -c bioconda rename
16S rRNA遺伝子リファレンスデータベースのダウンロード
qiime2のData resources
wget https://data.qiime2.org/2019.7/common/gg-13-8-99-nb-classifier.qza
code:bash
#確認
$ ls -lah gg-13-8-99-nb-classifier.qza
-rw-rw-r-- 1 kyamada kyamada 100M 6月 11 02:39 gg-13-8-99-nb-classifier.qza
ファイルのサイズの確認。
解析用データのダウンロード
＜データについて＞
MiSeqでシーケンスされたマウス糞便試料のデータ。
野生型と遺伝性肥満マウス(ob/obマウス)を用いて、乳酸菌 L. paracasei K71株摂取による腸内細菌叢への影響を解析。
SRA→PRJDB7113：リンク
https://gyazo.com/a65752a34dc5c98a52a56688b523428c
code:bash
#SRAを取得
$ prefetch --option-file SRR_Acc_List.txt
$ mkdir fastq
＊prefetchの保存先は、~/ncbi/public/sra/！
（SRA_Toolkit 使い方 データ保存フォルダの変更）
↓
https://gyazo.com/01d3883a018d4a5650cd2de2f5b1c88f
fastq-dump
fastq-dumpコマンドは一度に一つのRunIDしか引数として受け付けない。
catとxargsコマンドと組み合わせることでRunIDリストの値を一行ずつfastq-dumpへ渡す処理を繰り返し実行。
code:bash
$ time cat ../SRR_Acc_List.txt | xargs -n1 fastq-dump --gzip --split-files
Read 119773 spots for DRR138864
Written 119773 spots for DRR138864
Read 216616 spots for DRR138865
Written 216616 spots for DRR138865
Read 143745 spots for DRR138866
Written 143745 spots for DRR138866
Read 116416 spots for DRR138867
Written 116416 spots for DRR138867
Read 103126 spots for DRR138868
Written 103126 spots for DRR138868
Read 94873 spots for DRR138869
Written 94873 spots for DRR138869
Read 108061 spots for DRR138870
Written 108061 spots for DRR138870
Read 89155 spots for DRR138871
Written 89155 spots for DRR138871
Read 89115 spots for DRR138872
Written 89115 spots for DRR138872
Read 98527 spots for DRR138873
Written 98527 spots for DRR138873
Read 130627 spots for DRR138874
Written 130627 spots for DRR138874
Read 114420 spots for DRR138875
Written 114420 spots for DRR138875
Read 121840 spots for DRR138876
Written 121840 spots for DRR138876
Read 125897 spots for DRR138877
Written 125897 spots for DRR138877
Read 84334 spots for DRR138878
Written 84334 spots for DRR138878
Read 107963 spots for DRR138879
Written 107963 spots for DRR138879
Read 122359 spots for DRR138880
Written 122359 spots for DRR138880
Read 114831 spots for DRR138881
Written 114831 spots for DRR138881
Read 93020 spots for DRR138882
Written 93020 spots for DRR138882
Read 152211 spots for DRR138883
Written 152211 spots for DRR138883
Read 143313 spots for DRR138884
Written 143313 spots for DRR138884
Read 124645 spots for DRR138885
Written 124645 spots for DRR138885
Read 134223 spots for DRR138886
Written 134223 spots for DRR138886
real 3m56.600s
user 3m52.726s
sys 0m2.293s
確認
ls * | wc -l
46
ファイル名の変換
FASTQファイル名をMiSeq形式のSampleName_S1_L001_R1_001.fastq.gzのようなファイル名に変換
（Sはサンプル番号、Rはリード番号を意味）
renameコマンド。
code:bash
$ rename 's/_1/_S1_L001_R1_001/' *.fastq.gz
$ rename 's/_2/_S1_L001_R2_001/' *.fastq.gz
$ ls
メタデータの作成
サンプルの属性や群分けなどの情報を記述したメタデータを作成。
SraRunTable.txtを以下のように変換。
タブ区切りテキスト (.txt)、ファイル名をmetadata.txt
↓
code:bash
#SampleID host_diet host_genotype group
DRR138864 AIN ob/ob ob-AIN
DRR138865 AIN ob/ob ob-AIN
DRR138866 AIN ob/ob ob-AIN
DRR138867 AIN ob/ob ob-AIN
DRR138868 AIN ob/ob ob-AIN
DRR138869 AIN ob/ob ob-AIN
DRR138870 AIN ob/ob ob-AIN
DRR138871 K71 ob/ob ob-K71
DRR138872 K71 ob/ob ob-K71
DRR138873 K71 ob/ob ob-K71
DRR138874 K71 ob/ob ob-K71
DRR138875 K71 ob/ob ob-K71
DRR138876 AIN ob/ob WT-AIN
DRR138877 AIN ob/ob WT-AIN
DRR138878 AIN ob/ob WT-AIN
DRR138879 K71 wild-type ob-K71
DRR138880 K71 wild-type ob-K71
DRR138881 K71 wild-type ob-K71
DRR138882 AIN wild-type WT-AIN
DRR138883 AIN wild-type WT-AIN
DRR138884 AIN wild-type WT-AIN
DRR138885 AIN wild-type WT-AIN
DRR138886 AIN wild-type WT-AIN
【解析】
FASTQファイルのインポート
code:bash
$ time qiime tools import \
--type 'SampleDataPairedEndSequencesWithQuality' \
--input-path fastq \
--input-format CasavaOneEightSingleLanePerSampleDirFmt \
--output-path demux.qza
Imported fastq as CasavaOneEightSingleLanePerSampleDirFmt to demux.qza
real 1m11.349s
user 1m6.146s
sys 0m3.231s
https://gyazo.com/315166d0e1088688cab9c23fc345b442
サマリーファイルを作成
code:bash
$ time qiime demux summarize \
--i-data demux.qza \
--o-visualization demux.qzv
Saved Visualization to: demux.qzv
real 0m52.795s
user 0m49.545s
sys 0m2.295s
https://gyazo.com/d5dc0c5c8087741b805ebf1d82c1a1f5
結果ファイルの閲覧方法
QIIME 2から出力されるファイル
①拡張子が.qzaのファイル
②データを図表として可視化した.qzvファイル
の2種類。
＊.qzaファイルはQIIME 2でツール間の入出力ファイルとして使われユーザが直接扱う場面は少ない。
コマンドver
code:bash
$ qiime tools view ファイル名.qzv
（例）$ qiime tools view demux.qzv
ブラウザver
QIIME2 view
解析処理をすべて行うと次の視覚化データファイルが得られる。
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 - 試料の系統構成比率
インポート結果の確認
ブラウザにて。
demux.qzv
Overviewタブには、サンプル数、サンプルごとのリード数、統計値が表示。
https://gyazo.com/56b33179e8fbf11734684910df97ad2d
https://gyazo.com/86c11c9577ad5b62eaf13afe6922374b
23サンプル。
https://gyazo.com/5516e5fe158b98ba0baca48302d0de9f
トリミングの明確な基準は示されていないが、中央値が20未満になる塩基以降をトリミングするのが一般的
今回はForwardリードは最後の10塩基以降を、Reverseリードは最後の50塩基以降を次工程でトリミング。
Miseqだから300くらいまで読める。シーケンスは、先にいけばいくほどクオリティが悪くなる。
Hiseqは150程度。
シーケンスQCとFeature tableの構築
code:bash
#説明
$ qiime dada2 denoise-paired \
--verbose \
--p-n-threads 0 \（使用可能な全CPUコアを使う）
--p-trim-left-f 17 \（Forwardの5’側からトリミングする塩基数=PCRのプライマー長）
--p-trim-left-r 21 \（Reverseの5’側からトリミングする塩基数=PCRのプライマー長）
--p-trunc-len-f 290 \（Forwardの5’側から使用する塩基数、以降トリミングされる）
--p-trunc-len-r 250 \（Reverseの5’側から使用する塩基数、以降トリミングされる）
--i-demultiplexed-seqs demux.qza \
--o-table table.qza \
--o-representative-sequences rep-seqs.qza \
--o-denoising-stats stats-dada2.qza
#実行
$ time qiime dada2 denoise-paired \
--verbose \
--p-n-threads 0 \
--p-trim-left-f 17 \
--p-trim-left-r 21 \
--p-trunc-len-f 290 \
--p-trunc-len-r 250 \
--i-demultiplexed-seqs demux.qza \
--o-table table.qza \
--o-representative-sequences rep-seqs.qza \
--o-denoising-stats stats-dada2.qza
Running external command line application(s). This may print messages to stdout and/or stderr.
The command(s) being run are below. These commands cannot be manually re-run as they will depend on temporary files that no longer exist.
Command: run_dada_paired.R /tmp/tmpi44h9i02/forward /tmp/tmpi44h9i02/reverse /tmp/tmpi44h9i02/output.tsv.biom /tmp/tmpi44h9i02/track.tsv /tmp/tmpi44h9i02/filt_f /tmp/tmpi44h9i02/filt_r 290 250 17 21 2.0 2.0 2 consensus 1.0 0 1000000
R version 3.5.1 (2018-07-02)
要求されたパッケージ Rcpp をロード中です
DADA2: 1.10.0 / Rcpp: 1.0.2 / RcppParallel: 4.4.3
1) Filtering .......................
2) Learning Error Rates
284179896 total bases in 1040952 reads from 15 samples will be used for learning the error rates.
238378008 total bases in 1040952 reads from 15 samples will be used for learning the error rates.
3) Denoise remaining samples .......................
4) Remove chimeras (method = consensus)
6) Write output
Saved FeatureTableFrequency to: table.qza
Saved FeatureDataSequence to: rep-seqs.qza
Saved SampleDataDADA2Stats to: stats-dada2.qza
real 8m34.270s
user 175m48.642s
sys 8m18.681s
https://gyazo.com/48c990bd9394522c333e0330841c0d6b
筆者の環境で3.5時間程度（MacBook Pro (13-inch, 2016)、Memory: 8 GB）
さすがtriton。。。
補足
・DADA2はIllumina のアンプリコン配列を検出し修正するパイプラインである。Phix配列やキメラ配列をフィルタリングしてくれる。
https://qiita.com/kohei-108/items/547b2fbdf28fb04c28ea
・Qiime1とQiime2の違い
・サマリーを作成
code:bash
$ time qiime metadata tabulate \
--m-input-file stats-dada2.qza \
--o-visualization stats-dada2.qzv
Saved Visualization to: stats-dada2.qzv
real 0m4.489s
user 0m2.958s
sys 0m0.767s
↓
stats-dada2.qzvファイルをQIIME2 viewで開いて各種統計値を確認
https://gyazo.com/22bba4ff5da63b782266d30c2869750f
input（入力リード数）が各工程で段階的に減少し、nonchim（キメラ配列除去後）の値が最終的にFeature構築に使用されたリード数。
今回は平均して約半数の入力リードがnonchimリードとして残った。
＊キメラ配列
https://gyazo.com/7865e6ec6936d01a0a0d9990504fa681
https://gyazo.com/93b8538e3fab2f99aacc1d5cf6d7944d
https://gyazo.com/e74e862dca816b52edfaea7341d55a7f
https://www.drive5.com/usearch/manual/chimeras.html
メタゲノムは、16S用のプライマーでPCRにかけているから、キメラ配列が起きている。
Feature tableとFeatureDataの集計
code:bash
#サマリーの作成
$ time qiime feature-table summarize \
--i-table table.qza \
--o-visualization table.qzv \
--m-sample-metadata-file metadata.txt
＜エラー＞
code:bash
There was an issue with loading the file metadata.txt as metadata:
Metadata file path doesn't exist, or the path points to something other than a file. Please check that the path exists, has read permissions, and points to a regular file (not a directory): metadata.txt
There may be more errors present in the metadata file. To get a full report, sample/feature metadata files can be validated with Keemei: https://keemei.qiime2.org
Find details on QIIME 2 metadata requirements here: https://docs.qiime2.org/2019.7/tutorials/metadata/
https://gyazo.com/e072f0f4f56ab908480a6885dae8162e
解決
https://gyazo.com/a6defe2bb6a6e8e1c9371bc30183ceba
ファイル名を変更していなかった。。。
エラー②
code:bash
There was an issue with loading the file metadata.txt as metadata:
Metadata row contains more cells than are declared by the header. The row has 6 cells, while the header declares 4 cells.
There may be more errors present in the metadata file. To get a full report, sample/feature metadata files can be validated with Keemei: https://keemei.qiime2.org
Find details on QIIME 2 metadata requirements here: https://docs.qiime2.org/2019.7/tutorials/metadata/
変更前
https://gyazo.com/11154137871cdb18c08f37d5901d59bf
↓
変更後
https://gyazo.com/e0791603e17e56c7f6c882562bd44b3e
code:bash
1#SampleID host_diet host_genotype group
2DRR138864 AIN ob/ob ob-AIN
3DRR138865 AIN ob/ob ob-AIN
4DRR138866 AIN ob/ob ob-AIN
5DRR138867 AIN ob/ob ob-AIN
6DRR138868 AIN ob/ob ob-AIN
7DRR138869 AIN ob/ob ob-AIN
8DRR138870 AIN ob/ob ob-AIN
9DRR138871 K71 ob/ob ob-K71
10DRR138872 K71 ob/ob ob-K71
11DRR138873 K71 ob/ob ob-K71
12DRR138874 K71 ob/ob ob-K71
13DRR138875 K71 ob/ob ob-K71
14DRR138876 AIN ob/ob WT-AIN
15DRR138877 AIN ob/ob WT-AIN
16DRR138878 AIN ob/ob WT-AIN
17DRR138879 K71 wild-type ob-K71
18DRR138880 K71 wild-type ob-K71
19DRR138881 K71 wild-type ob-K71
20DRR138882 AIN wild-type WT-AIN
21DRR138883 AIN wild-type WT-AIN
22DRR138884 AIN wild-type WT-AIN
23DRR138885 AIN wild-type WT-AIN
24DRR138886 AIN wild-type WT-AIN
↓
実行
code:bash
#サマリーの作成
$ time qiime feature-table summarize \
--i-table table.qza \
--o-visualization table.qzv \
--m-sample-metadata-file metadata.txt
Saved Visualization to: table.qzv
real 0m5.413s
user 0m3.770s
sys 0m0.772s
code:bash
$ time qiime feature-table tabulate-seqs \
--i-data rep-seqs.qza \
--o-visualization rep-seqs.qzv
Saved Visualization to: rep-seqs.qzv
real 0m3.948s
user 0m3.110s
sys 0m0.657s
table.qzvファイルをQIIME2 viewで開いて、Feature tableの各種統計値を確認
https://gyazo.com/f4c2093817d9ee5409726fc0e2ef0a44
https://gyazo.com/7e8b554616487c472a42d42a217bbb65
https://gyazo.com/9f244c301e8b3b06d6a8d0144ecb8fae
https://gyazo.com/4fe93d30880a5c8b453d090d519e5b08
https://gyazo.com/d4a4eccc95990ae0b3d17c524e93c0ac
Number of featuresがいわゆるOTU数で、今回のデータセットでは663である。Interactive Sample Detailタブではレアファクション処理の際に設定するリード数（サンプリング深度）の検証が行える。今回のデータセットではリード数が最小のサンプルは43,591リードであることがわかる。後ほど設定するサンプリング深度をこの値に設定する。Feature Detailタブでは、FeatureのIDと属するリード数、出現するサンプル数が表示される。
・featureは、OTU、シーケンスバリアント、遺伝子、代謝産物などの観測の任意の単位。
https://forum.qiime2.org/t/what-is-a-feature-exactly/2201
・OTU･･･operational taxonomic unit。細菌の必須遺伝子（一般に，16SリボソームRNA遺伝子）の塩基配列をコンピュータ上でその類似度を指標に分類したときに得られる単位。
・depth
https://gyazo.com/779a104d282acf2e73db1413832dc029
rep-seqs.qzvファイルをQIIME2 viewで開いて、Featureの代表配列を確認。
塩基配列をクリックすると、直接NCBIのBLASTで相同性検索が行われる。
https://gyazo.com/a689652ea4e2b0e9c376b89592401941
分子系統樹の計算
code:bash
$ time 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
Saved FeatureDataAlignedSequence to: aligned-rep-seqs.qza
Saved FeatureDataAlignedSequence to: masked-aligned-rep-seqs.qza
Saved PhylogenyUnrooted to: unrooted-tree.qza
Saved PhylogenyRooted to: rooted-tree.qza
real 0m20.193s
user 0m19.106s
sys 0m0.866s
α多様性とβ多様性の解析
α多様性（サンプル内での菌種数）とβ多様性（サンプル間距離）の解析を行う。計算の過程でサンプル間のリード数を揃えるレアファクション処理が行われる。サンプリング深度は「Feature tableとFeatureDataの集計」で求めたリード数が最小のサンプルの値を指定する。サンプリング深度未満のリード数のサンプルは多様性解析結果の出力からは除かれる。
code:bash
#説明
$ qiime diversity core-metrics-phylogenetic \
--i-phylogeny rooted-tree.qza \
--i-table table.qza \
--p-sampling-depth 43591 \（レアファクション処理のサンプリング深度）
--m-metadata-file metadata.txt \
--output-dir core-metrics-results
#実行
$ time qiime diversity core-metrics-phylogenetic \
--i-phylogeny rooted-tree.qza \
--i-table table.qza \
--p-sampling-depth 43591 \
--m-metadata-file metadata.txt \
--output-dir core-metrics-results
Saved FeatureTableFrequency to: core-metrics-results/rarefied_table.qza
Saved SampleDataAlphaDiversity % Properties('phylogenetic') to: core-metrics-results/faith_pd_vector.qza
Saved SampleDataAlphaDiversity to: core-metrics-results/observed_otus_vector.qza
Saved SampleDataAlphaDiversity to: core-metrics-results/shannon_vector.qza
Saved SampleDataAlphaDiversity to: core-metrics-results/evenness_vector.qza
Saved DistanceMatrix % Properties('phylogenetic') to: core-metrics-results/unweighted_unifrac_distance_matrix.qza
Saved DistanceMatrix % Properties('phylogenetic') to: core-metrics-results/weighted_unifrac_distance_matrix.qza
Saved DistanceMatrix to: core-metrics-results/jaccard_distance_matrix.qza
Saved DistanceMatrix to: core-metrics-results/bray_curtis_distance_matrix.qza
Saved PCoAResults to: core-metrics-results/unweighted_unifrac_pcoa_results.qza
Saved PCoAResults to: core-metrics-results/weighted_unifrac_pcoa_results.qza
Saved PCoAResults to: core-metrics-results/jaccard_pcoa_results.qza
Saved PCoAResults to: core-metrics-results/bray_curtis_pcoa_results.qza
Saved Visualization to: core-metrics-results/unweighted_unifrac_emperor.qzv
Saved Visualization to: core-metrics-results/weighted_unifrac_emperor.qzv
Saved Visualization to: core-metrics-results/jaccard_emperor.qzv
Saved Visualization to: core-metrics-results/bray_curtis_emperor.qzv
real 0m5.714s
user 0m4.441s
sys 0m1.013s
各種α多様性指数とメタデータとの関連解析
code:bash
$ time qiime diversity alpha-group-significance \
--i-alpha-diversity core-metrics-results/observed_otus_vector.qza \
--m-metadata-file metadata.txt \
--o-visualization core-metrics-results/observed_otus-group-significance.qzv
Saved Visualization to: core-metrics-results/observed_otus-group-significance.qzv
real 0m3.569s
user 0m2.898s
sys 0m0.618s
$ time qiime diversity alpha-group-significance \
--i-alpha-diversity core-metrics-results/shannon_vector.qza \
--m-metadata-file metadata.txt \
--o-visualization core-metrics-results/shannon-group-significance.qzv
Saved Visualization to: core-metrics-results/shannon-group-significance.qzv
real 0m4.339s
user 0m2.917s
sys 0m0.700s
$ time qiime diversity alpha-group-significance \
--i-alpha-diversity core-metrics-results/faith_pd_vector.qza \
--m-metadata-file metadata.txt \
--o-visualization core-metrics-results/faith-pd-group-significance.qzv
Saved Visualization to: core-metrics-results/faith-pd-group-significance.qzv
real 0m3.485s
user 0m2.882s
sys 0m0.553s
各種β多様性指数とメタデータとの関連解析
code:bash
#説明
$ qiime diversity beta-group-significance \
--i-distance-matrix core-metrics-results/unweighted_unifrac_distance_matrix.qza \
--m-metadata-file metadata.txt \
--m-metadata-column group \（群分けに使うメタデータの変数名）
--o-visualization core-metrics-results/unweighted-unifrac-group-significance.qzv \
--p-pairwise
#実行
$ time qiime diversity beta-group-significance \
--i-distance-matrix core-metrics-results/unweighted_unifrac_distance_matrix.qza \
--m-metadata-file metadata.txt \
--m-metadata-column group \
--o-visualization core-metrics-results/unweighted-unifrac-group-significance.qzv \
--p-pairwise
Saved Visualization to: core-metrics-results/unweighted-unifrac-group-significance.qzv
real 0m5.562s
user 0m3.855s
sys 0m0.853s
#説明
$ qiime diversity beta-group-significance \
--i-distance-matrix core-metrics-results/weighted_unifrac_distance_matrix.qza \
--m-metadata-file metadata.txt \
--m-metadata-column group \（群分けに使うメタデータの変数名）
--o-visualization core-metrics-results/weighted-unifrac-group-significance.qzv \
--p-pairwise
#実行
$ time qiime diversity beta-group-significance \
--i-distance-matrix core-metrics-results/weighted_unifrac_distance_matrix.qza \
--m-metadata-file metadata.txt \
--m-metadata-column group \
--o-visualization core-metrics-results/weighted-unifrac-group-significance.qzv \
--p-pairwise
Saved Visualization to: core-metrics-results/weighted-unifrac-group-significance.qzv
real 0m4.469s
user 0m3.777s
sys 0m0.636s
各種α/β多様性指数の解析結果がcore-metrics-resultsフォルダ内に出力されるのでブラウザのQIIME2 viewで閲覧
α多様性の解析では各種の多様性指数をメタデータの変数で群間比較ができる。β多様性の解析ではサンプル間の類似性を3次元プロットで視覚的に表示したり、サンプル間距離の検定結果が得られる。
https://gyazo.com/36978d5798f78e700679bc43e481baf5
たとえば、
observed_otus-group-significance.qzv
https://gyazo.com/77a7bbff0c3c4e5a292300cfdd6d3861
weighted_unifrac_emperor.qzv
https://gyazo.com/cf90ab2a453d3e545474f3a20ee7024c
α-レアファクションカーブの作図
横軸にサンプリングしたリード深度、縦軸にα多様性指数をプロットしたカーブ、α-レアファクションカーブ。
code:bash
#説明
$ qiime diversity alpha-rarefaction \
--i-table table.qza \
--i-phylogeny rooted-tree.qza \
--p-max-depth 114909 \（最大サンプリング深度、最大リード数のサンプルの値）
--m-metadata-file metadata.txt \
--o-visualization alpha-rarefaction.qzv
#実行
$ time qiime diversity alpha-rarefaction \
--i-table table.qza \
--i-phylogeny rooted-tree.qza \
--p-max-depth 114909 \
--m-metadata-file metadata.txt \
--o-visualization alpha-rarefaction.qzv
＜エラー＞
code:bash
Plugin error from diversity:
Provided max_depth of 114909 is greater than the maximum sample total frequency of the feature_table (114279).
Debug info has been saved to /tmp/qiime2-q2cli-err-ud1m56i9.log
↓
数字を変更する。
code:bash
$ time qiime diversity alpha-rarefaction \
--i-table table.qza \
--i-phylogeny rooted-tree.qza \
--p-max-depth 114279 \
--m-metadata-file metadata.txt \
--o-visualization alpha-rarefaction.qzv
Saved Visualization to: alpha-rarefaction.qzv
real 0m8.753s
user 0m7.080s
sys 0m0.813s
alpha-rarefaction.qzvファイルをQIIME2 viewで開いて確認
https://gyazo.com/421e79c69ba9b0c1d5816672252b6e41
・Sequencing Depthとは
http://www.iu.a.u-tokyo.ac.jp/~kadota/20140625_kadota.pdf
Coverage（カバレッジ）
ゲノム解読したいときなどに、解読するために必要とされる指標となる数値。ゲノムサイズ(X)に
対する、sequencerで読んだ塩基配列長の和のこと。一般に、この数値が高いほどよい。
Sequence depthという表現と実質的に同じような指標。
https://gyazo.com/52adcca4620795277cfe90655a550791
https://bi.biopapyrus.jp/rnaseq/exp/replicate-sequencing-depth.html
系統解析
Featureの代表配列を16S rRNA遺伝子リファレンスデータベースに対して相同性検索を行い菌種同定を行う。準備段階でQIIME 2のWebサイトからダウンロードしたGreengenesのNaive Bayes classifier用トレーニングファイルを用いる。
code:bash
#説明
$ qiime feature-classifier classify-sklearn \
--p-n-jobs -1 \（使用可能な全CPUコアを使う）
--i-classifier gg-13-8-99-nb-classifier.qza \
--i-reads rep-seqs.qza \
--o-classification taxonomy.qza
#実行
$ time qiime feature-classifier classify-sklearn \
--p-n-jobs -1 \
--i-classifier gg-13-8-99-nb-classifier.qza \
--i-reads rep-seqs.qza \
--o-classification taxonomy.qza
Saved FeatureDataTaxonomy to: taxonomy.qza
real 0m25.452s
user 8m18.832s
sys 1m30.597s
#閲覧できるように変換
$ time qiime metadata tabulate \
--m-input-file taxonomy.qza \
--o-visualization taxonomy.qzv
Saved Visualization to: taxonomy.qzv
real 0m3.822s
user 0m2.946s
sys 0m0.674s
$ time qiime taxa barplot \
--i-table table.qza \
--i-taxonomy taxonomy.qza \
--m-metadata-file metadata.txt \
--o-visualization taxa-bar-plots.qzv
Saved Visualization to: taxa-bar-plots.qzv
real 0m4.675s
user 0m2.994s
sys 0m0.826s
taxa-bar-plots.qzvファイルをQIIME2 viewで確認。
菌種同定結果が棒グラフで表示される。
Taxonomic Levelで分類階層を選択できる。数値データをCSVファイルに出力することもできる。
https://gyazo.com/7149d035f21da5a75ad50d7300e0ebd9
Taxonomic Levelを変更して、どれが一番上手くグラフとして差がでているかを手動で調べる。
ヒートマップの作図
Feature tableから任意の分類階級（今回は科レベル）でヒートマップを描いてみる。
code:bash
#説明
$ qiime taxa collapse \
--i-table table.qza \
--i-taxonomy taxonomy.qza \
--p-level 5 \（分類階級 1:Kingdom; 2:Phylum; 3:Class; 4:Order; 5:Family; 6:Genus）
--o-collapsed-table table-l5.qza
#実行
$ time qiime taxa collapse \
--i-table table.qza \
--i-taxonomy taxonomy.qza \
--p-level 5 \
--o-collapsed-table table-l5.qza
Saved FeatureTableFrequency to: table-l5.qza
real 0m3.243s
user 0m2.743s
sys 0m0.453s
#説明
$ qiime feature-table heatmap \
--i-table table-l5.qza \
--m-metadata-file metadata.txt \
--m-metadata-column group \（ヒートマップに表示するメタデータの変数名）
--o-visualization heatmap_l5_group.qzv
#実行
$ time qiime feature-table heatmap \
--i-table table-l5.qza \
--m-metadata-file metadata.txt \
--m-metadata-column group \
--o-visualization heatmap_l5_group.qzv
Saved Visualization to: heatmap_l5_group.qzv
real 0m5.903s
user 0m4.922s
sys 0m0.921s
heatmap_l5_group.qzvファイルをQIIME2 viewで開くとヒートマップが表示
https://gyazo.com/dcd6dc338bb9e9a0216e44bb854d592c
変動菌種の検出
Feature tableから2群間で有意に存在量が変動した菌種を指定した分類階級で検出してみる。QIIME 2にはANCOMとgneissの二つのアルゴリズムが実装されているが本稿ではANCOMを用いる方法を紹介する。ANCOMは群間で変動するFeatureが少数（25%未満）であることを前提としているので、より多くのFeatureが変動が予想される場合は他の手法を検討すべきである。今回はob/obマウスにおいて乳酸菌摂取により有意に変動した菌種を科レベルで同定してみる。最初に以下のコマンドを入力し、Feature tableからob/obマウスだけを抽出する。
code:bash
$ qiime feature-table filter-samples \
--i-table table.qza \
--m-metadata-file metadata.txt \
--p-where "host_genotype='ob'" \
--o-filtered-table table_ob.qza
Saved FeatureTableFrequency to: table_ob.qza
code:bash
#ANCOMによる統計解析
#説明
$ qiime taxa collapse \
--i-table table_ob.qza \
--i-taxonomy taxonomy.qza \
--p-level 5 \（分類階級 1:Kingdom; 2:Phylum; 3:Class; 4:Order; 5:Family; 6:Genus）
--o-collapsed-table table_ob_l5.qza
#実行
$ time qiime taxa collapse \
--i-table table_ob.qza \
--i-taxonomy taxonomy.qza \
--p-level 5 \
--o-collapsed-table table_ob_l5.qza
＜エラー＞
Please provide a DataFrame with a string-based Index
↓
ob/obではなく、obで！
再度
code:bash
$ qiime feature-table filter-samples \
--i-table table.qza \
--m-metadata-file metadata.txt \
--p-where "host_genotype='ob'" \
--o-filtered-table table_ob.qza
#実行
$ time qiime taxa collapse \
--i-table table_ob.qza \
--i-taxonomy taxonomy.qza \
--p-level 5 \
--o-collapsed-table table_ob_l5.qza
Saved FeatureTableFrequency to: table_ob_l5.qza
real 0m3.216s
user 0m2.689s
sys 0m0.476s
code:bash
$ time qiime composition add-pseudocount \
--i-table table_ob_l5.qza \
--o-composition-table comp_table_ob_l5.qza
Saved FeatureTableComposition to: comp_table_ob_l5.qza
real 0m3.458s
user 0m2.714s
sys 0m0.547s
#説明
$ qiime composition ancom \
--i-table comp_table_ob_l5.qza \
--m-metadata-file metadata.txt \
--m-metadata-column host_diet \（比較したいメタデータの変数名）
--o-visualization ancom_table_ob_l5_diet.qzv
#実行
$ time qiime composition ancom \
--i-table comp_table_ob_l5.qza \
--m-metadata-file metadata.txt \
--m-metadata-column host_diet \
--o-visualization ancom_table_ob_l5_diet.qzv
Saved Visualization to: ancom_table_ob_l5_diet.qzv
real 0m3.604s
user 0m2.950s
sys 0m0.590s
ancom_table_ob_l5_diet.qzvファイルをQIIME2 viewで開くとVolcano Plot、変動菌種とその存在量が表示
https://gyazo.com/b653b3d869a5318f10cd4516d1022d89
まとめ
QIIME 2はチュートリアル等の公式ドキュメントが充実。