使用ReporterScore包进行富集分析

最新推荐文章于 2024-05-22 00:00:00 发布
最新推荐文章于 2024-05-22 00:00:00 发布
阅读量678 收藏 6
点赞数
分类专栏： R语言 功能富集分析 文章标签： r语言 数据分析 信息可视化
版权声明：本文为博主原创文章，遵循 CC 4.0 BY-SA 版权协议，转载请附上原文出处链接和本声明。
本文链接：https://blog.csdn.net/Asa12138/article/details/132216910
版权

Introduction

上次已经在一篇推文中介绍过了微生物组分析常用的功能富集分析方法以及reporter score方法的原理，并且也介绍了我开发的R包ReporterScore。

但时那个时候R包写的还比较粗糙，功能也不多，最近进一步优化了这个包的各种功能，支持两种模式，多种统计检验方法，支持两组甚至更多组的实验设计（这个挺好用的），KEGG数据库的同步做的也比较好了，还增加了一些可视化形式。

以下是我给这个R包github主页（https://github.com/Asa12138/ReporterScore）下写的介绍和用法，这里简单翻译一下贴过来了。但其实里面还有不少功能没在下面写出，可以在R包里探索一下，哈哈。

Citation

这个包暂时还没有在刊物上发表，要在出版物中引用 ReporterScore，请使用：

Chen Peng, Chao Jiang. ReporterScore: an R package for Reporter Score-based Microbial Enrichment Analysis. R package (2023), https://github.com/Asa12138/ReporterScore

🤗也欢迎在GitHub上点个star⭐️

Install

if(!require("devtools"))install.packages("devtools")
devtools::install_github('Asa12138/pcutils',dependencies=T)
devtools::install_github('Asa12138/ReporterScore',dependencies=T)
library(ReporterScore)

Usage

1. Inputdata (KO abundance table and metadata)

通常，我们可以使用KEGG数据库来注释我们的微生物组测序数据，特别是环境微生物组，因为KEGG相对来说更全面（当然大部分还是unknown）。

具体方法包括直接使用blast对KEGG序列数据库进行比对，使用KEGG官方mapper软件，使用EggNOG数据库并将结果转化为KEGG注释。

这样我们就可以得到一个KO丰度表（行是KO，列是样本）用于我们的富集分析：

data("KO_abundance_test")
head(KO_abundance[,1:6])

##                WT1         WT2         WT3         WT4         WT5         WT6
## K03169 0.002653545 0.005096380 0.002033923 0.000722349 0.003468322 0.001483028
## K07133 0.000308237 0.000280458 0.000596527 0.000859854 0.000308719 0.000878098
## K03088 0.002147068 0.002030742 0.003797459 0.004161979 0.002076596 0.003091182
## K03530 0.003788366 0.000239298 0.000445817 0.000557271 0.000222969 0.000529624
## K06147 0.000785654 0.001213630 0.001312569 0.001662740 0.002387006 0.001725797
## K05349 0.001816325 0.002813642 0.003274701 0.001089906 0.002371921 0.001795214

还需要提供实验设计的元数据metadata（行是样本，列是组）。

head(metadata)

##     Group Group2
## WT1    WT     G3
## WT2    WT     G3
## WT3    WT     G3
## WT4    WT     G3
## WT5    WT     G3
## WT6    WT     G1

2. Pathway database

ReporterScore内置了KEGG 通路和模块数据库（2023-08 版）用于KO 丰度表分析。

你可以使用 load_KOlist() 查看并使用 update_KO_file() 更新这些数据库（通过 KEGG API），因为使用最新的数据库非常重要。

或者你可以使用custom_modulelist()自定义你自己的通路数据库（感兴趣的基因集）。

load_KOlist()
head(KOlist$pathway)

3. One step enrichment

使用函数reporter_score可以一步得到reporter score结果。

有一些重要的参数可供调节：

  • mode: “mixed” 或 “directed”（仅适用于两组差异分析或多组相关分析），详细信息参见pvalue2zs。
  • 方法：统计检验类型。 默认为”wilcox.test”：
    • t.test （参数）和 wilcox.test （非参数）。 对两组样品进行比较。 如果分组变量包含两个以上水平，则执行成对比较。 - anova（参数）和 kruskal.test（非参数）。 执行比较多个组的单向方差分析测试。
    • “pearson”、“kendall”或”spearman”（相关分析），请参见cor。
  • p.adjust.method：用于测试结果的p.adjust.method，参见p.adjust。
  • type/modulelist：选择通路数据库，默认数据库为”pathway”或”module”，或使用自定义的模块列表。

group作为因子变量，第一个水平将被设置为对照组，你可以更改因子水平来改变你的比较。

例如，我们想要比较两组”WT-OE”，并使用”directed”模式，因为我们只想知道 OE 组 中哪些通路被富集或耗尽：

cat("Comparison: ",levels(factor(metadata$Group)))
## Comparison:  WT OE

reporter_score_res=reporter_score(KO_abundance,"Group",metadata,mode="directed")

结果是一个”reporter_score”对象：

elementsdescription
kodf你的输入 KO_abundance 表
ko_pvalueko 统计结果包含 p.value
ko_statko统计结果包含p.value和z_score
reporter_s在每个途径中的reporter score
modulelist默认 KOlist 或自定义模块列表数据框
group你的数据中的比较组
metadata示例信息数据框包含组

4. Visualization

绘制最显着富集的通路：

#View(reporter_score_res$reporter_s)
plot_report(reporter_score_res,rs_threshold = c(-2,2))

当我们聚焦于一条通路时，例如 “map00780”：

plot_KOs_in_pathway(reporter_score_res,map_id = "map00780")

或者显示为网络：

plot_KOs_network(reporter_score_res,map_id = "map00780",main="")

我们也可以查看通路中每个 KO 的丰度：

plot_KOs_box(reporter_score_res,map_id = "map00780",only_sig = TRUE)

或者热图形式呈现：

plot_KOs_heatmap(reporter_score_res,map_id = "map00780",only_sig = TRUE,heatmap_param = list(cutree_rows=2))

example for “mixed”

如果我们的实验设计超过两组，我们可以选择多组比较和“mixed”模式：

cat("Comparison: ",levels(factor(metadata$Group2)))
## Comparison:  G1 G2 G3

reporter_score_res2=reporter_score(KO_abundance,"Group2",metadata,mode="mixed",
      method = "kruskal.test",p.adjust.method1 = "none")

plot_KOs_in_pathway(reporter_score_res2,map_id = "map00541")

Details

Step by step

一步函数 reporter_score 由三部分组成：

  1. ko.test：此函数有助于通过各种内置方法计算 KO_abundance 的 p-value，例如差异分析（t.test、wilcox.test、kruskal.test、anova）或相关分析（pearson 、spearman、kendall）。 你还可以通过其他方法计算 KO_abundance 的 p-value，例如“DESeq2”、“Edger”、“Limma”、“ALDEX”、“ANCOM”，并自行进行 p值矫正，然后跳过ko.test 步骤转到步骤2…
  2. pvalue2zs：该函数将 KO 的 p-value 转换为 Z-score（选择模式: “mixed” 或 “directed”）。
  3. get_reporter_score 该函数计算特定数据库中每个通路的reporter score。 你可以在此处使用自定义数据库。

这样你就可以一步一步得到reporter score。

Other enrichment methods

ReporterScore 还提供了其他丰富方法，如 KO_fisher(fisher.test)、KO_enrich(fisher.test, from clusterProfiler)、KO_gsea (GSEA, from clusterProfiler)，输入数据来自 reporter_score，并且也支持自定义数据库，因此你可以轻松比较各种富集方法的结果并进行全面分析：

data("KO_abundance_test")
reporter_score_res2=reporter_score(KO_abundance,"Group",metadata,mode="mixed")
#View(reporter_score_res2$reporter_s)
#reporter_score
reporter_score_res2$reporter_s$p.adjust=p.adjust(reporter_score_res2$reporter_s$p.value,"BH")
filter(reporter_score_res2$reporter_s,(ReporterScore)>1.64,p.adjust<0.05)%>%pull(ID)->RS
#fisher
ko_pvalue=reporter_score_res2$ko_pvalue
fisher_res=KO_fisher(ko_pvalue)
filter(fisher_res,p.adjust<0.05)%>%pull(ID)->Fisher
#enricher
enrich_res=KO_enrich(ko_pvalue)
filter(enrich_res,p.adjust<0.05)%>%pull(ID)->clusterProfiler
#GESA
set.seed(1234)
gsea_res=KO_gsea(ko_pvalue)
filter(gsea_res@result,p.adjust<0.05)%>%pull(ID)->GSEA

venn_res=list(RS=RS,Fisher=Fisher,CP=clusterProfiler,GSEA=GSEA)
library(pcutils)
venn(venn_res)

venn(venn_res,"network",vertex.label.cex=c(rep(1,4),rep(0.5,22)))

KO levels

KEGG BRITE 是一个层次分类系统的集合，捕获各种生物对象的功能层次结构，特别是那些表示为 KEGG 对象的功能层次结构。

我们收集了 k00001 KEGG Orthology (KO) 表，以便你可以总结每个级别的丰度。 使用 load_KO_htable 获取 KO_htable 并使用 update_KO_htable 进行更新。 使用 up_level_KO 可以升级到“pathway”、“module”、“level1”、“level2”、“level3”、“module1”、“module2”、“module3”之一中的特定级别。

load_KO_htable()
head(KO_htable)

## # A tibble: 6 × 8
##   level1_id level1_name level2_id level2_name        level3_id level3_name KO_id
##   <chr>     <chr>       <chr>     <chr>              <chr>     <chr>       <chr>
## 1 A09100    Metabolism  B09101    Carbohydrate meta… map00010  Glycolysis… K008…
## 2 A09100    Metabolism  B09101    Carbohydrate meta… map00010  Glycolysis… K124…
## 3 A09100    Metabolism  B09101    Carbohydrate meta… map00010  Glycolysis… K008…
## 4 A09100    Metabolism  B09101    Carbohydrate meta… map00010  Glycolysis… K250…
## 5 A09100    Metabolism  B09101    Carbohydrate meta… map00010  Glycolysis… K018…
## 6 A09100    Metabolism  B09101    Carbohydrate meta… map00010  Glycolysis… K068…
## # ℹ 1 more variable: KO_name <chr>

plot_KO_htable()

KO_level1=up_level_KO(KO_abundance,level = "level1",show_name = TRUE)
pcutils::stackplot(KO_level1[-which(rownames(KO_level1)=="Unknown"),])

Reference

  1. Patil, K. R.
    & Nielsen, J. Uncovering transcriptional regulation of metabolism by using metabolic network topology.
    Proc Natl Acad Sci U S A 102, 2685–2689 (2005).

  2. L. Liu, R. Zhu, D. Wu, Misuse of reporter score in microbial enrichment analysis. iMeta. 2, e95 (2023).

  3. https://github.com/wangpeng407/ReporterScore

关注公众号，获取最新推送

关注公众号 ‘biollbug’,获取最新推送。

Asa12138
关注
  • 0
    点赞
  • 踩
  • 6
    收藏
    觉得还不错? 一键收藏
  • 4
    评论
专栏目录
Reporter Score 微生物功能富集分析
Asa12138的博客
05-15 1024
Reporter score是一种改良的微生物富集分析的新方法，这里分享其原理和一个实现的R包。
GSEAmining | 来看看你的GSEA结果是不是需要瘦身啦！~
m0_72224305的博客
02-23 300
1写在前面 最近真是累的不行，今天抽空写一下新的教程，关于人人都会做的GSEA（Gene Set Enrichment Analysis）。 但有时候我们做完GSEA后结果实在太多，无法确定其中重要的生物学意义，难以解释。🤨 本期我们介绍一下GSEAmining包，对我们的GSEA结果做一个瘦身吧，基本原理是：👇 1️⃣ 对参与类似生物过程的基因集应该有共同的基因。 2️⃣ 对拥有一定数量的共同基因的相似基因集进行功能聚类。 2用到的包 rm(list = ls())# if (!requireName
r语言进行go富集分析_GO富集柱状图
weixin_39923262的博客
12-20 3542
分析模块，输入差异基因GO富集分析结果，由分析模块“GO Enrichment Analysis”生成，输出GO富集结果柱状图。输入：差异基因GO富集分析结果文件，由分析模块“GO Enrichment Analysis”生成。示例：id enrichment description ratio_in_study ratio_in_pop p_un...
python实现富集分析（Enrichment Analysis）
qq_39579290的博客
09-19 4700
python实现KEGG富集分析（Enrichment Analysis）
功能的显著性分析——GO Enrichment Analysis
weixin_30300523的博客
12-15 4710
Gene Ontology（GO）是基因功能国际标准分类体系。GO富集分析是对差异基因等按GO分类，并对分类结果进行基于离散分布的显著性分析、错判率分析、富集度分析，得到与实验目的有显著联系的、低误判率的、靶向性的基因功能分类，该分类即导致样本性状差异的最重要的功能差别。在芯片的数据分析中，研究者可以找出哪些变化基因属于一个共同的GO功能分支，并用统计学方法检定结果是否具有统计学意义，从而...
ReporterScore 包的一些更新
Asa12138的博客
11-09 117
最近ReporterScore 包又新加了一些好用的功能
植物MWAS研究—小米产量与微生物组关联分析
刘永鑫的博客——宏基因组公众号
09-10 7188
MWAS简介微生物组关联分析(Microbiome/Metagenome-wide association studies , MWAS)是指捕获多维尺度上的互作作用，从而提供捕获复杂作用关系的方法，该方法切实可行的预测微生物组和疾病状态的关系。就是不做任何假设的分析整个微生物组的结构变化，目的是鉴定与疾病、表型相关的微生物物种、基因或代谢物质等。全微生物组关联分析 (MWAS) 与全基因组关联分析
iMeta | 浙大吴顶峰和同济刘蕾等对微生物富集分析中Reporter Score的误用研究
悟道西方
03-05 387
点击蓝字 关注我们微生物富集分析中Reporter Score的误用iMeta主页：http://www.imeta.science研究论文●原文链接DOI: https://doi.org/10.1002/imt2.95●2023年2月23日，浙江大学吴顶峰和同济大学刘蕾等在iMeta在线发表了题为“Misuse of reporter score in microbial enric...
生信学习之通路富集一（GO分析）
热门推荐
coffeeii的博客
04-25 2万+
生信学习之通路富集（一）
Cell文章最喜欢用的差异基因GO、KEGG富集分析工具
qq_42090739的博客
01-04 7536
对于转录组分析，差异基因筛选完成后，就算是成功一大半了。接下来就可进行下游分析了。首先最常见的就是对上下调基因的富集分析了，富集分析我们之前出过一期R语言版本的： clusterProfiler：基因功能富集分析的惊喜之作 这里就不再赘述了，我们介绍一款NCS文章比较喜欢的、出现率较高的功能富集工具---Metascape (https://metascape.org/gp/index.html)。 图片 Metascape是一款免费的在线工具，使用简便易操作。更重要的是其包含多个数据库，结果可
GSEA - Gene set enrichment analysis 基因集富集 | ORA - Over-Representation Analysis 分析原理与应用...
weixin_30294709的博客
02-13 2326
R批量做GSEA分析还没有官方的包，但是clusterprofiler可以做，它调用了最新的gfsea包。Gene Set Testing for RNA-seq - fgsea教程 RNA-seq是利器，大部分做实验的老板手下都有大量转录组数据，所以RNA-seq的分析需求应该是很大的（大部分的生信从业人员应该都差不多要沾边吧）。 普通的转录组套路并不多，差异表达基因、富集分析、WGC...
r包调用legend函数_你还在愁毕业？隔壁实验室的小哥从网上抄了几十行代码打了个R包，发了SCI，毕业了！...
weixin_30120827的博客
12-25 231
你是否对《R包辣鸡之CorMut》之篇还有印象？没有的话，快点打开阅读一下，这文章让作者毕业了，毕竟发表的是一篇Bioinformatics。《[连载3]：辣眼睛，一篇抄袭引发的系列血案！》，而这一篇文章中揭露了某讲师抄袭了两个R包，晋升副教授了！文章中还稍带了另外一个学生也是抄袭，当然也发表了SCI，也毕业了。今天讲另外一个R包，它做下面的事情：当然前面的事情是goseq干的，它的功能...
proliferfer各种富集分析大杂烩 各种多种数据集来源的富集分析
生信小博士的博客
04-16 304
proliferfer各种富集分析大杂烩 各种多种数据集来源的富集分析 12 Universal enrichment analysis The clusterProfiler package (Yu et al. 2012) supports both > hypergeometric test and gene set enrichment analyses of many ontology/pathway, > but it’s still not enough for users may
R语言lavaan结构方程模型（SEM）
weixin_xiao5kou4chang6kai4的博客
05-21 388
结构方程模型（Sructural Equation Modeling，SEM）是分析系统内变量间的相互关系的利器，可通过图形化方式清晰展示系统中多变量因果关系网，具有强大的数据分析功能和广泛的适用性，是近年来生态、进化、环境、地学、医学、社会、经济等众多领域应用十分广泛的统计方法。在R语言结构方程程序包中，lavaan具有简洁的语法结构、成熟模型构建和调整过程和稳定可靠的结果等特点，使其不亚于收费商业软件，是最受欢迎的结构方程模型程序包之一。1）R及Rstudio介绍：背景、软件及程序包安装、基本设置等。
C++ RBTree
m0_74922218的博客
05-21 543
C++实现 RBTree，没有网传的那么恐怖
R语言数据分析案例框架
Dalao_zzl的博客
05-21 176
。假设我们是一家电商公司的数据分析师，公司最近推出了一系列促销活动，我们希望通过分析销售数据来评估这些活动的效果。 请注意，以上只是一个简化的案例框架和示例代码片段，你需要根据具体的数据和业务需求来扩展和完善这个案例。
【启明智显技术分享】SOM2D02-2GW核心板适配ALSA（适用Sigmastar ssd201/202D）
ami82的博客
05-21 93
本文针对sigmastat ssd202d/SSD201方案适配ALSA的方法介绍
R实验 参数估计
最新发布
Jiangxia13的博客
05-22 809
实验目的： 掌握矩法估计与极大似然估计的求法； 了解估计量的优良性准则：无偏性、有效性、相合性（一致性）； 学会利用R软件完成一个正态总体均值和两个正态总体均值差的区间估计； 学会利用R软件完成两个成对数据均值差的区间估计； 学会利用R软件完成一个总体比例和两个总体比例差的区间估计； 掌握大样本数据关于单个总体均值和总体比例的样本容量的确定方法。
R语言：多因素cox
weixin_69558614的博客
05-14 805
head(rt)#cox模型建立head(cox)

“相关推荐”对你有帮助么？

  • 非常没帮助
  • 没帮助
  • 一般
  • 有帮助
  • 非常有帮助
提交
评论 4
添加红包

请填写红包祝福语或标题

个

红包个数最小为10个

元

红包金额最低5元

当前余额3.43元 前往充值 >
需支付：10.00元
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付元
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明：

1.余额是钱包充值的虚拟货币，按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载，可以购买VIP、付费专栏及课程。

余额充值