搜索
查看: 2431|回复: 0

[mRNA-seq] 转录组入门分析8_差异表达基因富集分析

[复制链接]

9

主题

15

帖子

132

积分

注册会员

Rank: 2

积分
132
发表于 2017-9-3 22:46:24 | 显示全部楼层 |阅读模式
本帖最后由 laofuzi 于 2017-9-3 22:47 编辑

                                
[AppleScript] 纯文本查看 复制代码
########################
#select diff_gene_id
diff_gene_id<-row.names(subset(res,padj < 0.05 & (log2FoldChange > 1| log2FoldChange < -1)))

                               
#install and load librarys
source("https://bioconductor.org/biocLite.R")
biocLite("clusterProfiler",dependencies=T)
library(clusterProfiler)
biocLite("org.Mm.eg.db")
library(org.Mm.eg.db)

##GOidmapping + GO enrich
ego_BP<- enrichGO(gene = diff_gene_id,
                    OrgDb=org.Mm.eg.db,
                    keytype= "ENSEMBL",
                    ont= "BP",
                    pAdjustMethod= "BH",
                    pvalueCutoff= 0.01,
                    qvalueCutoff= 0.05
                    )
ego_CC<- enrichGO(gene = diff_gene_id, OrgDb = org.Mm.eg.db, keytype ="ENSEMBL", ont = "CC", pAdjustMethod = "BH",pvalueCutoff = 0.05, qvalueCutoff = 0.05)
ego_MF<- enrichGO(gene = diff_gene_id, OrgDb = org.Mm.eg.db, keytype ="ENSEMBL", ont = "MF", pAdjustMethod = "BH",pvalueCutoff = 0.05, qvalueCutoff = 0.05)

#viewstructure of ego_*
ego_BP
ego_CC
ego_MF
##output the result of GO enrich
write.csv(as.data.frame(ego_BP),file= "/home/R/AKAP95/matrix/GO_enriched_BP.csv",row.names= F)
write.csv(as.data.frame(ego_CC),file= "/home/R/AKAP95/matrix/GO_enriched_CC.csv",row.names= F)
write.csv(as.data.frame(ego_MF),file= "/home/R/AKAP95/matrix/GO_enriched_MF.csv",row.names= F)

################Draw figures
##barplot
pdf(file="/home/R/AKAP95/matrix/barplot_BP.pdf",onefile = FALSE, width=12, height=8)
barplot(ego_BP,showCategory=10, font.size =5)
dev.off()

pdf(file="/home/R/AKAP95/matrix/barplot_CC.pdf",onefile = FALSE, width=12, height=8)
barplot(ego_CC,showCategory=10, font.size =5)
dev.off()

pdf(file="/home/R/AKAP95/matrix/barplot_MF.pdf",onefile = FALSE, width=12, height=8)
barplot(ego_MF,showCategory=10, font.size =5)
dev.off()


##dotplot
pdf(file="/home/R/AKAP95/matrix/dotplot_BP.pdf",onefile = FALSE, width=12, height=8)
dotplot(ego_BP,showCategory=10, font.size =5)
dev.off()

pdf(file="/home/R/AKAP95/matrix/dotplot_CC.pdf",onefile = FALSE, width=12, height=8)
dotplot(ego_CC,showCategory=10, font.size =5)
dev.off()

pdf(file="/home/R/AKAP95/matrix/dotplot_MF.pdf",onefile = FALSE, width=12, height=8)
dotplot(ego_MF,showCategory=10, font.size =5)
dev.off()


##enrichMap
pdf(file="/home/R/AKAP95/matrix/enrichMap_BP.pdf",onefile = FALSE, width=12, height=8)
enrichMap(ego_BP,vertex.label.cex=1.2, layout=igraph::layout.kamada.kawai)
dev.off()

pdf(file="/home/R/AKAP95/matrix/enrichMap_CC.pdf",onefile = FALSE, width=12, height=8)
enrichMap(ego_CC,vertex.label.cex=1.2, layout=igraph::layout.kamada.kawai)
dev.off()

pdf(file="/home/R/AKAP95/matrix/enrichMap_MF.pdf",onefile = FALSE, width=12, height=8)
enrichMap(ego_MF,vertex.label.cex=1.2, layout=igraph::layout.kamada.kawai)
dev.off()


##netplot
pdf(file="/home/R/AKAP95/matrix/plotGOgraph_BP.pdf",onefile = FALSE, width=12, height=8)
plotGOgraph(ego_BP)
dev.off()

pdf(file="/home/R/AKAP95/matrix/plotGOgraph_CC.pdf",onefile = FALSE, width=12, height=8)
plotGOgraph(ego_CC)
dev.off()

pdf(file="/home/R/AKAP95/matrix/plotGOgraph_MF.pdf",onefile = FALSE, width=12, height=8)
plotGOgraph(ego_MF)
dev.off()

################KEGG
##GeneID transfer
ids<- bitr(diff_gene_id, fromType = "ENSEMBL", toType ="ENTREZID", OrgDb = "org.Mm.eg.db")
##KEGGenrichment
kk<- enrichKEGG(gene = ids[,2], organism = "mmu", keyType= "kegg", pvalueCutoff = 0.05, pAdjustMethod = "BH",qvalueCutoff = 0.1)
#viewresults
enrich_kegg<- as.data.frame(kk)
head(summary(enrich_kegg))
write.csv(as.data.frame(enrich_kegg),file= "/home/R/AKAP95/matrix/KEGG-enrich.csv", row.names=F)
#plot
pdf(file="/home/R/AKAP95/matrix/KEGG-enrich.pdf",onefile = FALSE, width=12, height=8)
dotplot(enrich_kegg,font.size=5)
dev.off()
后记:
跌跌撞撞的基本上把RNA-seq分析做完了,学习了很多。但目前还有很多尚未完全理解。然而,这并不是结束。有的时候,我们知道的越多,我们了解的就更少。(Sometime the more we know, the less we understand.)对于我来说,这仅仅是开始(This is just a end of beginning). Jimmy大神已经发布了ChiP-seq学习的号召。我准备开始了,你呢?




上一篇:将肿瘤进化作为下一个靶标 tumor evolution as a therapeutic target
下一篇:研究完了基因,还有其编码蛋白的domain
回复

使用道具 举报

您需要登录后才可以回帖 登录 | 立即注册

本版积分规则

QQ|手机版|小黑屋|生信技能树 ( 粤ICP备15016384号  

GMT+8, 2019-7-24 09:19 , Processed in 0.029066 second(s), 28 queries .

Powered by Discuz! X3.2

© 2001-2013 Comsenz Inc.