分析前准备
# 进入工作目录 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=" ";FS=" "} NR==FNR {a[$1]=$0} NR>FNR {print a[$1]}' result/rep_seqs_tax_assignments.txt temp/tax_rep_seqs_clus.id|sed 's/; / /g'|cut -f 1-5 |sed 's/p__//g;s/c__//g;s/o__//g' > result/tax_rep_seqs.tax
21. 其它
其它都是一些简单的格式转换，为后面统计分析而准备文件。
# 将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/;/ /g;s/ //g' result/rep_seqs_tax_assignments.txt > result/rep_seqs_tax.txt