scripts/topGO_2.4.R at master · PBR/scripts · GitHub

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
getwd()
setwd('/media/bulk_01/users/goryu003/samples')
ls()


#load required packages

library(BiocGenerics, lib.loc="/media/bulk_01/users/goryu003/R/x86_64-pc-linux-gnu-library/3.1/")

library(AnnotationDbi, lib.loc="/media/bulk_01/users/goryu003/R/x86_64-pc-linux-gnu-library/3.1/")

library(S4Vectors, lib.loc="/media/bulk_01/users/goryu003/R/x86_64-pc-linux-gnu-library/3.1/")

library(graph, lib.loc="/media/bulk_01/users/goryu003/R/x86_64-pc-linux-gnu-library/3.1/")

library(GO.db, lib.loc="/media/bulk_01/users/goryu003/R/x86_64-pc-linux-gnu-library/3.1/")

library(Rgraphviz, lib.loc="/media/bulk_01/users/goryu003/R/x86_64-pc-linux-gnu-library/3.1/")

library(topGO, lib.loc="/media/bulk_01/users/goryu003/R/x86_64-pc-linux-gnu-library/3.1/")


#Create list of files to use in loop

list_files=dir (pattern='.*_F-Change_signif.txt')

#Create topGO objects

for (i in 1:length( list_files))

    {
        geneID2GO <- readMappings(file = "GO_2-4.txt",sep="")

        L=lapply(geneID2GO,function(x) {gsub('\"',"",x)})

        geneID2GO=L

        names(geneID2GO)<- gsub('\"',"",names(geneID2GO))

        geneUniverse <- names(geneID2GO)

        tableOfInterest=read.table(list_files[[i]],as.is=T,header=T)

        genesOfInterest=tableOfInterest$gene

        geneList <- factor(as.integer(geneUniverse %in% genesOfInterest))

        names(geneList) <- geneUniverse

        myGOdata_BP <- new("topGOdata", description="My project", ontology="BP", allGenes=geneList,  annot = annFUN.gene2GO, gene2GO = geneID2GO, nodeSize=5)

        print(myGOdata_BP)

        myGOdata_CC <- new("topGOdata", description="My project", ontology="CC", allGenes=geneList,  annot = annFUN.gene2GO, gene2GO = geneID2GO, nodeSize=5)

        print(myGOdata_CC)

        myGOdata_MF <- new("topGOdata", description="My project", ontology="MF", allGenes=geneList,  annot = annFUN.gene2GO, gene2GO = geneID2GO, nodeSize=5)

        print(myGOdata_MF)

        #Perform the analysis

        #Biological process

        resultFisher <- runTest(myGOdata_BP, algorithm="weight01", statistic="fisher")

        allRes <- GenTable(myGOdata_BP, weight01Fisher = resultFisher, topNodes = length(resultFisher@score))

        print(allRes)

        allRes_2=allRes[as.numeric(allRes$weight01Fisher)<0.05,]

        print(allRes_2)

        write.table(allRes_2,file=paste(sep="",strsplit(list_files[[i]], "\\.")[[1]][1],"_BP_enrichment.txt"))

        showSigOfNodes(myGOdata_BP, score(resultFisher), firstSigNodes = 5, useInfo ='all')

        printGraph(myGOdata_BP, resultFisher, firstSigNodes = nrow(allRes_2), fn.prefix = paste(sep="",strsplit(list_files[[i]], "\\.")[[1]][1],"_BP_enrichment", useInfo = "all", pdfSW = TRUE))

        myterms=allRes_2[,1]

        mygenes=genesInTerm(myGOdata_BP, myterms)

        mygenes_2=unlist(mygenes)

        df=data.frame(gene_name=mygenes_2,GO_term=names(mygenes_2))

        write.table(df, paste(sep="",strsplit(list_files[[i]], "\\.")[[1]][1],"genes_corresponding_overrepresent_BP_GO.txt"))

        rm(resultFisher)

        #Cellular component

        resultFisher <- runTest(myGOdata_CC, algorithm="weight01", statistic="fisher")

        allRes <- GenTable(myGOdata_CC, weight01Fisher = resultFisher, topNodes = length(resultFisher@score))

        print(allRes)

        allRes_2=allRes[as.numeric(allRes$weight01Fisher)<0.05,]

        print(allRes_2)

        write.table(allRes_2,file=paste(sep="",strsplit(list_files[[i]], "\\.")[[1]][1],"_CC_enrichment.txt"))

        showSigOfNodes(myGOdata_CC, score(resultFisher), firstSigNodes = 5, useInfo ='all')

        printGraph(myGOdata_CC, resultFisher, firstSigNodes = nrow(allRes_2), fn.prefix = paste(sep="",strsplit(list_files[[i]], "\\.")[[1]][1],"_CC_enrichment", useInfo = "all", pdfSW = TRUE))

        myterms=allRes_2[,1]

        mygenes=genesInTerm(myGOdata_CC, myterms)

        mygenes_2=unlist(mygenes)

        df=data.frame(gene_name=mygenes_2,GO_term=names(mygenes_2))

        write.table(df, paste(sep="",strsplit(list_files[[i]], "\\.")[[1]][1],"genes_corresponding_overrepresent_CC_GO.txt"))

        rm(resultFisher)


        #Molecular function

        resultFisher <- runTest(myGOdata_MF, algorithm="weight01", statistic="fisher")

        allRes <- GenTable(myGOdata_MF, weight01Fisher = resultFisher, topNodes = length(resultFisher@score))

        print(allRes)

        allRes_2=allRes[as.numeric(allRes$weight01Fisher)<0.05,]

        print(allRes_2)

        write.table(allRes_2,file=paste(sep="",strsplit(list_files[[i]], "\\.")[[1]][1],"_MF_enrichment.txt"))

        showSigOfNodes(myGOdata_MF, score(resultFisher), firstSigNodes = 5, useInfo ='all')

        printGraph(myGOdata_MF, resultFisher, firstSigNodes = nrow(allRes_2), fn.prefix = paste(sep="",strsplit(list_files[[i]], "\\.")[[1]][1],"_MF_enrichment", useInfo = "all", pdfSW = TRUE))

        myterms=allRes_2[,1]

        mygenes=genesInTerm(myGOdata_MF, myterms)

        mygenes_2=unlist(mygenes)

        df=data.frame(gene_name=mygenes_2,GO_term=names(mygenes_2))

        write.table(df, paste(sep="",strsplit(list_files[[i]], "\\.")[[1]][1],"genes_corresponding_overrepresent_MF_GO.txt"))

        rm(resultFisher)


    }