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

本网对Markdown排版支持较差，请跳转“宏基因组”公众号阅读；

写在前面

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

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

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

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

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

本课程代码的运行，至少需要Linux平台+安装QIIME1.9.1，我之前发布过三种安装QIIME的方法详见文章目录，总有一款适合你。

第五节. 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程序 python是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

写在后面

今天先到这里，本文已经讲了太多了，够大家学习一会的了。要想了解这些程序的更多功能，一定要阅读程序的帮助全文，才能有更深入的理解。

下节预告：进化树构建、Alpha和Beta多样性。

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

Reference

http://biom-format.org/
http://qiime.org/scripts/assign_taxonomy.html
http://qiime.org/scripts/filter_samples_from_otu_table.html
http://qiime.org/scripts/filter_taxa_from_otu_table.html
http://qiime.org/scripts/filter_fasta.html

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

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