扩增子分析解读7物种分类统计,筛选进化树和其它

写在前面

之前发布的《扩增子图表解读》系列，相信关注过我的朋友大部分都看过了(链接直达7月文章目录)。这些内容的最初是写本实验室的学生们学习的材料，加速大家对同行文章的解读能力。

《扩增子分析解读》系列文章介绍

扩增子分析是目前宏基因组研究中最常用的技术，由于微生物组受环境影响大，实验间重复较差，更需要更多的实验重复和分析技术来保证结果的准确性、可重复性。

本系统文章叫分析解读，即有详细的扩增子分析流程代码，又有本人对使用参数、备选参数意义的解读，可以让大部分零基础的人，更好的理解数据分析过程，并可亲自实践在自己的课题上，获得更好、更合理的实验结果。

本文采用目前最主流的扩增子测序数据类型HiSeq2500 PE250类型数据为例，结合目前主流方法QIIME+USearch优点组合定制的分析流程。本课程中所需的测序数据、实验设计和课程分析生成的中间文件，均可以直去百度云下载。链接：http://pan.baidu.com/s/1hs1PXcw 密码：y33d。

本课程代码的运行，至少需要Linux平台+安装QIIME1.9.1，我之前发布过QIIME1.9.1安装的三种方法：1虚拟机、 2Docker、 3管理员，总有一款适合你。

声明：本套流程主要依赖QIIME1。之前发布的QIIME2不是QIIME的升级版，而是完全独立的分析系统，两者没有任何通用的地方，而且现在还不成熟，明年才有稳定版。请读者千万别混淆。不要再犯用QIIME2系统运行本套扩增子分析解读，无法找到相关程序的错误。

第七节. 物种分类统计,筛选进化树和其它

本节课程，需要完成扩增子分析解读之前的任务：
1质控,实验设计,双端序列合并
2提取barcode,质控及样品拆分,切除扩增引物
3格式转换,去冗余,聚类
4去嵌合体,非细菌序列,生成代表性序列和OTU表
4去嵌合体,非细菌序列,生成代表性序列和OTU表
5物种注释,OTU表操作
6进化树,Alpha,Beta多样性

先看一下扩增子分析的整体流程，从下向上逐层分析。

分析前准备

# 进入工作目录
cd example_PE250

上一节回顾：我们获得了OTU序列的进化分析、同时计算Alpha和Beta多样性值。

本节是最后一节，我们对物种进行分类统计,筛选高丰度结果用于进化树展示，和其它用于R统计分析的结果生成

19. 按物种分类级别分类汇总

OTU表中最重要的注释信息是物种注释信息。通常的物种注释信息分为7个级别：界、门、纲、目、科、属、种。种是最小的级别，和OTU类似但有不相同。
我们除了可以比较样品和组间OTU水平差异外，还可以研究不同类似级别上的差异，它们是否存在那些共同的变化规律。

按照注释的级别进行分类汇总，无论是Excel还R操作起来，都是很麻烦的过程。这里我们使用QIIME自带 的脚本summarize_taxa.py。

# 结果按门、纲、目、科、属五个级别进行分类汇总，对应结果的L2-L6
summarize_taxa.py -i result/otu_table4.biom -o result/sum_taxa # summary each level percentage
# 修改一下文本表头，适合R读取的表格格式
sed -i '/# Const/d;s/#OTU ID.//g' result/sum_taxa/* # format for R read
# 以门为例查看结果
less -S result/sum_taxa/otu_table4_L2.tx

以门为例，我们看到样品的OTU分布在19个门，及每个门在各样品中的相对比例。其它的各级别，用户自己看吧。

这步的结果将用于后期统计和绘图。

20. 筛选可展示的进化树

我们在文章中看到几种漂亮的进化树，但是OTU通常成百上千，如果直接展示是根本看不清也是极丑的。
下面教大家一些通常的方江来筛选数据，用于用成漂亮的进化树。

# 选择OTU表中丰度大于0.1%的OTU
filter_otus_from_otu_table.py --min_count_fraction 0.001 -i result/otu_table4.biom -o temp/otu_table_k1.biom
# 获得对应的fasta序列
filter_fasta.py -f result/rep_seqs.fa -o temp/tax_rep_seqs.fa -b temp/otu_table_k1.biom 
# 统计序列数量，104条，一般100条左右即有大数据的B格，又能读懂和更清规律和细节
grep -c '>' temp/tax_rep_seqs.fa # 104
# 多序列比对
clustalo -i temp/tax_rep_seqs.fa -o temp/tax_rep_seqs_clus.fa --seqtype=DNA --full --force --threads=30
# 建树
make_phylogeny.py -i temp/tax_rep_seqs_clus.fa -o temp/tax_rep_seqs.tree
# 格式转换为R ggtree可用的树
sed "s/'//g" temp/tax_rep_seqs.tree > result/tax_rep_seqs.tree # remove '
# 获得序列ID
grep '>' temp/tax_rep_seqs_clus.fa|sed 's/>//g' > temp/tax_rep_seqs_clus.id
# 获得这些序列的物种注释，用于树上着色显示不同分类信息
awk 'BEGIN{OFS="\t";FS="\t"} NR==FNR {a[$1]=$0} NR>FNR {print a[$1]}' result/rep_seqs_tax_assignments.txt temp/tax_rep_seqs_clus.id|sed 's/; /\t/g'|cut -f 1-5 |sed 's/p__//g;s/c__//g;s/o__//g' > result/tax_rep_seqs.tax

21. 其它

其它都是一些简单的格式转换，为后其统计分析而准备文件。以后再碰到什么问题，在这里补充，在CSDN的文章我可以修改，并把大家反馈的问题和答案放在这部分。

# 将mappingfile转换为R可读的实验设计
sed 's/#//' mappingfile.txt > result/design.txt
# 转换文本otu_table格式为R可读
sed '/# Const/d;s/#OTU //g;s/ID.//g' result/otu_table4.txt > result/otu_table.txt
# 转换物种注释信息为制表符分隔，方便R读取
sed 's/;/\t/g;s/ //g' result/rep_seqs_tax_assignments.txt > result/rep_seqs_tax.txt

写在后面

在Shell中分析的部分就到这里了，其实QIIME也可出一些图，但是由于结果难看且不容易修改，我很少用。接下来的教程，我会基于分析的结果在R语言中进行统计分析，并绘制出各种各样的图，大家期待吧。要想了解这些程序的更多功能，一定要阅读程序的帮助全文，才能有更深入的理解。

(宏基因组7月文章目录，更多精彩等你读)

Reference

1. http://qiime.org/scripts/summarize_taxa.html
2. http://qiime.org/scripts/filter_otus_from_otu_table.html
3. http://qiime.org/scripts/make_phylogeny.html

想了解更多16S/ITS/18S扩增子、宏基因组、宏转录组文献阅读和分析相关文章，快关注“宏基因组”公众号，干货第一时间推送。

系统学习生物信息，快关注“生信宝典”，那里有几千志同道合的小伙伴一起学习。