GO enrichment with all significant GOs and GOSemSim script

MikheyevLab · Jun 25, 2019 · d788911 · d788911
1 parent 8f1a66e
commit d788911
Show file tree

Hide file tree

Showing 54 changed files with 61,625 additions and 148 deletions.
diff --git a/.DS_Store b/.DS_Store
diff --git a/GOSemSim.R b/GOSemSim.R
@@ -0,0 +1,35 @@
+library(GOSemSim)
+library(org.Dm.eg.db)
+#Reference for Fly
+d = godata('org.Dm.eg.db', ont = "BP", computeIC = FALSE)
+contrast = read.csv("contrasts_mite_life_stages.csv", header = F)
+similarities_upregulated = c()
+similarities_downregulated = c()
+
+for (id in 1:length(contrast$V1))
+{
+  first_contrast = toString(contrast$V1[id])
+  second_contrast = toString(contrast$V2[id])
+
+  #UP-REGULATED
+  filename = paste("./GO_enrichment_microbiome_output/contrast_", first_contrast, "_", second_contrast, "_upregulated", ".csv", sep = "")
+  GO1 = read.csv(filename, header = T)
+  GO1ID = as.character(GO1[,1])
+
+  filename = paste("./GO_enrichment_output/contrast_", first_contrast, "_", second_contrast, "_upregulated", ".csv", sep = "")
+  GO2 = read.csv(filename, header = T)
+  GO2ID = as.character(GO2[,1])
+
+  similarities_upregulated[id] = mgoSim(GO1ID, GO2ID, semData = d, measure = "Wang", combine = "BMA")
+
+  #DOWN-REGULATED
+  filename = paste("./GO_enrichment_microbiome_output/contrast_", first_contrast, "_", second_contrast, "_downregulated", ".csv", sep = "")
+  GO1 = read.csv(filename, header = T)
+  GO1ID = as.character(GO1[,1])
+
+  filename = paste("./GO_enrichment_output/contrast_", first_contrast, "_", second_contrast, "_downregulated", ".csv", sep = "")
+  GO2 = read.csv(filename, header = T)
+  GO2ID = as.character(GO2[,1])
+
+  similarities_downregulated[id] = mgoSim(GO1ID, GO2ID, semData = d, measure = "Wang", combine = "BMA")
+}
diff --git a/GO_enrichment_Microbiome.R b/GO_enrichment_Microbiome.R
@@ -30,7 +30,7 @@ samples <- tibble(id = dir(file.path("mondet/kallisto"))) %>% mutate(path = past
 contrast = read.csv("contrasts_mite_life_stages.csv", header = F)
 gene2GO = readMappings2("GOterm/vdesgoassoc.csv")
 t2g = read.csv("IDrna_gene_map.csv", header=T)
-for (id in 4:length(contrast$V1))
+for (id in 1:length(contrast$V1))
 {
   print(paste(id, "id"))
   first_contrast = toString(contrast$V1[id])
@@ -44,6 +44,7 @@ for (id in 4:length(contrast$V1))
   so = sleuth_wt(so, cond)
   results_table = sleuth_results(so, cond)
   #DE_genes = subset(results_table, qval <= 0.05)
+
   #UP-REGULATED
   up_regulated = subset(results_table, b > 0) #beta postive is up-regulated
 
@@ -62,24 +63,12 @@ for (id in 4:length(contrast$V1))
   resultKS <- runTest(sampleGOdata, algorithm = "classic", statistic = "ks")
   resultKS.elim <- runTest(sampleGOdata, algorithm = "elim", statistic = "ks")
 
-  allRes <- GenTable(sampleGOdata, classicFisher = resultFisher, classicKS = resultKS, elimKS = resultKS.elim, orderBy = "elimKS", ranksOf = "classicFisher", topNodes = 5)
+  allGO = usedGO(sampleGOdata)
+  allRes <- GenTable(sampleGOdata, classicFisher = resultFisher, classicKS = resultKS, elimKS = resultKS.elim, orderBy = "elimKS", ranksOf = "classicFisher", topNodes = length(allGO))
 
   filename = paste("./GO_enrichment_microbiome_output/contrast_", first_contrast, "_", second_contrast, "_upregulated", ".csv", sep = "")
   write.csv(allRes, file = filename, row.names = F)
 
-  pValue.classic <- score(resultKS)
-  pValue.elim <- score(resultKS.elim)[names(pValue.classic)]
-  gstat <- termStat(sampleGOdata, names(pValue.classic))
-  gSize <- gstat$Annotated / max(gstat$Annotated) * 4
-  gCol <- colMap(gstat$Significant)
-
-  plotfile = paste("./GO_enrichment_microbiome_output/contrast_", first_contrast, "_", second_contrast, "_upregulated", ".png", sep = "")
-  png(filename=plotfile)
-  plot(pValue.classic, pValue.elim, xlab = "p-value classic", ylab = "p-value elim",
-       pch = 19, cex = gSize, col = gCol, main=paste(first_contrast, second_contrast, "up regulated"))
-  dev.off()
-
-  showSigOfNodes(sampleGOdata, score(resultKS.elim), firstSigNodes = 5, useInfo = 'all')
 
   #DOWN-REGULATED
   down_regulated = subset(results_table, b < 0) #beta negative is down-regulated
@@ -99,23 +88,12 @@ for (id in 4:length(contrast$V1))
   resultKS <- runTest(sampleGOdata, algorithm = "classic", statistic = "ks")
   resultKS.elim <- runTest(sampleGOdata, algorithm = "elim", statistic = "ks")
 
-  allRes <- GenTable(sampleGOdata, classicFisher = resultFisher, classicKS = resultKS, elimKS = resultKS.elim, orderBy = "elimKS", ranksOf = "classicFisher", topNodes = 5)
+  allGO = usedGO(sampleGOdata)
+  allRes <- GenTable(sampleGOdata, classicFisher = resultFisher, classicKS = resultKS, elimKS = resultKS.elim, orderBy = "elimKS", ranksOf = "classicFisher", topNodes = length(allGO))
 
   filename = paste("./GO_enrichment_microbiome_output/contrast_", first_contrast, "_", second_contrast, "_downregulated", ".csv", sep = "")
   write.csv(allRes, file = filename, row.names = F)
 
-  pValue.classic <- score(resultKS)
-  pValue.elim <- score(resultKS.elim)[names(pValue.classic)]
-  gstat <- termStat(sampleGOdata, names(pValue.classic))
-  gSize <- gstat$Annotated / max(gstat$Annotated) * 4
-  gCol <- colMap(gstat$Significant)
-
-  plotfile = paste("./GO_enrichment_microbiome_output/contrast_", first_contrast, "_", second_contrast, "_downregulated", ".png", sep = "")
-  png(filename=plotfile)
-  plot(pValue.classic, pValue.elim, xlab = "p-value classic", ylab = "p-value elim",
-       pch = 19, cex = gSize, col = gCol, main=paste(first_contrast, second_contrast, "down regulated"))
-  dev.off()
 
-  showSigOfNodes(sampleGOdata, score(resultKS.elim), firstSigNodes = 5, useInfo = 'all')
 
 }
diff --git a/GO_enrichment_comparison_mite.xlsx b/GO_enrichment_comparison_mite.xlsx
diff --git a/GO_enrichment_loop.R b/GO_enrichment_loop.R
@@ -50,7 +50,7 @@ for (id in 1:length(contrast$V1))
   so = sleuth_wt(so, cond)
   results_table = sleuth_results(so, cond)
   DE_genes = subset(results_table, qval <= 0.05)
-
+#UP-REGULATED
   up_regulated = subset(results_table, b > 0) #beta postive is up-regulated
 
   #Preapre inputs for sampleOGData of class topGOdata
@@ -69,27 +69,14 @@ for (id in 1:length(contrast$V1))
   #Kolmogorov_smirnov - based on gene scores
   resultKS <- runTest(sampleGOdata, algorithm = "classic", statistic = "ks")
   resultKS.elim <- runTest(sampleGOdata, algorithm = "elim", statistic = "ks")
-
-  allRes <- GenTable(sampleGOdata, classicFisher = resultFisher, classicKS = resultKS, elimKS = resultKS.elim, orderBy = "elimKS", ranksOf = "classicFisher", topNodes = 5)
+  allGO = usedGO(sampleGOdata)
+  allRes <- GenTable(sampleGOdata, classicFisher = resultFisher, classicKS = resultKS, elimKS = resultKS.elim, orderBy = "elimKS", ranksOf = "classicFisher", topNodes = length(allGO))
 
   filename = paste("./GO_enrichment_output/contrast_", first_contrast, "_", second_contrast, "_upregulated", ".csv", sep = "")
   write.csv(allRes, file = filename, row.names = F)
 
-  pValue.classic <- score(resultKS)
-  pValue.elim <- score(resultKS.elim)[names(pValue.classic)]
-  gstat <- termStat(sampleGOdata, names(pValue.classic))
-  gSize <- gstat$Annotated / max(gstat$Annotated) * 4
-  gCol <- colMap(gstat$Significant)
-
-  plotfile = paste("./GO_enrichment_output/contrast_", first_contrast, "_", second_contrast, "_upregulated", ".png", sep = "")
-  png(filename=plotfile)
-  plot(pValue.classic, pValue.elim, xlab = "p-value classic", ylab = "p-value elim",
-       pch = 19, cex = gSize, col = gCol, main=paste(first_contrast, second_contrast, "up regulated"))
-  dev.off()
-
-  showSigOfNodes(sampleGOdata, score(resultKS.elim), firstSigNodes = 5, useInfo = 'all')
-
 
+  #DOWN-REGULATED
   down_regulated = subset(results_table, b < 0) #beta postive is up-regulated
 
     #Preapre inputs for sampleOGData of class topGOdata
@@ -108,25 +95,13 @@ for (id in 1:length(contrast$V1))
     #Kolmogorov_smirnov - based on gene scores
     resultKS <- runTest(sampleGOdata, algorithm = "classic", statistic = "ks")
     resultKS.elim <- runTest(sampleGOdata, algorithm = "elim", statistic = "ks")
-
-    allRes <- GenTable(sampleGOdata, classicFisher = resultFisher, classicKS = resultKS, elimKS = resultKS.elim, orderBy = "elimKS", ranksOf = "classicFisher", topNodes = 5)
+    allGO = usedGO(sampleGOdata)
+    allRes <- GenTable(sampleGOdata, classicFisher = resultFisher, classicKS = resultKS, elimKS = resultKS.elim, orderBy = "elimKS", ranksOf = "classicFisher", topNodes = length(allGO))
 
     filename = paste("./GO_enrichment_output/contrast_", first_contrast, "_", second_contrast, "_downregulated", ".csv", sep = "")
     write.csv(allRes, file = filename, row.names = F)
 
-    pValue.classic <- score(resultKS)
-    pValue.elim <- score(resultKS.elim)[names(pValue.classic)]
-    gstat <- termStat(sampleGOdata, names(pValue.classic))
-    gSize <- gstat$Annotated / max(gstat$Annotated) * 4
-    gCol <- colMap(gstat$Significant)
-
-    plotfile = paste("./GO_enrichment_output/contrast_", first_contrast, "_", second_contrast, "_downregulated", ".png", sep = "")
-    png(filename=plotfile)
-    plot(pValue.classic, pValue.elim, xlab = "p-value classic", ylab = "p-value elim",
-         pch = 19, cex = gSize, col = gCol, main=paste(first_contrast, second_contrast, "down regulated"))
-    dev.off()
 
-    showSigOfNodes(sampleGOdata, score(resultKS.elim), firstSigNodes = 5, useInfo = 'all')
 
 }