GO富集出来的结果可以非常多,如果没有一个层级显示,将会非常眼花缭乱。TopGO则通过DAG(Directed acyclic graph)解决了这个问题。
TopGO的详细信息,都可以在bioconductor中的TopGO页面找到,网站中附有详细说明。下面代码不涉及其他统计分析,只涉及出DAG图所需要的的最简单的代码。
安装:
if (!requireNamespace("BiocManager", quietly = TRUE))
install.packages("BiocManager")
BiocManager::install("topGO")
文件准备:
-
所有基因列表。这个可以在excel中做成两列数据,这里是RNA-seq数据,因此第一列为基因ID,第二列可以为0或者1,其中“0”表示非差异表达,“1”表示差异表达。如果是做芯片之类的,也可以是基因分值或者p-value值,具体视实验而定。接着将内容复制到txt文件中。将数据框先转变成矩阵,再转变成因子即可。
基因列表 - 差异基因列表。这个可以通过对“0”和“1”进行筛选,或者直接对所有基因列表进行排序,将差异基因,也就是“1”的基因排在前面,直接截取即可。
-
基因和GO对应的关系列表。这个可以是Gene-GO,或者GO-Gene。以GO-Gene为例,这个文件还是两列,第一列是GO terms,第二列是用逗号分开的Genes,如下图:
GO-Gene
这个文件不能直接读入,而是用readMappings
函数进行解释。TopGO里面有一个很有用的function,就是inverseList
,可以把Gene-GO和GO-Gene两种文件相互转换。
代码
setwd("D:/BioInfo/R/Package information/TopGO/Test")
library(topGO)
#01. read the GeneList file and convert it to 2 levels factor
GL = read.table("01.GL.txt")
GL1 = as.matrix(GL)
rownames(GL1) = GL1[,1]
GL2 = GL1[,2]
GL3 = as.factor(GL2)
#02. Directly generate DEGs,DEGs with value of "1" were at the top of the list.
DEG = GL3[1:200]
#03. parse the annotation file, "readMappings" function was used to parse both of the GO-Gene or Gene-GO files
GeneGO = readMappings("03.GeneGO.txt")
## inverse the direction
## Notably, use the readMappings function firstly, then inverseList function.
GOGene = inverseList(GeneGO)
#04. Generate GOdata, three parts should be run separately, including "BP","MF" and "CC"
GOdata = new("topGOdata",
ontology = "BP",
allGenes = GL3,
geneSelectionFun = DEG,
annot = annFUN.gene2GO,
nodeSize = 5,
gene2GO = GeneGO
)
#05. statistic
resultFis=runTest(GOdata,algorithm = "classic",statistic = "fisher")
#06. Visualization, firstSigNodes could be changed.
showSigOfNodes(GOdata,
score(resultFis),
firstSigNodes = 5,
useInfo = 'all')
#07. Print
printGraph(GOdata,
resultFis,
firstSigNodes = 5,
fn.prefix = "tGO",
useInfo = "all",
pdfSW = TRUE)
结果如下:
Top5的MF结果显示
直接在Rstudio显示不清晰,但输出成EPS或者PDF就清晰了。
另外一个问题就是默认输出全部信息的时候,GO terms太长的会显示不全,这个时候我们可以用AI打开文件,并进行相应修改,再导出tiff文件即可。
网友评论