本示例的的數據來自文章《Moving pictures of the human microbiome》，Genome Biology 2011，取樣來自兩個人身體四個部位五個時間點
進入環境
source activate qiime 2-2017.8
退出環境
source deactivate
准備數據
# 創建並進入工作目錄
mkdir -p qiime2-moving-pictures-tutorial
cd qiime2-moving-pictures-tutorial
# 下載實驗設計（-O 重命名下載的文件）
wget -O sample-metadata.tsv https://data.qiime2.org/2017.6/tutorials/moving-pictures/sample_metadata.tsv
## 上面一步下載失敗，可嘗試刪除空文件並使用我建立的備份鏈接下載；否則跳過下面兩行命令
rm sample_metadata.tsv
wget http://bailab.genetics.ac.cn/markdown/sample-metadata.tsv
# 下載實驗測序數據
mkdir -p emp-single-end-sequences
wget -O emp-single-end-sequences/barcodes.fastq.gz https://data.qiime2.org/2017.6/tutorials/moving-pictures/emp-single-end-sequences/barcodes.fastq.gz
wget -O emp-single-end-sequences/sequences.fastq.gz https://data.qiime2.org/2017.6/tutorials/moving-pictures/emp-single-end-sequences/sequences.fastq.gz
# 生成qiime需要的artifact文件(qiime文件格式，將原始數據格式標准化)
qiime tools import \
--type EMPSingleEndSequences \
--input-path emp-single-end-sequences \
--output-path emp-single-end-sequences.qza
拆分樣品
# 按barcode拆分樣品Demultiplexing sequences
qiime demux emp-single \
--i-seqs emp-single-end-sequences.qza \
--m-barcodes-file sample-metadata.tsv \
--m-barcodes-category BarcodeSequence \
--o-per-sample-sequences demux.qza
# 結果統計
qiime demux summarize \
--i-data demux.qza \
--o-visualization demux.qzv
# 查看結果 (依賴XShell+XManager或其它ssh終端和圖形界面軟件)
qiime tools view demux.qzv
序列質控和生成OTU表
此步主要有DADA2和Deblur兩種方法可選，
推薦使用DADA2，去年發表在Nature Method上，比較同類方法優於其它OTU聚類結果；相較QIIME的UPARSE聚類方法，目前DADA2方法僅去噪去嵌合，不再按相似度聚類。比上一代分析結果更准確。
DADA2
主要作用是去除低質量序列、嵌合體；再生成
OTU表，現在叫Feature表，因為不再使用聚類方法，相當於QIIME時代100%相似度的OTU表。
讀者思考時間：基於上面對拆分樣品的統計結果，如何設置下面生成OTU表的參數。
–p-trim-left 截取左端低質量序列，我們看圖中箱線圖，左端質量都很高，無低質量區，設置為0；
–p-trunc-len 序列截取長度，也是為了去除右端低質量序列，我們看到大於120以后，質量下降極大，甚至中位數都下降至20以下，需要全部去除。
# 單端序列去噪, 去除左端0bp(--p-trim-left用於切除邊緣低質量區)，序列切成120bp長；生成代表序列和OTU表；並重命名用於下游分析
qiime dada2 denoise-single \
--i-demultiplexed-seqs demux.qza \
--p-trim-left 0 \
--p-trunc-len 120 \
--o-representative-sequences rep-seqs-dada2.qza \
--o-table table-dada2.qza
mv rep-seqs-dada2.qza rep-seqs.qza
mv table-dada2.qza table.qza
Deblur
與DADA2二選一，用戶可自行比較結果的差異，根據喜好選擇
# 按測序質量過濾序列
qiime quality-filter q-score \
--i-demux demux.qza \
--o-filtered-sequences demux-filtered.qza \
--o-filter-stats demux-filter-stats.qza
# 去冗余生成OTU表和代表序列；結果文件名有deblur，沒有用於下游分析，請讀者想測試的自己嘗試
qiime deblur denoise-16S \
--i-demultiplexed-seqs demux-filtered.qza \
--p-trim-length 120 \
--o-representative-sequences rep-seqs-deblur.qza \
--o-table table-deblur.qza \
--o-stats deblur-stats.qza
Feature表統計
qiime feature-table summarize \
--i-table table.qza \
--o-visualization table.qzv \
--m-sample-metadata-file sample-metadata.tsv
qiime tools view table.qzv
代表序列統計
qiime feature-table tabulate-seqs \
--i-data rep-seqs.qza \
--o-visualization rep-seqs.qzv
qiime tools view rep-seqs.qzv
建樹：用於多樣性分析
# 多序列比對
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
# 建樹
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
Alpha多樣性
讀者思考時間：下面多樣性分析，需要基於標准化的OTU表，標准化采用重抽樣至序列一致，如何設計樣品重抽樣深度參數。–p-sampling-depth
如是數據量都很大，選最小的即可。如果有個別數據量非常小，去除最小值再選最小值。比如此分析最小值為917，我們選擇1080深度重抽樣，即保留了大部分樣品用於分析，又去除了數據量過低的異常值。
注：本示例為454時代的測序，數據量很小。現在一般采用HiSeq PE250測序，數據量都非常大，通常可以采用3萬或5萬的標准篩選，仍可保留90%以上樣本。過低或過高一般結果也會異常，不建議放在一起分析。
# 計算多樣性(包括所有常用的Alpha和Beta多樣性方法)，輸入有根樹、Feature表、樣本重采樣深度(一般為最小樣本數據量，或覆蓋絕大多數樣品的數據量)
qiime diversity core-metrics \
--i-phylogeny rooted-tree.qza \
--i-table table.qza \
--p-sampling-depth 1080 \
--output-dir core-metrics-results
# 輸出結果包括多種多樣性結果，文件列表和解釋如下：
# beta多樣性bray_curtis距離矩陣 bray_curtis_distance_matrix.qza
# alpha多樣性evenness(均勻度，考慮物種和豐度)指數 evenness_vector.qza
# alpha多樣性faith_pd(考慮物種間進化關系)指數 faith_pd_vector.qza
# beta多樣性jaccard距離矩陣 jaccard_distance_matrix.qza
# alpha多樣性observed_otus(OTU數量)指數 observed_otus_vector.qza
# alpha多樣性香農熵(考慮物種和豐度)指數 shannon_vector.qza
# beta多樣性unweighted_unifrac距離矩陣，不考慮豐度 unweighted_unifrac_distance_matrix.qza
# beta多樣性unweighted_unifrac距離矩陣，考慮豐度 weighted_unifrac_distance_matrix.qza
# faith_pd算法統計Alpha多樣性組間差異是否顯著，輸入多樣性值、實驗設計，輸出統計結果
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
# 統計evenness組間差異是否顯著
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
# 網頁展示結果，只要是qzv的文件，均可用qiime tools view查看或在線
https://view.qiime2.org/查看，以后不再贅述
qiime tools view core-metrics-results/evenness-group-significance.qzv
讀者思考時間：實驗設計中的那一種分組方法，與微生物群體的豐富度差異相關，這些差異顯著嗎？
解答：圖中可按Catalogy選擇分類方法，查看不同分組下箱線圖間的分布與差別。圖形下面的表格，詳細詳述組間比較的顯著性和假陽性率統計。
結果我們會看到本實驗設計的分組方式有Bodysite, Subject, ReportAntibioticUse，只有身體位置各組間差異明顯，且下面統計結果也存在很多組間的顯著性差異
Beta多樣性
# 按BodySite分組，統計unweighted_unifrace距離的組間是否有顯著差異
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 BodySite \
--o-visualization core-metrics-results/unweighted-unifrac-body-site-significance.qzv \
--p-pairwise
# 按Subject分組，統計unweighted_unifrace距離的組間是否有顯著差異
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 Subject \
--o-visualization core-metrics-results/unweighted-unifrac-subject-group-significance.qzv \
--p-pairwise
# 可視化三維展示unweighted-unifrac的主坐標軸分析
qiime emperor plot \
--i-pcoa core-metrics-results/unweighted_unifrac_pcoa_results.qza \
--m-metadata-file sample-metadata.tsv \
--p-custom-axis DaysSinceExperimentStart \
--o-visualization core-metrics-results/unweighted-unifrac-emperor.qzv
# 可視化三維展示bray-curtis的主坐標軸分析
qiime emperor plot \
--i-pcoa core-metrics-results/bray_curtis_pcoa_results.qza \
--m-metadata-file sample-metadata.tsv \
--p-custom-axis DaysSinceExperimentStart \
--o-visualization core-metrics-results/bray-curtis-emperor.qzv
# 網頁展示結果，或下載在線查看
qiime tools view core-metrics-results/bray-curtis-emperor.qzv
讀者思考時間：按subject分組有顯著區別嗎？按body-site分組有顯著區別嗎？那些body-site組間存在區別？
按其它距離計算的結果，讀者可以仔細看看不同距離矩陣計算結果的區別。個人感覺，一般比較好解釋科學問題的方法就是適合的方法
物種分類
# 下載物種注釋
wget -O gg- 13-8-99-515-806-nb-classifier.qza https://data.qiime2.org/2017.6/common/gg-13-8-99-515-806-nb-classifier.qza
# 物種分類
qiime feature-classifier classify-sklearn \
--i-classifier gg- 13-8-99-515-806-nb-classifier.qza \
--i-reads rep-seqs.qza \
--o-classification taxonomy.qza
# 物種結果轉換表格，可用於查看
qiime taxa tabulate \
--i-data taxonomy.qza \
--o-visualization taxonomy.qzv
# 物種分類柱狀圖
qiime taxa barplot \
--i-table table.qza \
--i-taxonomy taxonomy.qza \
--m-metadata-file sample-metadata.tsv \
--o-visualization taxa-bar-plots.qzv
# 網頁展示結果，或下載在線查看
qiime tools view taxa-bar-plots.qzv
讀者思考時間1：代表序列文件rep-seqs.qzv可視化結果中，可以下載fasta文件采用NCBI進行blast注釋物種信息，與我們目前的結果比較，看看有什么不同，各分類級別的注釋定義的相似程度是什么？
讀者思考時間2：查看門水平(level2)分類結果柱狀圖，看每一類body-site中主要豐度的門類是什么？
差異豐度分析
差異豐度分析采用ANCOM (analysis of composition of microbiomes)，是2015年發布在Microb Ecol Health Dis上的方法，文章稱在微生物組方面更專業，但不接受零值(零在二代測序結果表中很常見)。我個人一直用edgeR，感覺靠譜，因為高通量測序本質上是相同的
差異Features/OTUs分析
# OTU表添加假count，因為ANCOM不允許有零
qiime composition add-pseudocount \
--i-table table.qza \
--o-composition-table comp-table.qza
# 采用ancon，按BodySite分組進行差異統計
qiime composition ancom \
--i-table comp-table.qza \
--m-metadata-file sample-metadata.tsv \
--m-metadata-category BodySite \
--o-visualization ancom-BodySite.qzv
# 查看結果
qiime tools view ancom-BodySite.qzv
讀者思考時間：不同身體部分有那些Features存在豐度差異？那一組是最高或最低豐度？這此差異的Features屬那些分類單元？
++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
差異分類學級別分析：以按門水平合並再統計差異
# 按門水平進行合並，統計各門的總reads
qiime taxa collapse \
--i-table table.qza \
--i-taxonomy taxonomy.qza \
--p-level 2 \
--o-collapsed-table table-l2.qza
# 同理去除零
qiime composition add-pseudocount \
--i-table table-l2.qza \
--o-composition-table comp-table-l2.qza
# 在門水平按取樣部分分析
qiime composition ancom \
--i-table comp-table-l2.qza \
--m-metadata-file sample-metadata.tsv \
--m-metadata-category BodySite \
--o-visualization l2-ancom-BodySite.qzv
讀者思考時間：不同身體部分有那些Features存在豐度差異？那一組是最高或最低豐度？這此差異的Features屬那些分類單元？
結果描述：結果的可視化(Visual)頁面，一共分為三部分。
第一個表為ANCOM statistical results，只列出組間存在顯著差異的門，其統計值W的計算及解釋尚不清楚，查原始文章也沒有找到。有待更新版中解釋。
第二個表為各組的豐度分位數，就是箱線圖的原始數據，為什么作者沒有直接出圖，我將與作者溝通討論；目前可以比較各組的分布，來具體分析組間的差異，但不夠直觀；
統計各門類的火山圖，坐標軸還沒有詳細解釋，但其意思是越靠上越顯著差異。此圖采用Python的bokeh庫生成的交互式圖形，可以點擊圖中的點來查看具體的詳細，如具體的分類學信息。相當於表1的可視化。
結果的網頁還有其它頁面，如peek頁面可以查看此文件的基本信息，Provenance頁面顯示當前結果的生成過程圖，點擊過程中的點可以查看具體的程序和參數；鏈接按扭可以生成共享鏈接；下載按扭可以下載原始文件。