你有没有想过差异基因做富集分析的时候，究竟该输入什么基因？
全体？上调？下调？
如果搞不清的时候就都做一下，看看有什么差别，然后就知道以后该怎么做了。
正好Y叔的compareCluster 函数就可以做这个，简直是神器。
rm(list = ls())
library(clusterProfiler)
library(dplyr)
library(tibble)
load(file = 'output/DEseq2_ESR1_Diff.Rdata')
数据的名称是res，数据的内容就是Deseq2差异分析的结果，也可以是limma芯片的结果
只要有基因和排序的标准就行，这里我们筛选p值小于0.05的，基因变化倍数不做要求，只用它来排序
diffgene <- res %>%
filter(gene !='') %>%
filter(adj.P.Val < 0.05) %>%
arrange(desc(logFC))
现在留下的都是p值小于0.05，按照logFC从高到低排序的表格了
分别提取升高的基因和降低的基因
upgeneall <- diffgene$gene[diffgene$logFC>0]
dngeneall <- diffgene$gene[diffgene$logFC<0]
因为这些基因被排过序，升高的就是按照从高到底排序
降低的基因，logfc是负的，他的绝对值logFC是越来越大的
分别选取，需要输入的差异基因
## 选取一定数目上调基因
upgene = head(upgeneall,3000)
## 选取一定数目下调基因
dngene = tail(dngeneall,3000)
compareCluster 的输入文件可以是个list，里面放上不同组的基因
我们这里就输入上调，下调，以及全部基因
gcSample = list(up = upgene,
down = dngene,
all = c(upgene,dngene))
加载基因集,这里用hallmarks作为测试，可以在board下载得到
h_df <- read.gmt('resource/h.all.v2022.1.Hs.symbols.gmt')
正式运行,
dd <- compareCluster(gcSample, fun='enricher',TERM2GENE = h_df)
此处就是三个参数，
第一，gcSample就是我们给的多组数据，也可以是单细胞每个亚群的marker基因
第二，富集分析函数enriche,Y 叔的神奇之处在于，他创造了这个通用的富集分数
你不用再究竟KEGG分析，GO分析，只要是个基因集，他都能做，真正的解放思想。
第三，是enriche函数的第二个参数（可以逗号隔开更多的参数哦），就是这里的基因集。
以上的代码，可以实现任意基因集的ORA分析，可谓是真正的万能代码。
接下来就可以用极其简单熟悉的代码来实现作图啦
dotplot(dd,label_format = 60)
这是一个ESR1敲减的数据，所以在这个图上我们能看到两个雌激素相关的通路被富集了。
而且，是在下调的基因中富集到的，他十分符合生物学背景。
如果我们只用上调的基因是得到不到结果的，用全部的基因，可以富集但是不知道是激活和抑制。
因此，我觉得以后，这串代码应该成为差异基因富集分析的标配。
如果我们换一个数据集，比如KEGG的，结果就会有点不一样，
h_df <- read.gmt('resource/c2.cp.kegg.v2022.1.Hs.symbols.gmt')
dd <- compareCluster(gcSample, fun='enricher',TERM2GENE = h_df)
dotplot(dd,label_format = 60)
图出来了，但是我们发现少了下调的富集结果
这个结果本身没有问题，但是你放在文章中，需要额外解释，
我们也是做了下调基因的富集分析的，就是没有富集到。
这个还能理解，如果6个大群的单细胞数据进行这样的分群富集分析，缺少2-3个后，解释起来就不直观。
那么，现在的需求是，能不能把缺失的组给我显示出来，即使没有富集到任何结果？
答案是肯定的啊。
我们创建一个数据，把其中一个分组信息的因子水平改成3个
df <- data.frame(type=c('A', 'A', 'A', 'B', 'B'), group=rep('group1', 5))
df$type <- factor(df$type, levels=c('A','B', 'C'))
就是这样的数据，type总共就呈现两类数据，因子的水平是三个，ABC
直接画图
library(ggplot2)
ggplot(df, aes(x=type)) + geom_bar()
他只会显示两组
重点来了，X轴的分类，是由因子的水平决定的，只是在画图过程中，缺少内容的因子水平被去掉了
我们可以通过这句代码，把他保留下来。
scale_x_discrete(drop=FALSE)
重新画图看看
ggplot(df, aes(x=type)) + geom_bar() + scale_x_discrete(drop=FALSE)
看看C组，虽然没有东西，但是组别还是显示出来了，这就是我们想要的效果
知道了实现方法，接下来就是需要改造compareCluster的富集结果
图就是数据画出来的，所以我们要看看数据在哪里
Y叔这一系列的富集数据，都可以通过as.data.frame来获取
h_df <- read.gmt('resource/c2.cp.kegg.v2022.1.Hs.symbols.gmt')
dd <- compareCluster(gcSample, fun='enricher',TERM2GENE = h_df)
data <- as.data.frame(dd)
我们在数据中看到了X轴的数据，就是Cluster这一类，确实只有up和all这两类
那么我们只要改一下这一列的因子水平就可以了。
但是实际上最终不能实现，因为这还不是最终的画图数据。
因为Y叔在画图前，用了fortify这个函数对数据进行了预先处理，添加了额外的信息。
比如在这里，cluster实际上并不单纯是up和all
他还添加了有效输入基因数目(782和1429)，提供了更加丰富的信息。
所以现在不应该修改这个数据，而是最终画图的数据。
我曾经说过，Y叔对我最大的影响就是下面这句话
图就是数据，数据就是图。
因此我们现在图画出来，再把它变成数据，就可以啦
p <- dotplot(dd,label_format = 60)
此时p是一个数据，是个拥有10个元素的list
而画图数据data就藏在里面
这个cluster列才是我们真正需要的
提取这个数据，发现他是一个factor，因子水平只有2个
这就好办了，我们给他添加进去就行啦。
p$data$Cluster <- factor(p$data$Cluster,levels =c('up\n(782)','down\n(0)','all\n(1429)') )
此刻再来作图
p + scale_x_discrete(drop=FALSE)
完美呈现，我们可以感受到富集出来的通路都是上调基因的结果，可以认为是激活的，下调的基因当前没有贡献。
这个如果展示的是单细胞群marker的富集结果，可以表现出通路的特异性。
但是这个图是由瑕疵的，因为我在down那里给的数值是0，不应该是0，而应该是合法的输入基因数目。
这些数据，在Y叔的结果中都是提供的。
我们只需要让他不过滤，我们自己能看一眼就行啦
dd1 <- compareCluster(gcSample, fun='enricher',TERM2GENE = h_df,pvalueCutoff = 1,qvalueCutoff = 1)
data1 <- as.data.frame(dd1)
通过结果我们可以看到，上调的是782，下调的是647
此外，上调加上下调应该等于总体。所以，即使不做这个操作
我们做个减法也能得到结果
1429-782=647
此刻我们重新修改一下图片
p <- dotplot(dd)
p$data$Cluster <- factor(p$data$Cluster,levels =c('up\n(782)','down\n(647)','all\n(1429)') )
p + scale_x_discrete(drop=FALSE)
dd1 <- compareCluster(gcSample, fun='enricher',TERM2GENE = h_df,pvalueCutoff = 1,qvalueCutoff = 1)
data1 <- as.data.frame(dd1)
data1 <- data1 %>%
select(Cluster,GeneRatio) %>%
mutate(GeneRatio = gsub('\\d+/','',GeneRatio)) %>%
distinct(Cluster,.keep_all = T)
看看这个结果，再多的组都不怕了
如果你还不满意，不想在自定义因子水平的时候手工写，我还可以把代码再调整一下
dd1 <- compareCluster(gcSample, fun='enricher',TERM2GENE = h_df,pvalueCutoff = 1,qvalueCutoff = 1)
data1 <- as.data.frame(dd1)
data2 <- data1 %>%
select(Cluster,GeneRatio) %>%
mutate(GeneRatio = gsub('\\d+/','',GeneRatio)) %>%
distinct(Cluster,.keep_all = T) %>%
mutate(level = paste0(Cluster,'\n','(',GeneRatio,')'))
dput(data2$level)
他会自动生成这样的字符串,堪称懒癌福利
c('up\n(782)', 'down\n(647)', 'all\n(1429)')
我写这么多，实际上只要Y叔愿意，他可以轻松的加上一些代码，在神不知鬼不觉的情况下完成这些操作，
但是，当你看到这里的时候，你难道不觉得，自己能写代码，能自由探索更有趣么。
而这才是我真正的目的啊，激发你的兴趣，让你觉得，我也可以。
在这个例子中，我们的ESR1敲减数据。
使用hallmarks 基因集能够富集过雌激素通路，而且我们看到是抑制的。
但是使用广为人知的KEGG却没有能够得到有效信息，你认为是什么原因呢？
是不是应该不要把分析局限在特定的基因集，而是应该跳出去接受更多的信息呢。
而这个就是enricher这个函数的使命啊，他跟GSEA的用法几乎一致。
富集分析背后的支撑就是基因集，
而基因集是前人经验的积累和总结，
我们在探索数据时，应该光谱的把能做的基因集都尝试来一遍，然后再尝试整合这些前人的知识。
做到在没有开展实验前，胸有成竹。
而这个整合能力，又是功夫在诗外了啊。
下一次，我再来展示一下compareCluster 函数更加神奇的使用。
联系客服
