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

Tax4Fun基于16S物种丰度预测细菌群落代谢功能及统计作图的一个简单实例

Tax4Fun是一个R包，可根据16S rRNA数据集预测微生物群落的功能（Aßhauer et, 2015）。Tax4Fun基于最小16S rRNA序列相似性的最近邻识别实现16S rRNA基因序列与原核生物基因组功能注释的连接。OTU丰度表必须来自SILVA特定版本的注释，随后Tax4Fun将OTU物种丰度转换为功能配置文件，该过程分为三个步骤。

首先，将基于SILVA的16S rRNA图谱转换为原核生物KEGG的分类学图谱，线性变换是通过预先计算的关联矩阵实现的。然后，根据NCBI基因组注释获得的16S rRNA拷贝数对KEGG生物的估计丰度进行归一化处理。最后，利用归一化的分类学丰度，将预测得到的KEGG生物的功能图谱线性组合，预测微生物群落的功能图谱。 

本篇使用一个简单示例，介绍使用R的Tax4Fun包，基于16S物种丰度数据预测细菌群落中代谢通路或功能基因的丰度，并在得到初步的预测结果之后，继续使用R或STAMP分别对预测结果进行简单的统计检验及作图展示，以探索群落功能的方法。

了解Tax4Fun更多详情请参见官网（Tax4Fun包也需在其官网中下载获得）：http://tax4fun.gobics.de/

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

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

示例文件

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

此处期望基于16S数据，根据群落物种组成使用Tax4Fun预测群落的功能，并想得知经过试验处理后的细菌群落在功能上与对照组相比发生了怎样的改变。 

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

 文件“data/group.txt”为样本分组信息，第一列（names）为各样本名称；第二列（group）为各样本的分组信息，即这些样本分别属于未添加化学物质的对照组（control组）或添加了化学物质的处理组（treat组）。

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

获得Tax4Fun输入文件

因Tax4Fun要求输入的OTU丰度表中，各OTU必须来自SILVA数据库的注释；并且Tax4Fun目前所具有3个参考版本，分别对应了3个SILVA物种注释库版本，最好二者相互对应。因此，为了预测结果的准确性，首先需要保证我们的OTU注释信息来自于SILVA数据库的注释。若最初的OTU注释信息非SILVA注释（如来自Greengenes、RDP等），此时需要使用SILVA数据库重新注释。SILVA注释时，最好使用与所将使用Tax4Fun参考版本所对应的SILVA物种注释数据集版本。

例如，以下分析示例中挑选了Tax4Fun的SILVA123版本进行细菌群落功能预测，因此此处最好使用SILVA数据库中的SILVA123物种注释库对我们的OTU进行物种注释。  

使用qiime为OTU添加注释（SILVA123物种注释）

对于非SILVA注释的OTU注释信息，可以使用这些OTU的代表序列重新进行注释。这里使用到无注释信息的OTU丰度表（otu_table.txt）以及包含各OTU的代表序列文件（otu.fasta），使用16S数据常用分析工具qiime来完成。

使用各OTU的代表序列，将它们与SILVA123数据库中的代表物种序列进行比对，找到每个OTU序列的最相似序列，以期得到每个OTU所对应的物种分类单元（OTU注释信息）。同时将OTU丰度表转化为便于qiime识别的biom样式，并使用qiime将各OTU注释信息添加至OTU表的最后一列。最终得到具有SILVA123物种注释信息的OTU丰度表格。

##注：以下步骤均在 shell 命令行中完成，运行程序或脚本均来自 qiime
#在 otu_table.tsv 开头添加一行“# Constructed from biom file”，以便将 otu_table.tsv 转为 qiime 可识别样式
cp otu_table.txt otu_table.tsv
sed -i '1i\# Constructed from biom file' otu_table.tsv
 
#otu_table.tsv 转换为 otu_table.biom
biom convert -i otu_table.tsv -o otu_table.biom --table-type="OTU table" --to-json
 
#OTU 注释，输出结果 otu_table_tax_assignments.txt 即为注释文件
assign_taxonomy.py -i otu.fasta -r SILVA_123_SSURef_Nr99_tax_silva.fasta -t SILVA_123_SSURef_Nr99_tax_silva.tax -o ./
 
#添加 otu 注释信息至 biom 格式的 otu 表（otu_table.biom ）的最后一列，并将列名称命名为 taxonomy
biom add-metadata -i otu_table.biom --observation-metadata-fp otu_table_tax_assignments.txt -o otu_table.silva.biom --sc-separated taxonomy --observation-header OTUID,taxonomy 
 
#otu_table.silva.biom 转换为 otu_table.silva.tsv
biom convert -i otu_table.silva.biom -o otu_table.silva.tsv --to-tsv --header-key taxonomy

网盘附件“qiime_result/otu_table.silva.tsv”为最后获得的OTU丰度表，其第一行为“# Constructed from biom file”；最后一列为OTU注释信息，且注释来源为SILVA数据库；其余行列即为各OTU名称、各样本名称、以及各OTU在各样本中的丰度分布。

Tax4Fun预测细菌群落功能（默认流程）

以下为Tax4Fun进行功能预测的默认流程，功能预测结果将得到KEGG功能基因及代谢通路的丰度，并进行了KEGG功能基因及代谢注释。

输入了OTU丰度表，包含了每个测序样本中物种类型组成及其丰度信息。Tax4Fun会根据这些信息，在其数据库中找到对应的物种类别，根据16S拷贝数将物种数据标准化后，统计所对应物种类别其本身所具有的功能基因即所含代谢通路的数量，预测细菌群落中所具有的功能基因丰度以及代谢通路丰度。因此可知，OTU注释结果中，若注释信息越详细（分类水平越细，越接近种水平），则理论上功能注释结果的准确度越高。

library(Tax4Fun)
 
##Tax4Fun 默认流程
#此处使用“otu_table.silva.tsv”作为输入文件
QIIMESingleData <- importQIIMEData('otu_table.silva.tsv')
folderReferenceData <- 'SILVA123'
 
#KEGG 功能基因（KO 第四级）丰度预测
Tax4FunOutput <- Tax4Fun(QIIMESingleData, folderReferenceData, fctProfiling = TRUE, refProfile = 'UProC', shortReadMode = TRUE, normCopyNo = TRUE)
tax4fun_gene <- as.data.frame(t(Tax4FunOutput$Tax4FunProfile))
gene <- rownames(tax4fun_gene)
write.table(cbind(gene, tax4fun_gene), 'tax4fun.gene.txt', row.names = FALSE, sep = '\t', quote = FALSE)
 
#KEGG 代谢通路（KO 第三级）丰度预测
Tax4FunOutput <- Tax4Fun(QIIMESingleData, folderReferenceData, fctProfiling = FALSE, refProfile = 'UProC', shortReadMode = TRUE, normCopyNo = TRUE)
tax4fun_pathway <- as.data.frame(t(Tax4FunOutput$Tax4FunProfile))
pathway <- rownames(tax4fun_pathway)
write.table(cbind(pathway, tax4fun_pathway), 'tax4fun.pathway.txt', row.names = FALSE, sep = '\t', quote = FALSE)

两个流程里面，一个重要参数“fctProfiling”可分别设定为“TRUE”和“FALSE”。当“fctProfiling = TRUE”时，将输出功能基因丰度预测结果；当“fctProfiling = FALSE”时，将输出代谢通路预测结果。其他参数一般情况下无需更改，默认即可。

对于功能基因丰度预测结果“tax4fun.gene.txt”（下图1），第1列为各功能基因KO ID（可根据此ID在KEGG数据库中检索以查看其详细功能信息等）、名称及其对应的酶类型（酶命名规则可参考https://enzyme.expasy.org/enzyme_details.html），第2列及之后为各样本中各功能基因的相对丰度。

对于代谢通路预测结果“tax4fun.pathway.txt”（下图2），第1列为各KEGG代谢通路KO ID（可根据此ID在KEGG数据库中检索以查看其详细功能信息等）及其名称，第2列及之后为各样本中各代谢通路的相对丰度。

以上即通过Tax4Fun功能预测得到了一种类似宏基因组测序分析中的KEGG代谢注释结果。若其中有我们感兴趣的代谢通路或者功能基因，或者想着重具体地了解某些代谢通路或者功能基因，可以进入KEGG官方网站（https://www.genome.jp/kegg/pathway.html），在网站中输入这些通路或者基因的ID进行查询，获知其详情。

例如，查询代谢通路“Glycolysis / Gluconeogenesis”了解详情，即“00010”，过程和结果如下所示。

根据预测结果进行统计分析（一个简单的示例）

以及通过一些统计分析，观测所预测出的代谢通路或功能基因丰度在各分组中是否具有显著差异，以便后续有针对性地查阅相关文献等。

例如以下展示一个小示例。

为 KEGG代谢通路结果添加分类注释

例如，将KEGG代谢通路映射到更高级的KO（KEGG Orthology）分类单元（如第一、二级），这样，相较于一开始的输出结果（只包含KEGG代谢通路名称），添加了这些代谢通路所属的更高级的分类信息。

#为 KEGG 代谢通路映射 KO 分类
#“pathway.anno.txt”为已经整理好的 KEGG 注释文件（可见网盘附件）
kegg_anno <- read.delim('pathway.anno.txt', sep = '\t', colClasses = 'character', check.names = FALSE)
tax4fun_pathway$KO3_id <- t(data.frame(strsplit(rownames(tax4fun_pathway), ';')))[ ,1]
tax4fun_pathway$KO3_id <- gsub('ko', '', tax4fun_pathway$KO3_id)
tax4fun_pathway <- merge(tax4fun_pathway, kegg_anno, by = 'KO3_id')
write.table(tax4fun_pathway, 'tax4fun.pathway.anno.txt', row.names = FALSE, sep = '\t', quote = FALSE)

KO1_id：KO第一级分类的ID；

KO1：KO第一级分类的名称；

KO2_id：KO第二级分类的ID；

KO2：KO第二级分类的名称；

KO3_id：KO第三级分类的ID；

KO3：KO第三级分类的名称，也即具体的KEGG代谢通路；

其余列：预测得到的丰度信息。

KO二级分类单元的组间差异检验

继续关注这些不同的功能在两个不同分组（处理组与对照组）的差异。

以下不妨以KO第二级分类为例展示，计算映射到的丰度，并根据分组信息计算差异。

#在KO 第二级统计求和，并添加样本分组信息
group <- read.delim('group.txt', sep = '\t', stringsAsFactors = FALSE, check.names = FALSE)
names(group)[1] <- 'variable'
tax4fun_pathway <- tax4fun_pathway[c(group$variable, 'KO2_id')]
tax4fun_pathway <- reshape2::melt(tax4fun_pathway, id = 'KO2_id')
tax4fun_pathway <- doBy::summaryBy(value~variable+KO2_id, tax4fun_pathway, FUN = sum)
tax4fun_pathway <- merge(tax4fun_pathway, group, by = 'variable')
 
#计算均值（mean）±标准差（sd），或标准误差（se）
se <- function(x) sd(x) / (length(x))^0.5
pathway_stat <- doBy::summaryBy(value.sum~group+KO2_id, tax4fun_pathway, FUN = c(mean, sd, se))
 
#添加注释
kegg_anno_2 <- kegg_anno[!duplicated(kegg_anno$KO2), ][-c(5:7)]
pathway_stat <- merge(pathway_stat, kegg_anno_2, by = 'KO2_id', all.x = TRUE)
#write.table(pathway_stat, 'tax4fun.KO2.anno.txt', row.names = FALSE, sep = '\t', quote = FALSE)
 
#显著性差异分析，例如这里直接使用非参数的 wilcoxon 检验两组间差异
KO2_id <- unique(tax4fun_pathway$KO2_id)
for (i in KO2_id) {
        tax4fun_pathway_2_i <- subset(tax4fun_pathway, KO2_id == i)
        test <- wilcox.test(value.sum~group, tax4fun_pathway_2_i)
        line_t <- which(pathway_stat$KO2_id == i & pathway_stat$group == 'treat')
        pathway_stat[line_t,'p_value'] <- test$p.value
        if (test$p.value < 0.05 & test$p.value >= 0.01) {
                pathway_stat[line_t,'sign'] <- '*'
        }
        if (test$p.value < 0.01 & test$p.value >= 0.001) {
                pathway_stat[line_t,'sign'] <- '**'
        }
        if (test$p.value < 0.001) {
                pathway_stat[line_t,'sign'] <- '***'
        }
}
 
#write.table(pathway_stat, 'tax4fun.KO2.anno_stat.txt', row.names = FALSE, sep = '\t', quote = FALSE, na = '')

“tax4fun.KO2.anno_stat.txt”结果中包含了各KO功能分类单元的名称，在各组中丰度的均值、标准差、标准误差，以及wilcoxon检验的显著性p值等信息。

ggplot2作图展示功能预测及组间差异分析结果

接下来，使用ggplot2进行结果的可视化，展示本次Tax4Fun功能预测及显著性差异分析结果。

#使用 ggplot2 作图
library(ggplot2)

#给 KO 名称按丰度数值大小排个序
pathway_stat <- pathway_stat[order(pathway_stat$value.sum.mean), ]
pathway_stat$KO2 <- factor(pathway_stat$KO2, levels = unique(pathway_stat$KO2))
pathway_stat$KO1 <- factor(pathway_stat$KO1, levels = rev(unique(pathway_stat$KO1)))

#作图
pathway_stat$value.sum.mean <- 100 * pathway_stat$value.sum.mean
pathway_stat$value.sum.sd <- 100 * pathway_stat$value.sum.sd

ko2_plot <- ggplot(pathway_stat, aes(x = KO2, y = value.sum.mean, fill = group)) +
geom_col(position = 'dodge', width = 0.8, colour = 'black', size = 0.05) + #“dodge 柱状图”样式
scale_fill_manual(values = c('red', 'blue')) + #填充颜色
geom_errorbar(aes(ymin = value.sum.mean - value.sum.sd, ymax = value.sum.mean + value.sum.sd),
size = 0.05, width = 0.35, position = position_dodge(width = 0.8)) + #添加误差线（均值±标准差）
geom_text(aes(label = sign, y = value.sum.mean + value.sum.sd + 0.5), size = 4, position = position_dodge(0.8)) + #添加显著性标记“*”
facet_grid(KO1~., space = 'free', scale = 'free_y') + #分面图，按 KO1 分类
coord_flip() + #横、纵坐标轴反转
theme(panel.grid = element_blank(), panel.background = element_rect(fill = 'transparent', color = 'black'),
legend.title = element_blank(), legend.position = 'top') + #背景、图例主题等设置
labs(x = 'KEGG Orthology (KO)', y = 'Relative Abundance (%)') #坐标轴标题

ggsave('ko2_stat.pdf', ko2_plot, width = 8, height = 10)
ggsave('ko2_stat.png', ko2_plot, width = 8, height = 10)

参考文献

Aßhauer, Kathrin P, Wemheuer B , Daniel R , et al. Tax4Fun: predicting functional profiles from metagenomic 16S rRNA data. Bioinformatics, 2015, 31(17):2882-2884.

友情链接