上下调基因各自独立进行GO数据库的3分类富集

Posted on by

火山图大家应该是也基本上都没有问题，下面的MA图其实跟火山图非常的类似，两者都是log2FC信息，不同的是火山图展现P值，而MA图展现的是表达量情况！

  • 火山图是为了说明log2FC比较大的一般来说具有统计学显著性
  • 而MA图是为了说明log2FC无论大小，都不应该与表达量有相关性。

我们通常呢，挑选差异基因，会选择那些log2FC比较大而且具有统计学显著性的上下调基因，不过加上MA图，就可以进一步挑选那些表达量也比较高的，因为这样的基因呢，容易去实验验证。而且呢，通常情况下常识会告诉我们高表达量基因更容易发挥作用。

MA图加上GO数据库富集

图例：

  • (C) Differential gene expression between primary polyp calicoblasts and adult colony calicoblasts. Significant genes are shown in purple (adjusted p value <1e5, Fisher exact test).
  • (D) Gene Ontology analysis for differentially expressed genes between polyp and adult calicoblasts.

当然了，也可以对kegg进行同样的分析，上下调基因独立进行富集注释：

文献超级漂亮的kegg图

来源文献：Liu et al. J Hematol Oncol (2021) 14:6 https://doi.org/10.1186/s13045-020-01021-x ，不过要画出上面的图还是有一点难度的， 如果你一直看的是我的代码，你会发现出图简直是不忍直视：

我的kegg图

我的代码是：

## KEGG pathway analysis
### 做KEGG数据集超几何分布检验分析，重点在结果的可视化及生物学意义的理解。
if(T){
 ### over-representation test
 kk.up <- enrichKEGG(gene = gene_up,
 organism = 'hsa',
 universe = gene_all,
 pvalueCutoff = 0.9,#根据实际情况调整
 qvalueCutoff =0.9)#根据实际情况调整
 head(kk.up)[,1:6]
 dotplot(kk.up );ggsave('kk.up.dotplot.png')
 kk.down <- enrichKEGG(gene = gene_down,
 organism = 'hsa',
 universe = gene_all,
 pvalueCutoff = 0.9,#根据实际情况调整
 qvalueCutoff =0.9)#根据实际情况调整
 head(kk.down)[,1:6]
 dotplot(kk.down );ggsave('kk.down.dotplot.png')
 kk.diff <- enrichKEGG(gene = gene_diff,
 organism = 'hsa',
 pvalueCutoff = 0.05)#根据实际情况调整
 head(kk.diff)[,1:6]
 dotplot(kk.diff );ggsave('kk.diff.dotplot.png')

 kegg_diff_dt <- as.data.frame(kk.diff)
 kegg_down_dt <- as.data.frame(kk.down)
 kegg_up_dt <- as.data.frame(kk.up)
 down_kegg<-kegg_down_dt[kegg_down_dt$pvalue<0.05,];down_kegg$group=-1
 up_kegg<-kegg_up_dt[kegg_up_dt$pvalue<0.05,];up_kegg$group=1
 ## 原来的画图函数
 source('functions.R') 
 g_kegg=kegg_plot(up_kegg,down_kegg)
 print(g_kegg)

 ggsave(g_kegg,filename = 'kegg_up_down.png')

 ### GSEA 
 kk_gse <- gseKEGG(geneList = geneList,
 organism = 'hsa',
 nPerm = 1000,#根据实际情况调整
 minGSSize = 10,#根据实际情况调整
 pvalueCutoff = 0.9,##根据实际情况调整
 verbose = FALSE)
 head(kk_gse)[,1:6]
 gseaplot(kk_gse, geneSetID = rownames(kk_gse[1,]))

 down_kegg<-kk_gse[kk_gse$pvalue<0.05 & kk_gse$enrichmentScore < 0,];down_kegg$group=-1
 up_kegg<-kk_gse[kk_gse$pvalue<0.05 & kk_gse$enrichmentScore > 0,];up_kegg$group=1

 g_kegg=kegg_plot(up_kegg,down_kegg)
 print(g_kegg)
 ggsave(g_kegg,filename = 'kegg_up_down_gsea.png')


}

然后很久以前的一个学徒在我的绘图代码上面进行了修饰：


### Update: Xiaojie Liang
### Date: 2021-07-12 23:38:07
### Email: liangxiaojiecom@163.com
### ---------------
library(ggthemes)
library(ggplot2)
go.kegg_plot <- function(up.data,down.data){
 # up.data <- up.data
 # down.data <- down.data
 dat=rbind(up.data,down.data)
 colnames(dat)
 dat$p.adjust = -log10(dat$p.adjust)
 dat$p.adjust=dat$p.adjust*dat$group 

 dat=dat[order(dat$p.adjust,decreasing = F),]

 gk_plot <- ggplot(dat,aes(reorder(Description, p.adjust), y=p.adjust)) +
 geom_bar(aes(fill=factor((p.adjust>0)+1)),stat="identity", width=0.7, position=position_dodge(0.7)) +
 coord_flip() +
 scale_fill_manual(values=c("#0072B2", "#B20072"), guide=FALSE) +
 labs(x="", y="" ) +
 theme_pander() +
 theme(panel.grid.major = element_blank(),
 panel.grid.minor = element_blank(),
 #axis.ticks.x = element_blank(),
 axis.line.x = element_line(size = 0.3, colour = "black"),#x轴连线
 axis.ticks.length.x = unit(-0.20, "cm"),#修改x轴刻度的高度，负号表示向上
 axis.text.x = element_text(margin = margin(t = 0.3, unit = "cm")),##线与数字不要重叠
 axis.ticks.x = element_line(colour = "black",size = 0.3) ,#修改x轴刻度的线 
 axis.ticks.y = element_blank(),
 axis.text.y = element_text(hjust=0),
 panel.background = element_rect(fill=NULL, colour = 'white')
 )
}

效果如下：

学徒的kegg图

啊，扯远了，我们这个教程其实是想征集大家的代码，关于上下调基因各自独立进行GO数据库的3分类富集后的可视化：

首先进行上下调基因的ID转换

前面的表达量矩阵差异分析代码我这里就不再赘述啦，我们分析已经拿到了 gene_up 和 gene_down 这两个变量（简单的向量而已）代表的上下调差异基因哦！

library(org.Hs.eg.db)
#把SYMBOL转换ENTREZID，这里会损失一部分无法匹配到的
gene_up=as.character(na.omit(AnnotationDbi::select(org.Hs.eg.db,
 keys = gene_up,
 columns = 'ENTREZID',
 keytype = 'SYMBOL')[,2]))
gene_down=as.character(na.omit(AnnotationDbi::select(org.Hs.eg.db,
 keys = gene_down,
 columns = 'ENTREZID',
 keytype = 'SYMBOL')[,2]))

通常使用clusterProfiler包做GO和KEGG数据库富集需要的是 ENTREZID的ID系统，所以转换一下。

然后上下调基因独立进行GO数据库富集

代码仍然是很简单：

library(clusterProfiler) 
#对正相关基因进行富集分析
go <- enrichGO(gene_up, OrgDb = "org.Hs.eg.db", ont="all") 
library(ggplot2)
library(stringr)
barplot(go, split="ONTOLOGY")+ facet_grid(ONTOLOGY~., scale="free")+ 
 ggsave('gene_up_GO_all_barplot.png') 
go=DOSE::setReadable(go, OrgDb='org.Hs.eg.db',keyType='ENTREZID')
write.csv( go@result,file = 'gene_up_GO_all_barplot.csv')

#对负相关基因进行富集分析
go <- enrichGO(gene_down, OrgDb = "org.Hs.eg.db", ont="all") 
barplot(go, split="ONTOLOGY")+ facet_grid(ONTOLOGY~., scale="free")+ 
 ggsave('gene_down_GO_all_barplot.png') 
go=DOSE::setReadable(go, OrgDb='org.Hs.eg.db',keyType='ENTREZID')
write.csv( go@result,file = 'gene_down_GO_all_barplot.csv')

这样有了最原始的图，如何绘制成为文章开头的美图呢？

学徒作业

看我六年前的表达芯片的公共数据库挖掘系列推文 ：

完成数据集：https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE97251 的表达量矩阵差异分析 ：

GSM2560200 plasmaRNA_health H13
GSM2560201 plasmaRNA_health H5
GSM2560202 plasmaRNA_health H8
GSM2560203 plasmaRNA_OSCC X474
GSM2560204 plasmaRNA_OSCC X574
GSM2560205 plasmaRNA_OSCC X602

这个数据集是很简单的两个分组，芯片平台是 Agilent-079487 Arraystar Human LncRNA microarray V4

如果确实感兴趣也可以去和原文对比：Screening and validation of plasma long non-coding RNAs as biomarkers for the early diagnosis and staging of oral squamous cell carcinoma. Oncol Lett 2021 Feb;21(2):172. PMID: 33552289

我需要大家完成文章开头的美图，并且给出来代码！

​

Comments are closed.