上下调基因各自独立进行GO数据库的3分类富集(求美图代码) / 开普饭

火山图大家应该是也基本上都没有问题，下面的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

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

上下调基因各自独立进行GO数据库的3分类富集(求美图代码)

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

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

学徒作业

相关推荐