查看原文
其他

扩增子分析QIIME2. 2分析实战Moving Pictures

2017-07-26 刘永鑫 宏基因组

欢迎点击「宏基因组」关注我们!专业干货每日推送!

声明:本文为QIIME2官方帮助文档的中文版,由中科院遗传发育所刘永鑫博士翻译并亲测有效,文档翻译己获QIIME2团队官方授权。由于QIIME2更新频繁,如使用中遇到问题请访问QIIME2官方论坛阅读最新版中文帮助。
https://forum.qiime2.org/t/qiime2-1-chinese-manual/838
如中文翻译没有急时更新,新阅读英文原版 https://docs.qiime2.org

本人只习惯使用命令行模式分析数据,图形界面和Ipython模式下使用暂不介绍。本系列的教程主要以命令行方式为大家演示。在其它方面有使用的经验的朋友,欢迎共享您的使用笔记于“宏基因组”公众号,方便大家学习和使用。

扩增子分析QIIME2. 2分析实战:人体不同部分微生物组Moving Pictures

QIIME2版本 2017.6

本文的操作的前提是完成QIIME2的安装,想安装QIIME2请阅读《扩增子分析QIIME2. 1简介和安装》阅读;

本示例的的数据来自文章《Moving pictures of the human microbiome》,Genome Biology 2011,取样来自两个人,身体四个部位,五个时间点。

启动QIIME2运行环境

对于上文提到了两种常用安装方法,我们每次在分析数据前,需要打开工作环境,根据情况选择对应的打开方式。

# 创建工作目录并进入 mkdir -p qiime2 cd qiime2 # Miniconda安装方式的请运行如下命令加载工作环境 source activate qiime2-2017.6 # 如果是docker安装方式的请运行如下命令,否则跳过此行 docker run --rm -v $(pwd):/data --name=qiime -it  qiime2/core:2017.6

准备数据

# 创建并进入工作目录 mkdir -p qiime2-moving-pictures-tutorial cd qiime2-moving-pictures-tutorial # 下载实验设计,此文件位于google服务器,服务器无法翻墙的可以用windows翻墙下载并上传到工作目录;或尝试下载我建立的备份 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

使用qiime tools view demux.qzv或访问 https://view.qiime2.org/?src=https%3A%2F%2Fdocs.qiime2.org%2F2017.6%2Fdata%2Ftutorials%2Fmoving-pictures%2Fdemux.qzv 在线查看
下图展示拆分样品的图表结果,其中的图表示各样品测序数据柱状分布图,展示不同测序深度下样品数量分布,可下载PDF文件;表可下载CSV格式文件,即每个样品的数据量表,可用Excel打开。



序列质控和生成OTU表

此步主要有DADA2和Deblur两种方法可选,推荐使用DADA2,去年发表在Nature Method上,比较同类方法结果(基于分析mock community模拟重组数据)优于其它OTU聚类结果,结果更真实;相较QIIME的UPARSE聚类方法,目前DADA2方法仅去噪去嵌合,不再按相似度聚类,相当于100%相似度的OTU。比上一代分析结果更准确。

  1. DADA2
    主要作用是去除低质量序列、嵌合体;再生成OTU表,现在叫Feature表,因为不再使用聚类方法,相当于QIIME时代100%相似度的OTU表。

读者思考时间:基于上图对拆分样品的统计结果,如何设置下面生成OTU表的参数。

  1. —p-trim-left 截取左端低质量序列,我们看上图中箱线图,左端质量都很高,无低质量区,设置为0;

  2. —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

2. 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 # 展打开网页版统计结果,看下图

图中展示了Feature表的统计结果


代表序列统计

qiime feature-table tabulate-seqs \  --i-data rep-seqs.qza \  --o-visualization rep-seqs.qzv qiime tools view rep-seqs.qzv # 展打开网页版统计结果,看下图

下图展示代表序列网页版详细,点"Click here"可下载fasta文件

建树:用于多样性分析

# 多序列比对 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 evenness-group-significance.qzv

以Evenness为例,下图展示多样性的箱线图,可以下载svg矢量图



也可点击如下网址,查看可交互结果。https://view.qiime2.org/visualization/?type=html&src=https%3A%2F%2Fdocs.qiime2.org%2F2017.6%2Fdata%2Ftutorials%2Fmoving-pictures%2Fcore-metrics-results%2Fevenness-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

以Bray_curtis可视化结果为例,下图展示多样性的三维散点图




读者思考时间:按subject分组有显著区别吗?按body-site分组有显著区别吗?那些body-site组间存在区别?

大家仔细查看按subject进行主坐标轴分析的结果,如下链接,我是没看到区别;
https://view.qiime2.org/?src=https%3A%2F%2Fdocs.qiime2.org%2F2017.6%2Fdata%2Ftutorials%2Fmoving-pictures%2Fcore-metrics-results%2Funweighted-unifrac-subject-group-significance.qzv
大家仔细查看按body-site进行主坐标轴分析的结果,如下链接,只有左右手间无明显差别,其它各组间均有显著差别。
https://view.qiime2.org/?src=https%3A%2F%2Fdocs.qiime2.org%2F2017.6%2Fdata%2Ftutorials%2Fmoving-pictures%2Fcore-metrics-results%2Funweighted-unifrac-body-site-significance.qzv

按其它距离计算的结果,读者可以仔细看看不同距离矩阵计算结果的区别。个人感觉,一般比较好解释科学问题的方法就是适合的方法。

物种分类

# 下载物种注释 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 # 展示taxonomy.qzv结果如下: #Feature ID    Taxonomy #d12759fe8dda1d65fe9077cc1ca9cf28    k__Bacteria; p__Bacteroidetes; c__Flavobacteriia; o__Flavobacteriales; f__[Weeksellaceae]; g__Chryseobacterium; s__ #5ada68b9a081358e1a7d5f1d351e656a    k__Bacteria; p__Fusobacteria; c__Fusobacteriia; o__Fusobacteriales; f__Leptotrichiaceae; g__Leptotrichia; s__ #d9095748835ade1b8914c5f57b6acbcf    k__Bacteria; p__Proteobacteria; c__Gammaproteobacteria; o__Aeromonadales; f__Aeromonadaceae; g__Oceanisphaera; s__ # 物种分类柱状图 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

以堆叠柱状图,展示物种各组、各分类级别下的相对丰度,图形可交互查看每个柱的详细信息,也可以下载svg矢量图。

读者思考时间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属那些分类单元?

使用qiime tools view ancom-BodySite.qzv查看结果,或下载结果用qiime2网页在线打开查看,回答上面的问题。

差异分类学级别分析:以按门水平合并再统计差异

# 按门水平进行合并,统计各门的总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属那些分类单元?
使用qiime tools view ancom-BodySite.qzv查看结果,或点如下链接查看详细,回答上面的问题。
https://view.qiime2.org/?src=https%3A%2F%2Fdocs.qiime2.org%2F2017.6%2Fdata%2Ftutorials%2Fmoving-pictures%2Fl2-ancom-BodySite.qzv

结果描述:结果的可视化(Visual)页面,一共分为三部分。

  1. 第一个表为ANCOM statistical results,只列出组间存在显著差异的门,其统计值W的计算及解释尚不清楚,查原始文章也没有找到。有待更新版中解释。

  2. 第二个表为各组的丰度分位数,就是箱线图的原始数据,为什么作者没有直接出图,我将与作者沟通讨论;目前可以比较各组的分布,来具体看组间的差异,但不够直观;

  3. 统计各门类的火山图,坐标轴还没有详细解释,但其意思是越靠上越显著差异。此图采用Python的bokeh库生成的交互式图形,可以点击图中的点来查看具体的详细,如具体的分类学信息。相当于表1的可视化。

结果的网页还有其它页面,如peek页面可以查看此文件的基本信息,Provenance页面显示当前结果的生成过程图,点击过程中的点可以查看具体的程序和参数;链接按扭可以生成共享链接;下载按扭可以下载原始文件。

Reference

  1. https://docs.qiime2.org/2017.6/tutorials/moving-pictures/

  2. Microb Ecol Health Dis. 2015 May 29;26:27663. doi: 10.3402/mehd.v26.27663. eCollection 2015. https://www.ncbi.nlm.nih.gov/pubmed/26028277

  3. Nat Methods. 2016 Jul;13(7):581-3. doi: 10.1038/nmeth.3869. Epub 2016 May 23. DADA2: High-resolution sample inference from Illumina amplicon data. https://www.ncbi.nlm.nih.gov/pubmed/27214047

  4. Deblur Rapidly Resolves Single-Nucleotide Community Sequence Patterns http://msystems.asm.org/content/2/2/e00191-16

  5. Python交互可视化图形库 http://bokeh.pydata.org/en/latest/

科研经验系列文章

1云笔记积累个人知识体系

2云协作建立实验室工作总结和内部资料共享平台

3公众号建立实验室共享知识体系和宣传窗口

扩增子图表解读第一季

1箱线图:Alpha多样性,老板再也不操心的我文献阅读了

2散点图:组间整体差异分析(Beta多样性)

3热图:差异菌、OTU及功能

4曼哈顿图:差异OTU或Taxonomy

5火山图:差异OTU数量及变化规律

6韦恩图:比较组间共有和特有OTU或分类单元

7三元图:一图展示三组两两比较

8网络图:节点OTU或类Venn比较

扩增子分析

QIIME2. 1简介和安装

QIIME. 1虚拟机安装配置及挂载外部目录

QIIME. 2.使用Docker安装并运行

QIIME. 3以管理员安装1.9.1至Ubuntu16.04

微生物相关网络构建教程中文:MENA, LSA, SparCC和 CoNet四种网络

想了解更多宏基因组、16S分析相关文章,

快关注“宏基因组”公众号,干货第一时间推送。

系统学习生物信息,快关注“生信宝典”,

那里有几千志同道合的小伙伴一起学习。




您可能也对以下帖子感兴趣

文章有问题?点此查看未经处理的缓存