基于16S的细菌群落功能预测工具-PICRUSt

PICRUSt（Phylogenetic Investigation of Communities by Reconstruction of Unobserved States）是一种使用标记基因数据和参考基因组数据库来预测宏基因组功能组成的计算方法，可基于16S rRNA图谱预测微生物群落的基因组功能组成（Langille et al, 2013）。

PICRUSt使用一种扩展的祖先状态重建算法来预测存在哪些基因家族，然后结合基因家族来估计复合宏基因组，包含一个两步过程。在最初的“基因含量推断”步骤中，通过参考系统发育树预先计算每个生物体的基因含量，这将重建一个基于16S的系统发育中每种生物的预测基因家族丰度表。随后的“宏基因组推断”步骤将对所有微生物分类群的所得基因含量预测与16S rRNA基因的拷贝数进行的校正，最终获得整个微生物群落预期的基因家族丰度。

接下来简介PICRUSt工具的使用。

考虑到在线版PICRUSt的网上教程貌似挺多的，但本地版PICRUSt的使用很少有教程提到；并且PICRUSt的输入文件（GreenGene注释的OTU表）也需本地操作，为了连贯性，所以本篇只介绍本地版PICRUSt工具。

PICRUSt简介及官方教程：http://picrust.github.io/picrust/index.html

下文所有示例数据、运行结果文件及命令行等，可在百度盘获取（提取码：sech）：

https://pan.baidu.com/s/1VsCThQWPddCvng02Pdmv3A

PICRUSt安装

conda安装PICRUSt环境

PICRUSt已经打包在bioconda里面了，可直接通过conda安装环境。

#使用 bioconda 安装 PICRUSt 环境
conda create -n picrust1 -c bioconda -c conda-forge picrust
 
#激活使用
source activate picrust1
#用完记得退出
source deactivate picrust1

手动下载需要的配置库文件

PICRUSt根据16S物种信息，预测群落功能。在这一过程中，需要读取相关的数据库内容信息，执行预测。PICRUSt提供了相关的库文件，根据需要下载。

链接：http://picrust.github.io/picrust/picrust_precalculated_files.html#id1

下载的文件是“.gz”压缩格式的，但无需解压。将下载好的文件放入PICRUSt环境安装位置的data路径下。

#比方说这里使用 conda 安装的 picrust 环境
#那么，data 目录就位于 conda 中的 picrust 环境目录下
 
#具体而言，比方说我的 conda 安装在
#~/software/Miniconda3
#那么本示例中，通过 conda 安装的 picrust 环境的 data 目录，即位于
#~/software/Miniconda3/envs/picrust1/lib/python2.7/site-packages/picrust/data/
#下载好的一系列的库文件，如“16S_13_5_precalculated.tab.gz”等，就移至这里即可

PICRUSt预测16S群落数据

软件安装并下载好配置文件后，就可以根据的16S物种数据预测群落功能了。

1、准备GreenGene注释的OTU丰度表

对于输入的OTU丰度表，必须为GreenGene数据库注释的。并且，PICRUSt识别的并不是常规OTU表中最后一列的注释信息，而是OTU的GreenGene ID。具体而言，就是这种形式。

如果不清楚，可以通过如下方法获得。

*通过qiime程序实现OTU序列的GreenGene注释

OTU表不是GreenGene注释的，就需要使用GreenGene进行注释。这里就以16S群落分析常用程序qiime为例进行注释。

示例文件可见网盘附件“data”。示例数据集共有80个16S测序样本，均来自土壤。因试验需求，在土壤中添加了某化学物质，目的为探究该化学物质对土壤微生物群落的影响。这80个样本中，40个为不添加化学物质的对照组（control组），40个为添加化学物质的处理组（treat组）。

文件“data/otu_table.txt”为OTU丰度表格，该OTU丰度表中无注释列，即所有OTU的物种分类是未知的，在后续将使用这些OTU的代表序列（见下文）鉴定这些OTU的物种分类。

文件“data/otu.fasta”中包含了OTU丰度表中各OTU的代表序列。在后续我们将使用OTU代表序列与GreenGene数据库中的参考物种序列进行比对，以鉴定这些OTU所属的“界门纲目科属种”分类单元。

首先根据获得的OTU的代表序列，通过GreenGene注释，并将注释结果添加到不带注释但已聚类好的OTU丰度表中。随后再根据注释概要，映射GreenGene ID，并用GreenGene ID替换原OTU表的OTU ID。

##OTU 序列注释
#（1）首先解压 gg_13_5_otus.tar.gz，这个是 GreenGene 数据库 13.5 版本
#（2）借助 qiime 程序实现 OTU 序列注释
 
#如果没安装 qiime，在 linux 中，通过 conda 安装 qiime 环境
#conda create -n qiime1 qiime
#激活使用
#source activate qiime1
#用完记得退出
#source deactivate qiime1
 
#激活 qiime 环境
#注释完成后，先不着急退出 qiime，后面的某些命令还需用它，最后再退
source activate qiime1
 
#qiime 的 OTU 序列注释过程
assign_taxonomy.py -i otu.fasta -r gg_13_5_otus/rep_set/97_otus.fasta -t gg_13_5_otus/taxonomy/97_otu_taxonomy.txt --uclust_max_accepts 1 -o gg97 
 
#根据具体的注释细节，对应 greengene id（这里用了一个手写的 R 脚本实现转换，在网盘附件中）
grep 'H' gg97/otu_tax_assignments.log > gg97/otu_tax_assignments2.log
Rscript replace_id.r gg97/otu_tax_assignments2.log otu_table.txt otu_table2.txt
 
#转为 biom 格式
sed -i '1i\# Constructed from biom file' otu_table2.txt
biom convert -i otu_table2.txt -o otu_table2.biom --table-type="OTU table" --to-json

qiime的输出可参考网盘附件“qiime_output”，“otu_table2.txt”即为最后获得的以GreenGene ID为注释的纯文本的OTU表；“qiime_output/otu_table2.biom”为转换为biom格式的OTU表，用于程序识别。

2、可选校正16S拷贝数（推荐）

对于很多细菌而言，一个个体可能包含多条16S（多拷贝16S）。PICRUSt提供了校正16S拷贝数的方法。

#激活 picrust 环境
#已经激活了 qiime，再激活 picrust，功能不冲突
source activate picrust1
 
#校正 16S 拷贝数
normalize_by_copy_number.py -i otu_table2.biom -o normalized_otus.biom
#转为 txt 表格
biom convert -i normalized_otus.biom -o normalized_otus.txt --to-tsv

可参考网盘附件“picrust_output/normalized_otus.txt”，即为校正16S拷贝数后的OTU丰度表；“picrust_output/normalized_otus.biom”为转换为biom格式的OTU表，用于程序识别。

3、群落功能预测

接下来，就使用校正16S拷贝数（物种丰度）后的OTU数据，执行功能预测。

#预测 KEGG 功能
predict_metagenomes.py -i normalized_otus.biom -o ko.biom -a nsti_ko.txt
#转为 txt 表格（分别为 KO 层级和 KO 功能描述）
biom convert -i ko.biom -o ko_pathways.txt --to-tsv --header-key KEGG_Pathways
biom convert -i ko.biom -o ko_description.txt --to-tsv --header-key KEGG_Description
 
#预测 COG 功能
predict_metagenomes.py -i normalized_otus.biom -o cog.biom -a nsti_cog.txt --type_of_prediction cog
#转为 txt 表格（分别为 COG 层级和 COG 功能描述）
biom convert -i cog.biom -o cog_category.txt --to-tsv --header-key COG_Category
biom convert -i cog.biom -o cog_description.txt --to-tsv --header-key COG_Description
 
#预测 Rfam 功能
predict_metagenomes.py -i normalized_otus.biom -o rfam.biom -a nsti_rfam.txt --type_of_prediction rfam
#转为 txt 表格
biom convert -i rfam.biom -o rfam_description.txt --to-tsv --header-key RFAM_Description

上述结果均可在网盘附件“picrust_output”中找到。

PICRUSt提供了对KEGG、COG、Rfam等数据库功能的映射。默认情况下，仅预测KEGG；其它数据库功能的预测，通过“--type_of_prediction”参数指定类别。

PICRUSt默认输入输出文件均为biom格式，可再转换为纯文本样式，文件结构类似OTU丰度表。以COG预测结果为例，第一列为预测所得COG功能ID，最后一列记录了该COG所示分类或功能描述，中间列为该功能的丰度信息。

获得群落功能丰度表后，就可以按照OTU丰度表的统计分析方法，去执行类似的分析了。这点可以找一些文献作参考，看别人是怎样做的。例如，首先计算特定功能丰度在组间的显著性，获得组间差异显著的功能，然后再从数据库官网上查询该功能的细节，解释生物学现象等。

4、PICRUSt的一些其它实用功能

除了获得整个群落水平的功能丰度外，PICRUSt还能提供在不同功能分类水平上的统计，以及单个OTU对功能的贡献。

#分解 KO，比方说在 KO 第 2 级或第 3 级功能上统计
categorize_by_function.py -f -i ko.biom -c KEGG_Pathways -l 3 -o ko_L3.txt
categorize_by_function.py -f -i ko.biom -c KEGG_Pathways -l 2 -o ko_L2.txt
 
#分解 COG，比方说在 COG 第 2 级功能上统计
categorize_by_function.py -f -i cog.biom -c COG_Category -l 2 -o cog_L2.txt
 
#OTU 对功能的贡献，需指定特定 KO 功能分类 id 或 COG 功能分类 id
metagenome_contributions.py -i normalized_otus.biom -l K00001,K00002,K00004 -o ko_metagenome_contributions.txt
metagenome_contributions.py -i normalized_otus.biom -l COG0001,COG0002 -t cog -o cog_metagenome_contributions.txt
 
#用完记得退出 qiime 或 picrust 环境哦
source deactivate picrust1
source deactivate qiime

上述结果均可在网盘附件“picrust_output”中找到。

biom文件转换为纯文本文件才可方便查看，以下继续以COG的预测结果为例展示。总之PICRUSt使用起来挺方便的，特别是它自己就含有直接在各功能分类层级上的统计功能，这样倒是可以省去再使用R等工具去作分类合并了。如果期望关注某些OTU是否对群落功能是重要的，贡献度表可以提供参考。

随后可能的统计分析方法，参考实际的文献中的方法来就可以了。

不妨作个图观察下组间功能分类的差异吧，还是以COG为例。

#R 语言操作
 
#读取 COG 二级分类统计
#把“cog_L2.txt”第一行的“# Constructed from biom file”删除再读取
cog <- read.delim('cog_L2.txt', sep = '\t', stringsAsFactors = FALSE, check.names = FALSE)
cog <- cog[-ncol(cog)]
 
#整理分类注释
cog$category <- apply(cog[1], 1, function(x) substr(x, 2, 2))
names(cog)[1] <- 'description'
cog <- reshape2::melt(cog, id = c('category', 'description'))
 
#合并样本分组文件
group <- read.delim('group.txt', sep = '\t', stringsAsFactors = FALSE, check.names = FALSE)
names(group)[1] <- 'variable'
cog <- merge(cog, group, by = 'variable')
 
#ggplot2 作图，细节就不调整了
library(ggplot2)
 
ggplot(cog, aes(x = category, y = value)) +
geom_boxplot(aes(color = group))

参考文献

Langille M G I, Zaneveld J, Caporaso J G, et al. Predictive functional profiling of microbial communities using 16S rRNA marker gene sequences. Nature Biotechnology, 2013, 31(9):814-821.

友情链接