本节课程，需要先完成《扩增子分析解读》系列之前的操作
分析前准备
# 进入工作目录
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查看一下吧，内容如下：
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