接着上文得到差异基因之后,我们就会做一些常见的分析,在这里介绍的是基因富集分析。
首先,需要明白一下为什么要做富集分析?这可以帮助我们了解到这种差异主要参与了什么样的生物学功能中。那么何为富集呢?其实就是看某个通路中的基因在你的差异基因中的个数与背景相比是否有差异(超几何分布)。那具体怎么做呢?如果是人和小鼠的话,个人比较喜欢使用metascope数据库,这个还是很好使用的。如果要做其他的物种的话,我一般使用的是CluserProfile这个包(由余光创老师写的一个R包,听说这个老师还不错),接下来就是使用ClusterProfile进行富集的代码,在这里使用基因的SYMBOL进行富集:
ego <- enrichGO(diff$gene, OrgDb = "org.Mm.eg.db", ont="BP", keyType = "SYMBOL",)#只选择BP富集
但是很多人是不是也看到过文献中用GSEA做富集分析,其实GSEA本身也就是一个富集分析。不过它解决了普通富集通常会出现的一个问题,就是上调和下调的基因分别做富集分析,结果富集到了同一个通路。GSEA就解决了这个问题,因为普通的GO富集分析只考虑了差异基因的名字是否与某个term中的基因名字重叠;而GSEA还把相对表达量考虑了进去,就是用DESeq2的结果logFC;简单的理解就是把logFC作为了一个权重,这样的话就把上调下调的基因一起组做GSEA就会出现一个富集的结果。从原理上来讲GSEA是比GO富集更好的方式。接下来就是GSEA的富集代码:
1.由于GSEA考虑了相对的表达量,因此首先需要对logFC进行从大到小的排序:
diff<-diff[sort[diff$logFC,decreasing = TRUE],]
2.把排好序的基因表达矩阵进行GSEA富集分析:
diff<-data.frame(gene=diff$gene,logFC=diff$logFC)
geneList <- diff[,2]
names(geneList) <- as.character(diff[,1])
Go_gseresult <- gseGO(geneList, 'org.Hs.eg.db', keyType = "SYMBOL", ont="CC", nPerm = 1000,
minGSSize = 3, maxGSSize = 1000, pvalueCutoff=1) ##这里对CC进行GSEA富集
这样就完成了GSEA的富集分析了,拿到结果后就可以根据自己的情况进行画图了。 |