学徒和学员已经陆续出师，是时候把生信技能树的舞台交给后辈了！
第一种注方法：直接用别人已经做好的探针注释包来得到探针和基因的对应关系，在这种注释方法中，也可以直接得到探针和“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个或者更多基因，那就说明这些上调的基因富集到该通路
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数据集超几何分布检验分析，重点在于结果的可视化以及生物学意义的理解
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()
}
}
}
最后得到9个图，以上调基因为例子展示结果。
点阵图（dotplot)怎么看？
MF（molecular function）
BP（biological process）
CC(cellular components)
这些分析，基本上读一下我五年前在生信技能树的表达芯片的公共数据库挖掘系列推文 就明白了；
我把3年前的收费视频课程：3年前的GEO数据挖掘课程你可以听3小时或者3天甚至3个月，免费到B站:
然后马上就有了3千多学习量，而且有学员给出来了图文并茂版本万字笔记，让我非常感动！
扫描下面二维码马上就可以学习起来啦，笔记需要至少半个小时来阅读哦！
联系客服
