做通路富集分析的时候究竟输入多少个基因？

[ixxmu]

https://mp.weixin.qq.com/s/4cV8R4NxNklW4jIzsqLGbg

[ixxmu]

做通路富集分析的时候究竟输入多少个基因？ by 果子学生信

我在做差异基因富集分析的时候常有两个疑问？
第一，上调基因和下调基因，需要分开做么，分开做有优势么？
第二，富集分析输入的差异基因究竟是多少个？

关于第一个问题，我们已经解决了
把compareCluster富集后缺失的组显示出来
现在来尝试解决第二个问题，我们的方法还是试。
我们将会用不同数量的基因来做富集分析，然后把结果可视化，看看有什么区别，还是会用到神奇函数compareCluster

启发是Y叔的例子

在Y叔最新发表的clusterProfiler的文章中，Figure4十分亮眼

他优秀的展示了一个细胞在不同时间点，接受两种不同药物处理后下游通路的变化。
我们能看到有一些通路是药物特异的，有一些下游通路随着时间推移是药物共有的。
这样的图，十分容易俘获基础科研人员的心，我觉得是生信和基础结合的典范。

原文中提供了该图分析的代码，并且随着文章一起公布了。

data(DE_GSE8057)
xx <- compareCluster(Gene∼time+treatment, 
data=DE_GSE8057, fun = enricher,
TERM2GENE=wp[,c("wpid", "gene")], 
TERM2NAME=wp[,c("wpid", "name")])

接下来我们就来复现这个图，一起了解下他需要的数据格式。

library(clusterProfiler)
load(file = "data/DE_GSE8057.rda")
## downloaded from https://wikipathways-data.wmcloud.org/current/gmt/
gmt <- 'data/wikipathways-20210710-gmt-Homo_sapiens.gmt'
wp <- read.gmt.wp(gmt)

先看DE_GSE8057这个数据
这个数据是个数据框，有三列，分别是Gene，time，treatment
应该原始数据中有2个药物6个时间点，总共12组数据，通过跟0h比较可以得到每一组的差异基因。
(具体的需要看这个数据的原文)

另外一个数据就是基因集，这里下载的是Wikipathways的数据
这个数据有5列

接下来就是使用compareCluster来进行富集分析
总共就三类参数
第一是geneClusters和data，是输入的数据以及分组信息
Gene肯定是必要的，此外需要用~指示分组信息，用+好来连接多个分组信息
第二是enricher这个通用的ORA富集函数
第三是enricher是其他参数，就是基因集的信息

xx <- compareCluster(geneClusters=Gene~time+treatment, 
                     data=DE_GSE8057, 
                     fun = enricher,
                     TERM2GENE=wp[,c("wpid", "gene")], 
                     TERM2NAME=wp[,c("wpid", "name")])

做完了就是直接画图

dotplot(xx, x="time")

但是我们发现，横坐标的顺序，不是我们希望的，这里面可以通过调整因子的levels水平来控制
跟之前显示缺失的组那里是一种手法
调整后就可以作图了

## 富集分析
xx <- compareCluster(geneClusters=Gene~time+treatment, 
                     data=DE_GSE8057, 
                     fun = enricher,
                     TERM2GENE=wp[,c("wpid", "gene")], 
                     TERM2NAME=wp[,c("wpid", "name")])
## 调整因子水平
xx@compareClusterResult$time = factor(xx@compareClusterResult$time,
                                      levels =c("0h","2h","6h","24h"))
## 作图
dotplot(xx, x="time")

现在这个顺序就是正确舒适的。

此时已经能够展示各个时间点特异的通路了，但是我们还可以用分面，把结果分成两个药物来展示
代码也很简单

dotplot(xx, x="time") + ggplot2::facet_grid(~treatment)

只不过效果出众，原文的图就出现了。

如果我们想要分开的是时间来展示，那也可以的，
把代码稍微调整一下就行。

dotplot(xx, x="treatment") + ggplot2::facet_grid(~time)+
  scale_x_discrete(guide = guide_axis(angle = 45)) 

如果你注意到Y轴的那些term名称有重叠，那么可以通过增加label_format参数来在一行显示
这个自动折叠功能很有用，尤其是wikipathways这样的条目名称很长的基因集

dotplot(xx, x="time",label_format = 60) + ggplot2::facet_grid(~treatment)

富集分析放入多少基因合适？

在搞清楚以上作图所需要的数据后，我们就可以来操作了。
我们只要输入不同数目的基因来做富集分析，并且把结果放在一起展示就行。
跟之前一样，把差异分析的结果加载进来，区分上调和下调基因，不使用LogFC筛选

rm(list = ls())
library(clusterProfiler)
library(dplyr)
library(tibble)
load(file = "output/DEseq2_ESR1_Diff.Rdata")
diffgene <- res %>% 
  filter(gene !="") %>% 
  filter(adj.P.Val < 0.05) %>% 
  arrange(desc(logFC))

upgeneall <- diffgene$gene[diffgene$logFC>0]
dngeneall <- diffgene$gene[diffgene$logFC<0]

为了简化操作，先写一个函数，可以输入需要的基因数目，然后返回一个准备好的数据

generate_genes <- function(number){
  ## 选取一定数目上调基因
  upgene = head(upgeneall,number)
  ## 选取一定数目下调基因
  dngene = tail(dngeneall,number)
  ## 输入文件
  gcSample = data.frame(group= c(rep("up",number),rep("down",number),rep("all",number*2)),
                        genes = c(upgene,dngene,c(upgene,dngene)),
                        num = number)
  return(gcSample)
}

测试一下函数是够能够运行

test <- generate_genes(500)

会返回一个数据库，也有三列，group里面是up，down 和all

接下来批量产生不同数目的数据,此处就以500,1000，1500, 2000为例

mydata <- do.call("rbind",lapply(c(500,1000,1500,2000),generate_genes))

正式富集分析,代码也很好理解啦

h_df <- read.gmt("resource/h.all.v2022.1.Hs.symbols.gmt")
xx <- compareCluster(genes~group+num, 
                     data=mydata, 
                     fun="enricher",
                     TERM2GENE = h_df)

为了使得最后的结果有正确的顺序，我们再次使用因子调整一下顺序

xx@compareClusterResult$num <- factor(xx@compareClusterResult$num,
                                      levels = c("500","1000","1500","2000"))
xx@compareClusterResult$group <- factor(xx@compareClusterResult$group,
                                      levels = c("up","down","all"))

激动人心的时刻来了，作图代码类似，增加了一些参数

dotplot(xx,
        showCategory = 30,
        label_format = 60,
        font.size = 10,
        x = "group")+ 
  facet_grid(~num)

我们能够看到，有一些通路，比如雌激素相关的通路，无论你输入多少个基因，他都是能富集到的。
这些通路的可靠性很高，并且我们引入了上调基因和下调基因，还能看出来，这两个通路是被抑制的。
而这个与该数据的背景十分吻合，因为他就是ESR1敲减的数据。

同时我们能看到，基因的数目在1000个往后，通路多了起来，有些通路无论基因如何增加
都很坚挺的存在，他们也应该是我们关注的通路。这个例子说，你应该尝试用不同数量的基因去做分析。
尝试，并不花钱，但是可以让你更加了解这个数据。
这样的尝试性质的图，他可以很直白的告诉别人，我选择的通路是可信的。
此时换一种表达形式，也十分方便

dotplot(xx,
        showCategory = 30,
        label_format = 60,
        font.size = 10,
        x = "num")+ 
  facet_grid(~group)+
  scale_x_discrete(guide = guide_axis(angle = 45)) 

这时候，激活，抑制，各个基因数目的影响，看的更直观了。

当然你可以尝试更多的基因数目，但是要注意版面不一定放得下。
此刻这里面红框的通路，只在特定情况出现，他可能可以救你一命，也可能害你三年，所以慎重选择。

总结升华

对于差异分析的结果，这里的富集分析用的是over-representation analysis
这个通路富集分析的第一代技术。
后面为了解决这种不知道输入多少基因，不知道究竟是激活还是抑制的问题，
第二代通路富集技术产生了，他就是我十分喜欢的基因集富集分析(GSEA)，我们还做过专门的专题。

同时，在本文中我演示了如何学习并把技能化为己用的过程。他包括两个步骤。
第一，要首先复现出阳性的例子，搞清楚一个函数需要什么数据
第二，要有能力把自己的数据调整到他需要的样子
这样，所有的R包都能自由学习了。

最后感谢Y叔写出了clusterPofiler这样神奇的R包，让ORA和GSEA这样的富集分析和可视化变得简单易用。
当前通路富集分析已经发展到了第四代，这些新的技术不够普及的最大原因，
可能就是没有出现clusterPofiler这样简单易用的R包，能够让基础科研的工作者能够简单的上手用起来。