本示例的的数据来自文章《Moving pictures of the human microbiome》，Genome Biology 2011，取样来自两个人身体四个部位五个时间点
进入环境
source activate qiime2-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页面显示当前结果的生成过程图，点击过程中的点可以查看具体的程序和参数；链接按扭可以生成共享链接；下载按扭可以下载原始文件。