本节课程，需要先完成《扩增子分析解读》系列之前的操作
分析前准备
# 进入工作目录 cd example_PE250
上一节回顾：我们学习了嵌合体的形成，以及基于参考数据库去嵌合体；也学习了基于数据库比对来筛选细菌或真菌；最后基于最确定的OTU，我们生成代表性序列和OTU表，这是每种高通量测序都有的结果，后续的结果将全部基于这两个文件。
接下来我们学习对OTU进行物种注释；OTU的操作，包括格式转换、筛选添加物种信息、数据量筛选样品、筛选高丰度的OTU、物种筛选等。
OTU表常用的BIOM格式
主页：http://biom-format.org/ 。BIOM是英文The Biological Observation Matrix的缩写，中文翻译为生物观测矩阵，是一种通过格式，用于生物学样品对应观测值的表格。它主要采用json/HD5F文件格式标准，即多维散列结构，保存表格结构数据结果。目前主流的宏基因组软件均支持此格式文件，如QIIME、MG-RAST、PICRUSt、Mothur、phyloseq、MEGAN、VAMPS、metagenomeSeq、Phinch、RDP Classifier、USEARCH、PhyloToAST、EBI Metagenomics、GCModeller、MetaPhlAn 2。知道它有多重要了吧。
Biom文件处理系统biom程序是QIIME的必装包，如果没有安装好，可尝试下面步骤重装
# 安装依赖包 pip install numpy # 安装biom格式转换包 pip install biom-format # 安装2.0格式支持 pip install h5py # 测序程序是否安装成功 biom
13. 物种注释
对于扩增子分析，最重要的就是物种信息。我们基于上节分析得到的代表性序列，采用上次已经下载的greengene的参考序列和物种注释信息，比对软件选择rdp方法，进行注释。
# 物种注释 assign_taxonomy.py -i result/rep_seqs.fa -r gg_13_8_otus/rep_set/97_otus.fasta -t gg_13_8_otus/taxonomy/97_otu_taxonomy.txt -m rdp -o result
注：如果是ITS/18S数据，建议数据库更改为UNITE，方法改为blast。详细使用说明，请读官方文档http://qiime.org/scripts/assign_taxonomy.html
14. OTU表统计、格式转换、添加信息
将OTU表转换为Biom格式，这样便于其它软件对其操作。可添加上面获得的物种信息，这样表格的信息就更丰富了，再转换为文本，便于人类可读，同时使用summarize-table查看OTU表的基本信息。
# 文本OTU表转换为BIOM：方便操作 biom convert -i temp/otu_table.txt -o result/otu_table.biom --table-type="OTU table" --to-json # 添加物种信息至OTU表最后一列，命名为taxonomy biom add-metadata -i result/otu_table.biom --observation-metadata-fp result/rep_seqs_tax_assignments.txt -o result/otu_table_tax.biom --sc-separated taxonomy --observation-header OTUID,taxonomy # 转换biom为txt格式，带有物种注释：人类可读 biom convert -i result/otu_table_tax.biom -o result/otu_table_tax.txt --to-tsv --header-key taxonomy # 查看OTU表的基本信息：样品，OUT数量统计 biom summarize-table -i result/otu_table_tax.biom -o result/otu_table_tax.sum
现在我们获得了OTU表的基本统计信息，用less result/otu_table_tax.sum查看一下吧，内容如下：
Num samples: 27 # 样品数据Num observations: 975 # OTU数据Total count: 409647 # 总数据量Table density (fraction of non-zero values): 0.464 # 非零的单元格Counts/sample summary:Min: 2352.0 # 样品数据量最小值Max: 35955.0 # 样品数据量最大值Median: 14851.000 # 样品数据量中位数Mean: 15172.111 # 样品数据量平均数Std. dev.: 10691.823 # 样品数据量标准变异Sample Metadata Categories: None provided # 样品分类信息：末提供Observation Metadata Categories: taxonomy # 观察值分类：物种信息Counts/sample detail: # 每个样品的数据量OE4: 2352.0OE3: 2353.0OE8: 3091.0OE2: 3173.0OE1: 3337.0OE5: 3733.0OE6: 4289.0OE9: 4648.0OE7: 5185.0WT3: 10741.0WT8: 12117.0WT6: 14316.0WT2: 14798.0WT7: 14851.0KO1: 14926.0WT9: 15201.0WT1: 15422.0WT5: 15773.0WT4: 16708.0KO2: 17607.0KO6: 23949.0KO5: 26570.0KO8: 27250.0KO4: 32303.0KO7: 33086.0KO9: 35913.0KO3: 35955.0
biom的详细使用说明，可以biom查看具体的功能，如添加注释功能biom add-metadata --help可查看详细说明。也可阅读官网http://biom-format.org/
15. OTU表筛选
实验中会有各种影响因素，我们要综合各种背景知识来判断如何筛选数据表，起到去伪存真，去粗取粗，由此及彼，有表及理的来回答科学问题。数据筛选是会运行分析流程和数据分析师的分水岭。
看上面的的统计结果，样本数据量从2k-35k，我们应去除过小的数据量样本，提供更可能高的样品最低丰度的数据用于下游标准化分析。这里我们选择只保留数据量大于3000的样品。
# 按样品数据量过滤：选择counts>3000的样品 filter_samples_from_otu_table.py -i result/otu_table_tax.biom -o result/otu_table2.biom -n 3000 # 查看过滤后结果：只有25个样品，975个OTU biom summarize-table -i result/otu_table2.biom
同时还要过滤低丰度的OTU，一般低于万分之一丰度的菌，在功能研究可能还是比较困难的(早期文章454测序数据量少，通常只关注丰度千分之五以上的OTU)。
# 按OTU丰度过滤：选择相对丰度均值大于万分之一的OTU filter_otus_from_otu_table.py --min_count_fraction 0.0001 -i result/otu_table2.biom -o result/otu_table3.biom # 查看过滤后结果：只有25个样品，346个OTU biom summarize-table -i result/otu_table3.biom
有些研究手段在特定有实验中存在偏差，如2012Nature报导V5-V7在植物中扩增会偏好扩增Chloroflexi菌门，建议去除。
# 按物种筛选OTU表：去除p__Chloroflexi菌门 filter_taxa_from_otu_table.py -i result/otu_table3.biom -o result/otu_table4.biom -n p__Chloroflexi # 查看过滤后结果：只有25个样品，307个OTU biom summarize-table -i result/otu_table4.biom
以上过滤条件是根据经验、相关文献设计的，如果不清楚，也不要随便过滤，容易引起假阴性。
得到的最终结果，还要转换为文本格式，和提取OTU表对应的序列，用于下游分析。
# 转换最终biom格式OTU表为文本OTU表格 biom convert -i result/otu_table4.biom -o result/otu_table4.txt --table-type="OTU table" --to-tsv # OTU表格式调整方便R读取 sed -i '/# Const/d;s/#OTU //g;s/ID.//g' result/otu_table4.txt # 筛选最终OTU表中对应的OTU序列 filter_fasta.py -f result/rep_seqs.fa -b result/otu_table4.biom -o result/rep_seqs4.fa