lyao222lll的个人博客分享 http://blog.sciencenet.cn/u/lyao222lll

博文

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

已有 8079 次阅读 2019-2-16 18:24 |个人分类:R统计分析|系统分类:科研笔记| _R语言, _菌群功能预测

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

 

本文使用一个简单实例,介绍怎样使用RTax4Fun包,基于16S物种丰度数据预测细菌群落中代谢通路或功能基因的丰度,以探索群落功能的方法。在本文中,我们不使用Tax4Fun官网所提供的示例数据,而是使用我们自己的数据进行演示Tax4Fun功能预测流程。并在得到初步的预测结果之后,继续使用RSTAMP分别对预测结果进行简单的统计检验及作图展示,以提供一个生物学问题分析方法示例的简要补充。

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

 

我们此处共有8016S测序样本,均来自土壤。因试验需求,在土壤中添加了某化学物质,目的为探究该化学物质对土壤微生物群落的影响。这80个样本中,40个为不添加化学物质的对照组(control组),40个为添加化学物质的处理组(treat组)。

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

 

本文所涉及的数据已上传至百度盘,其中:

测试数据、脚本、简要结果等:https://pan.baidu.com/s/151AaSs61YNF7-zRXTFdnZg

SILVA123物种注释文件:https://pan.baidu.com/s/1OzJpZ5P25hgmcoBXrWsOnQ

KEGG代谢通路及功能基因注释文件:https://pan.baidu.com/s/18j8hZNTi9yCSDEpSeu-Tfw

 


示例文件简要


网盘附件“data”中,包含4个文件,内容分别如下。

 

文件“otu_table.txt”为OTU丰度表格,其内容展示如下。

每一列为一个样本,每一行为一种OTU,交叉区域为每种OTU在各样本中的丰度。该OTU丰度表中无注释列,即所有OTU的物种分类是未知的,在后续我们将使用这些OTU的代表序列(见下文)鉴定这些OTU的物种分类。

1.png

 

文件“group.txt”为样本分组信息,其内容展示如下。

第一列(names)为各样本名称;第二列(group)为各样本的分组信息,即这些样本分别属于未添加化学物质的对照组(control组)或添加了化学物质的处理组(treat组)。

2.png

 

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

3.png

 

文件“otu_table.not_sliva.tsv”同样为OTU丰度表格,每一列为一个样本,每一行为一种OTU,交叉区域为每种OTU在各样本中的丰度。

与文件“otu_table.txt”不同的地方在于,“otu_table.not_sliva.tsv”为qiime所得的标准输出样式,其第一行为“# Constructed from biom file”。此外,该文件的最后一列包含了各OTU的注释信息,但注释信息非SILVA数据库注释结果,且额外带有“p__”、“c__”、“g__”等分类单元前缀(这样无法被Tax4Fun所识别)。在后续我们将使用该文件,展示如何将这种注释样式转化为Tax4Fun可识别的注释样式。

4.png

 


获得Tax4Fun输入文件


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

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

 

使用qiimeOTU添加注释(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丰度表格。

备注:以下qiime步骤所得的中间结果文件,均放置在网盘附件“qiime_result”中,可点击查看。

##注:以下步骤均在 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

 

此外,在得到OTU注释文件“otu_table_tax_assignments.txt”之后,还可在R中,为无注释信息的OTU丰度表(otu_table.txt)添加注释信息。

除了使用 qiime,还可使用 R 读入无注释的 otu 表以及已经得到的 otu 注释文件(使用 qiime 注释得到),为 otu 表添加 sliva 注释
otu <- read.delim('otu_table.txt', sep = '\t', stringsAsFactors = FALSE, check.names = FALSE)
tax <- read.delim('otu_table_tax_assignments.txt', header = FALSE, sep = '\t', stringsAsFactors = FALSE, check.names = FALSE)[1:2]
names(otu)[1] <- 'otu_id'; names(tax) <- c('otu_id', 'taxonomy')
otu <- merge(otu, tax, by = 'otu_id')
write.table(otu, 'otu_table.silva.txt', row.names = FALSE, sep = '\t', quote = FALSE)
#继续使用的话,记得在“otu_table.silva.txt”第一行添加“# Constructed from biom file”

 

最终得到如下样式的OTU丰度表,其第一行为“# Constructed from biom file”;最后一列为OTU注释信息,且注释来源为SILVA数据库;其余行列即为各OTU名称、各样本名称、以及各OTU在各样本中的丰度分布。

1.png


SILVA123注释结果的转化(伪样式)


事实上,对于非SILVA注释的OTU结果,也可以使用Tax4Fun进行分析,只是实际效果会差那么一些。并且,Tax4Fun毕竟只是一种功能预测手段,其预测结果与真实结果肯定会存在差异,至于差异究竟多大也不好准确衡量。即便是使用SILVA注释的16S数据,也并不能完全代表真实的群落功能,一般也只是作为参考来使用。因此,若只是简单评估一下群落功能的话,也可使用非SILVA的注释结果来进行预测。

一般测序公司所提供的OTU注释结果是长这个样子的,即注释结果列中额外带有“p__”、“c__”、“g__”等分类单元前缀(参见示例文件“otu_table.not_sliva.tsv”)。在上述结果中可以看到,对于Tax4Fun能够识别的输入OTU表来讲,其OTU注释列那一栏中是不含有这些前缀的。若我们不在意测序公司是否使用了SILVA数据库进行OTU注释(尽管测序公司会在报告中给出),直接使用他们所提供现有的OTU表来进行Tax4Fun功能预测的话,首先需要将这种注释列中带分类单元前缀的OTU表进行转化(即去除注释前缀),得到Tax4Fun能够识别的样式。

读入非 sliva 注释的 otu 表,该表为 qiime 标准输出格式,需要将注释列转化为 Tax4Fun 能够识别的注释格式
otu <- read.delim('otu_table.not_sliva.txt', sep = '\t', stringsAsFactors = FALSE, check.names = FALSE)
otu$taxonomy <- gsub('[a-z]__', '', otu$taxonomy)
write.table(otu, 'otu_table.not_sliva.convert.txt', row.names = FALSE, sep = '\t', quote = FALSE)
#继续使用的话,记得在“otu_table.not_sliva.convert.txt”第一行添加“# Constructed from biom file”

 

最终得到如下样式的OTU丰度表,其第一行为“# Constructed from biom file”;最后一列为OTU注释信息,且注释内容中不再带有“p__”、“c__”、“g__”等分类单元前缀;其余行列即为各OTU名称、各样本名称、以及各OTU在各样本中的丰度分布。

1.png

 


Tax4Fun预测细菌群落功能(默认流程)


以下为Tax4Fun进行功能预测的默认流程,功能预测结果将得到KEGG功能基因及代谢通路的丰度。因此可以大致认为我们测了宏基因组,并进行了KEGG功能基因及代谢注释。

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

library(Tax4Fun)
 
##Tax4Fun 默认流程
##注:为了结果的准确可靠性,推荐使用 silva 数据库重新注释,尽量不要使用非 silva 数据库的结果
#此处使用“otu_table.silva.tsv”作为输入文件
QIIMESingleData <- importQIIMEData('otu_table.silva.tsv')
folderReferenceData <- 'SILVA123'
 
#功能基因丰度预测
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 代谢第三层级)丰度预测
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列为各功能基因登记ID(可根据此IDKEGG数据库中检索以查看其详细功能信息等)、名称及其对应的酶类型(酶命名规则可参考 https://enzyme.expasy.org/enzyme_details.html),第2列及之后为各样本中各功能基因的相对丰度。

对于代谢通路丰度预测结果“tax4fun.pathway.txt”(下图2),第1列为各KEGG代谢通路KO编号(可根据此IDKEGG数据库中检索以查看其详细功能信息等)及其名称,第2列及之后为各样本中各代谢通路的相对丰度。注意,此处的代谢通路为KEGG第三层级代谢通路结果。

5.png

6.png

 


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


以上即通过Tax4Fun功能预测得到了一种类似宏基因组测序分析中的KEGG代谢注释结果。那么,我们在得到了预测结果之后,应该怎样进行后续分析呢?

首先,若其中有我们感兴趣的代谢通路或者功能基因,或者我们想着重具体地了解某些代谢通路或者功能基因,可以进入KEGG官方网站(https://www.genome.jp/kegg/pathway.html),在网站中输入这些通路或者基因的ID进行查询,获知其详情。

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

0.png

1.png

2.png

3.png

 

或者,我们想要通过一些统计学手段,观测所预测出的代谢通路或功能基因丰度在各分组中是否具有显著差异,以帮助我们筛选出特定的信息出来以便后续有针对性地查阅相关文献等。下面给出了一个简单的示例作为参考。实际分析可能更为复杂,更需查阅大量文献中的方法作为参考,复杂类型的统计分析此处不再举例。

回顾本文一开始的问题,我们期望基于16S数据,得到群落的功能信息,并想得知经过试验处理后的细菌群落在功能上与对照组相比发生了怎样的改变。当然了,面对数量如此众多的功能基因,或者代谢通路,从中筛选出特殊功能的基因或代谢通路是件很不容易的事情,这里暂且不考虑。先考虑个简单的,我们从KEGG代谢通路高级层级(如第二层级,仅有几十个通路,便于统计)出发,关注这些代谢通路丰度在两个不同分组(处理组与对照组)的差异。

 

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


Tax4Fun预测出的KEGG代谢通路为第三层级结果。我们可以根据KEGG官网上所提供的代谢通路/功能基因分类注释文件,将低层级的通路(第三层级)归类到高层级的通路(第一、二层级)类型中,之后进行简要的统计分析。

在网盘附件中,提供了代谢通路/功能基因分类注释表(已在原列表文件的基础上转化为方便统计的表格样式,若有疑问可先查看附件中的说明)。通过该表,我们为一开始得到的预测结果“tax4fun.pathway.txt”添加注释。

##为 KEGG 第三层级代谢通路结果添加注释
#“pathway.anno.txt”为已经整理好的 KEGG 注释文件(可见网盘附件)
kegg_anno <- read.delim('pathway.anno.txt', sep = '\t', stringsAsFactors = FALSE, check.names = FALSE)
tax4fun_pathway$Pathway3_level <- t(data.frame(strsplit(rownames(tax4fun_pathway), ';')))[ ,1]
tax4fun_pathway <- merge(tax4fun_pathway, kegg_anno, by = 'Pathway3_level')
write.table(tax4fun_pathway, 'tax4fun.pathway.anno.txt', row.names = FALSE, sep = '\t', quote = FALSE)

 

这样,相较于一开始的输出结果(只包含KEGG第三层级代谢通路ID、名称及丰度),添加了这些代谢通路所属的更高级代谢层级的分类信息。

7.png

n列:各样本中各代谢通路的相对丰度信息,为KEGG第三层级代谢通路的统计结果;

Pathway1_levelKEGG代谢第一层级的通路编号;

Pathway1KEGG代谢第一层级的通路类型名称;

Pathway2_levelKEGG代谢第二层级的通路编号;

Pathway2KEGG代谢第二层级的通路类型名称;

Pathway3_levelKEGG代谢第三层级的通路编号;

Pathway3KEGG代谢第三层级的通路类型名称。

 

添加样本分组信息,并在KEGG代谢通路第二层级进行通路丰度的组间差异检验


由于Tax4Fun预测出的代谢通路丰度为KEGG第三层级代谢通路结果,将这些通路按其所对应的第二层级代谢分类求和后,即可得到KEGG第二层级各代谢通路在各样本中的丰度信息。之后根据分组文件为各样本添加分组信息。

library(reshape2)    #使用其中的 melt() 重排表格
library(doBy)  #使用其中的 summaryBy() 按分组(第二层级代谢通路)统计求和
 
#在 KEGG 第二层级统计求和,并添加样本分组信息
group <- read.delim('group.txt', sep = '\t', stringsAsFactors = FALSE, check.names = FALSE)
names(group)[1] <- 'variable'
tax4fun_pathway <- tax4fun_pathway[c(group$variable, 'Pathway2_level')]
tax4fun_pathway <- melt(tax4fun_pathway, id = 'Pathway2_level')
tax4fun_pathway <- summaryBy(value~variable+Pathway2_level, tax4fun_pathway, FUN = sum)
tax4fun_pathway <- merge(tax4fun_pathway, group, by = 'variable')

若此时查看结果“tax4fun_pathway”,其内容如下。

第一列为样本名称,第二列为KEGG第二层级各代谢通路的编号,第三列为KEGG第二层级各代谢通路在各样本中的相对丰度,第四列为各样本所属的分组信息。

8.png

 

然后我们继续使用doBy包中的summaryBy()命令,对第二层级各代谢通路类型的丰度根据样本分组求均值(mean)、标准差(sd)和标准误差(se)。此步骤用于得到作图文件。

#计算均值±标准差
se <- function(x) sd(x) / (length(x))^0.5
pathway_stat <- summaryBy(value.sum~group+Pathway2_level, tax4fun_pathway, FUN = c(mean, sd, se))

简单查看“pathway_stat”。第一列为两大样本分组(controltreat),第二列为KEGG第二层级各代谢通路的编号,第三、四、五列分别为各类型代谢通路在各样本分组中相对丰度发均值、标准差和标准误差。

9.png

 

以上结果中,对于所出现的KEGG第二层级各代谢通路,只包含了其编号。然后我们根据KEGG代谢通路注释文件,添加KEGG第二层级各代谢通路的名称,以及其所属第一(最高)层级通路分类。

#添加注释
kegg_anno_2 <- kegg_anno[!duplicated(kegg_anno$Pathway2), ][-c(5:7)]
pathway_stat <- merge(pathway_stat, kegg_anno_2, by = 'Pathway2_level', all.x = TRUE)
write.table(pathway_stat, 'tax4fun.pathway2_level.anno.txt', row.names = FALSE, sep = '\t', quote = FALSE)

输出文件“tax4fun.pathway2_level.anno.txt”中,与上述结果相比即添加了更详细的内容。

10.png

 

回顾我们所预期的问题,我们想观察经过试验处理后的细菌群落在功能上与对照组相比发生了怎样的改变。这里从KEGG代谢通路第二层级出发,关注这些代谢通路丰度在两个不同分组(处理组与对照组)的差异。

最容易想到的方法,也是最简单的方法,我们执行显著性检验。因为只有两个分组用于比较,因此我们考虑直接使用wilcoxon秩和检验。

#可选进行显著性差异分析,这里直接使用 wilcoxon 秩和检验
Pathway2_level <- unique(tax4fun_pathway$Pathway2_level)
for (i in Pathway2_level) {
    tax4fun_pathway_2_i <- subset(tax4fun_pathway, Pathway2_level == i)
    test <- wilcox.test(value.sum~group, tax4fun_pathway_2_i)
    line_t <- which(pathway_stat$Pathway2_level == 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'] <- '***'
    }
}

对于各类型的代谢通路,逐个执行wilcoxon秩和检验,比较每个类型的代谢通路在control组和treat组中是否存在相对丰度的差异。若存在擦怀疑,则根据显著性p值,在作图数据“pathway_stat”中添加显著性标记“*”以方便作图。

此时,可选查看“pathway_stat”。两两比较后(control组和treat组),我们将计算得到的显著性p值添加在二者之一的“treat”行中(“control”行不再添加重复的p值,R程序自动记为NA),同时在该行根据显著性添加“*”标记(同样地,“control”行不再添加)。其中,p0.05,不标记“*”;0.01p0.05,标记“*”;0.001p0.00,标记“**”;p0.001,标记“***”。

11.png

 

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


接下来,我们使用ggplot2进行结果的可视化,以直观表示本次的Tax4Fun功能预测及显著性差异分析结果。

library(ggplot2)
 
#使用 ggplot2 作图
pathway_stat$value.sum.mean <- 100 * pathway_stat$value.sum.mean
pathway_stat$value.sum.sd <- 100 * pathway_stat$value.sum.sd
 
pathway2_plot <- ggplot(pathway_stat, aes(Pathway2, value.sum.mean, fill = group)) +
geom_col(position = 'dodge', width = 0.8, colour = 'black', size = 0.05) +  #“dodge 柱状图”样式
geom_errorbar(aes(ymin = value.sum.mean - value.sum.sd, ymax = value.sum.mean + value.sum.sd), size = 0.05, width = .35, position = position_dodge(width = .8)) +  #添加误差线(均值±标准差)
scale_fill_manual(values = c('red', 'blue')) +  #填充颜色
theme(legend.title = element_blank(), legend.position = c(0.9, 0.95)) +  #去除图例标题,调整图例位置
coord_flip() +  #横、纵坐标轴反转
theme(panel.grid = element_blank(), panel.background = element_rect(fill = 'transparent',  color = 'black')) +  #去除默认的背景框
geom_text(aes(label = sign, y = value.sum.mean + value.sum.sd + 0.5), size = 4, position = position_dodge(0.8)) +  #添加显著性标记“*”
labs(x = 'KEGG pathway2', y = 'Relative Abundance (%)')  #坐标轴标题
 
ggsave('Kegg_pathway2.pdf', pathway2_plot, width = 8, height = 10)
ggsave('Kegg_pathway2.png', pathway2_plot, width = 8, height = 10)


可视化展示结果中,各类型代谢通路在各分组中多峰丰度信息一目了然,并且我们也很直观的知道了经过试验处理后的细菌群落在功能上与对照组相比发生了怎样的改变。

注:此处仅简单展示了统计KEGG第二层级代谢通路丰度在两大分组中的差异的思路和方法,其它的更深入、更有针对性的分析还需参阅更多的文献。

12.png

对于R的原始出图,我们可以在后期使用AIPS等加以调整,如多图组合等,以得到效果更好的图片,如下所示。

注:下图与本次的示例数据毫无关系,仅作为效果图展示(只因为本人实在懒得再对上图进行图操作了……直接使用了一张现有的其它的图)。

13.png

 

使用STAMP进行组间差异分析及作图


这类有关于代谢通路差异的分析,我们经常可以在文献中见到使用STAMP进行的分析结果,效果蛮不错的。因此接下来作为补充,展示STAMP进行这类代谢途径差异分析的流程,以及结果的可视化展示。

同样地,我们仍然在KEGG第二层级代谢通路水平上统计各类型丰度在两大分组中的差异。

 

首先需要得到STAMP输入文件,这里我们从得到Tax4Fun最初的预测结果开始,一步步处理得到。

library(Tax4Fun)
library(reshape2)
library(doBy)
 
##得到 STAMP 输入文件
#此处使用“otu_table.silva.tsv”作为输入文件
QIIMESingleData <- importQIIMEData('otu_table.silva.tsv')
folderReferenceData <- 'SILVA123'
 
#代谢通路(KEGG 代谢第三层级)丰度预测
Tax4FunOutput <- Tax4Fun(QIIMESingleData, folderReferenceData, fctProfiling = FALSE, refProfile = 'UProC', shortReadMode = TRUE, normCopyNo = TRUE)
tax4fun_pathway <- as.data.frame(t(Tax4FunOutput$Tax4FunProfile))
 
#添加 KEGG 代谢第二层级注释
kegg_anno <- read.delim('pathway.anno.txt', sep = '\t', stringsAsFactors = FALSE, check.names = FALSE)[c(4, 5)]
tax4fun_pathway$Pathway3_level <- t(data.frame(strsplit(rownames(tax4fun_pathway), ';')))[ ,1]
tax4fun_pathway <- merge(tax4fun_pathway, kegg_anno, by = 'Pathway3_level')[-1]
 
#在 KEGG 第二层级统计求和
tax4fun_pathway <- melt(tax4fun_pathway, id = 'Pathway2')
tax4fun_pathway <- summaryBy(value~variable+Pathway2, tax4fun_pathway, FUN = sum)
 
#转化为 STAMP 输入文件格式
tax4fun_pathway <- dcast(tax4fun_pathway, Pathway2~variable)
write.table(tax4fun_pathway, 'stamp.pathway2_level.txt', row.names = FALSE, sep = '\t', quote = FALSE)

输出结果文件“stamp.pathway2_level.txt”即为STAMP的输入文件。其内容如下。

第一列为KEGG第二层级各代谢通路名称,其后列为各样本中所包含各代谢通路的相对丰度。

14.png

 

然后我们运行STAMP,点击“File > Load data”将KEGG第二层级代谢通路丰度表格文件“stamp.pathway2_level.txt”以及样本分组信息文件“group.txt”导入至STAMP中。数据导入后,我们点击“View > Group legend”调出分组图例选项查看,这里只有两个分组control组和treat组。

15.png

之后在界面中选择显著性检验方式,并设置参数。软件自动统计分析并绘制图形。

可在图片区域下方选择图形展示样式,我们选择带误差棒的柱状图样式(STAMP特色作图类型)。

16.png

在图片下方的“Configure plot”选项中,可以调整作图参数,修改作图样式。例如,误差线展示类型、图片尺寸、文字大小等。

17.png

点击“View > Two group statistics table”,可查看统计检验的详细计算结果表格信息。包含了KEGG第二层级各代谢通路丰度的均值、标准差、差异检验p值等内容。若有需要,可点击“Save”将统计结果保存为表格样式。

18.png

此外,可以通过点击“File > Save plot...”选项,将作图结果导出。矢量图如pdf格式,位图如png格式。STAMP导出的图片也可以后期使用AIPS等进行调整,以得到更好的效果。下图为原始出图。

19.png

 


参考文献


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.




http://blog.sciencenet.cn/blog-3406804-1162514.html

上一篇:基因组分析简介之K-mer分析
下一篇:R包mixOmics进行微生物群落偏最小二乘判别分析(PLS-DA)

0

该博文允许注册用户评论 请点击登录 评论 (8 个评论)

数据加载中...
扫一扫,分享此博文

Archiver|手机版|科学网 ( 京ICP备07017567号-12 )

GMT+8, 2020-8-12 18:52

Powered by ScienceNet.cn

Copyright © 2007- 中国科学报社

返回顶部