本文将介绍RNA-seq基因富集分析的实战代码。
原文地址:https://bioinformatics-core-shared-training.github.io/RNAseq-R/rna-seq-gene-set-testing.nb.html
书接上文(用R也可以完成的RNA-Seq下游分析-4)。
在基因差异分析后,我们会得到一个包含了差异表达基因信息的表格。但有时过于冗长的基因列表或是繁杂的信息会阻碍我们从中提取有效的信息。因此,科学家们时常会采用基因富集分析的办法,去掌握差异表达基因真正影响到的信号通路或是代谢途径等。
例如,Gene Ontology(GO)或是 KEGG就是较为常用的基因富集分析,不仅可以帮助我们从冗长的数据中提取出有生物学意义的信息,还能有助于我们从系统的水平观察细胞/机体的变化。
接下来,我们就以GOseq
这个包向大家展示在R中处理的GO富集分析流程。当然我们只是简要地介绍一下,如果想要深入了解的朋友,可以查看Bioconductor上的帮助文档:GOseq
本次流程所需要的包和数据集
library(edgeR)
library(goseq)
# 此前差异分析的结果
load("DE.Rdata")
GOseq analysis pipeline
results <- as.data.frame(topTags(lrt.BvsL, n = Inf))
# 筛选显著的差异表达基因进行分析
genes <- results$FDR < 0.01
names(genes) <- rownames(results)
# Fit the Probability Weighting Function (PWF) ,以生成零分布
pwf <- nullp(genes, "mm10","knownGene")
# 富集分析
go.results <- goseq(pwf, "mm10","knownGene")
head(go.results,10)
category over_represented_pvalue under_represented_pvalue numDEInCat numInCat
16121 GO:0071944 4.062458e-63 1 933 3543
2723 GO:0005886 2.686421e-62 1 908 3420
2540 GO:0005576 4.492776e-57 1 469 1457
11253 GO:0044421 1.026496e-50 1 391 1175
11283 GO:0044459 2.158946e-49 1 563 1893
2562 GO:0005615 2.018827e-47 1 336 975
7483 GO:0031224 2.789374e-41 1 833 3414
7485 GO:0031226 3.115084e-40 1 323 923
2724 GO:0005887 7.600430e-40 1 309 867
4426 GO:0009653 1.471331e-39 1 592 2186
term ontology
16121 cell periphery CC
2723 plasma membrane CC
2540 extracellular region CC
11253 extracellular region part CC
11283 plasma membrane part CC
2562 extracellular space CC
7483 intrinsic component of membrane CC
7485 intrinsic component of plasma membrane CC
2724 integral component of plasma membrane CC
4426 anatomical structure morphogenesis BP
GO的分析也就到此为止了,至于后续的可视化的话就不先在此展开了。实际上还有许多R包也能在R中进行基因富集分析,例如Y叔的clusterProfiler
,fgsea
,camera
等等。
关于RNA-seq分析流程到这就先告一段落了,文章中有任何错误都请各位小伙伴指出,共同学习!
完。