三阴性乳腺癌表达矩阵探索笔记之差异基因富集分析 / 开普饭

学徒和学员已经陆续出师，是时候把生信技能树的舞台交给后辈了！

下面是学徒写的《GEO数据挖掘课程》的配套笔记（第3篇）

前面作为了差异基因，但是是以探针为单位，这个时候探针需要转换为ENTREZID或者SYMBOL

第一种注方法：直接用别人已经做好的探针注释包来得到探针和基因的对应关系，在这种注释方法中，也可以直接得到探针和“ENTREZID”的对应关系。

select(hgu133plus2.db, keys = keys(hgu133plus2.db), columns = c("SYMBOL","ENTREZID"))

但是，并不是每一个数据集都有现成的注释包

第二种方法：一种通用的基因名转“ENTREZID的方法

library(ggplot2) library(clusterProfiler) library(org.Hs.eg.db) df <-

富集分析

统计学原理：超几何分布检验

总共有20000个基因，如果通路A总共有200个基因，那么该通路基因的比例为1%；

如果有200个基因上调，有2个基因参与到该通路，可以看做是正常现象，但是如果有4个或者更多基因，那就说明这些上调的基因富集到该通路

KEGG富集

library(ggplot2) library(clusterProfiler) library(org.Hs.eg.db) df <- bitr(unique(deg_1$symbol),fromType = "SYMBOL", toType = c("ENTREZID"), OrgDb = org.Hs.eg.db) #文件df是一个SYMBOL和ENTREZID的对应表 DEG <- merge(deg,df, by.y="SYMBOL",by.x="symbol") #很多的包都是根据ENTREZID来进行注释的


gene_up <- DEG[DEG$g=="up","ENTREZID"]

gene_down <- DEG[DEG$g=="down","ENTREZID"]

gene_diff <- c(gene_up, gene_down)

gene_all <- as.character(DEG[,"ENTREZID"])
kk.up <- enrichKEGG(gene = gene_up,

                    organism = "hsa",

                    universe = gene_all,

                    pvalueCutoff = 0.9,

                    qvalueCutoff = 0.9)

kk.down <- enrichKEGG(gene = gene_down, organism = "hsa", universe = gene_all, pvalueCutoff = 0.9, qvalueCutoff = 0.9) kk.diff <- enrichKEGG(gene = gene_down, organism = "hsa", universe = gene_all, pvalueCutoff = 0.9) 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 g_kegg = kegg_plot(up_kegg,down_kegg) print(g_kegg)

创建一个绘图函数keggplot

kegg_plot <- function(up_kegg,down_kegg){ dat=rbind(up_kegg,down_kegg) colnames(dat) dat$pvalue = -log10(dat$pvalue) dat$pvalue=dat$pvalue*dat$group dat=dat[order(dat$pvalue,decreasing = F),] g_kegg<- ggplot(dat, aes(x=reorder(Description,order(pvalue, decreasing = F)), y=pvalue, fill=group)) + geom_bar(stat="identity") + scale_fill_gradient(low="blue",high="red",guide = FALSE) + scale_x_discrete(name ="Pathway names") + scale_y_continuous(name ="log10P-value") + coord_flip() + theme_bw()+theme(plot.title = element_text(hjust = 0.5))+ ggtitle("Pathway Enrichment") }

GO富集

做GO数据集超几何分布检验分析，重点在于结果的可视化以及生物学意义的理解

g_list=list(gene_up=gene_up, gene_down=gene_down, gene_diff=gene_diff) # 因为go数据库非常多基因集，所以运行速度会很慢。 if(F){ go_enrich_results <- lapply( g_list , function(gene) { lapply( c('BP','MF','CC') , function(ont) { cat(paste('Now process ',ont )) ego <- enrichGO(gene = gene, universe = gene_all, OrgDb = org.Hs.eg.db, ont = ont , pAdjustMethod = "BH", pvalueCutoff = 0.99, qvalueCutoff = 0.99, #p值个q值根据自己的需要更改 readable = TRUE) print( head(ego) ) return(ego) }) }) save(go_enrich_results,file = 'go_enrich_results.Rdata') } # 只有第一次需要运行，就保存结果哈，下次需要探索结果，就载入即可。 load(file = 'go_enrich_results.Rdata') #分别对每个GO条目进行绘图 n1= c('gene_up','gene_down','gene_diff') n2= c('BP','MF','CC') for (i in 1:3){ for (j in 1:3){ fn=paste0('dotplot_',n1[i],'_',n2[j],'.png') cat(paste0(fn,'\n')) png(fn,res=150,width = 1080) print( dotplot(go_enrich_results[[i]][[j]] )) dev.off() } } }