扩增子分析解读5物种注释 OTU表操作

本节课程，需要先完成《扩增子分析解读》系列之前的操作
1质控实验设计双端序列合并
 2提取barcode 质控及样品拆分切除扩增引物
 3格式转换去冗余聚类
 4去嵌合体非细菌序列生成代表性序列和OTU表

分析前准备

# 进入工作目录
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.0
OE3: 2353.0
OE8: 3091.0
OE2: 3173.0
OE1: 3337.0
OE5: 3733.0
OE6: 4289.0
OE9: 4648.0
OE7: 5185.0
WT3: 10741.0
WT8: 12117.0
WT6: 14316.0
WT2: 14798.0
WT7: 14851.0
KO1: 14926.0
WT9: 15201.0
WT1: 15422.0
WT5: 15773.0
WT4: 16708.0
KO2: 17607.0
KO6: 23949.0
KO5: 26570.0
KO8: 27250.0
KO4: 32303.0
KO7: 33086.0
KO9: 35913.0
KO3: 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