Added visualiziation of enriched GO terms.

Snakefile

@@ -42,6 +42,7 @@ onerror:
 rule all:
     """Define desired output from pipeline."""
     input:
+        "Settings/Report.txt",
         #"Samplesize/Report.txt",
         #"Reports/Tools.txt",
         "GO/comparison.txt"

config.yaml

@@ -42,5 +42,4 @@ comparison_settings: ["-ta 0.1 -td 0.1 -qvt 0.25 -brlen 0.7 -cap 1.5 -rx 1 -gene
 
 #Settings for sample size differences
 sizes: [20, 30, 40, 50, 60, 70, 80, 90]
-#sizes: [30, 40]
 repeats: [1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20]

rules/ComparisonSettings.smk

@@ -4,15 +4,21 @@ rule gistic_settings:
         gistic_directory=os.path.join(config["gisticdir"], "gistic2"),
         seg="Input/Segments_tumor.txt"
     output:
-        expand("Settings/GISTIC_{setting_nr}/all_lesions.conf_" + config["gistic_precision"] + ".txt", setting_nr=config["comparison_settings"],
-        expand("Settings/GISTIC_{setting_nr}/regions_track.conf_" + config["gistic_precision"] + ".bed", setting_nr=config["comparison_settings"]
+        expand("Settings/GISTIC_{setting_nr}/all_lesions.conf_" + config["gistic_precision"] + ".txt", setting_nr=range(len(config["comparison_settings"])),
+        expand("Settings/GISTIC_{setting_nr}/regions_track.conf_" + config["gistic_precision"] + ".bed", setting_nr=range(len(config["comparison_settings"]))
     params:
         cnv="",
         ref=config["reference"],
         ref_file="",
-        extra={setting_nr}
+        extra=config["comparison_settings"][{setting_nr}]
     wrapper:
         "file:" + workflow.basedir + "/wrappers/GISTIC2"
 
 rule compare_settings:
-    ##
+    input:
+        "Settings/GISTIC_{setting_nr}/all_lesions.conf_" + config["gistic_precision"] + ".txt"
+    output:
+        "Settings/Report.txt"
+    run:
+        with open(output[0], 'w') as out:
+            out.write(input[0])

wrappers/topgo/environment.yaml

@@ -4,5 +4,5 @@ channels:
 dependencies:
     - bioconductor-topgo =2.30.0
     - bioconductor-biomart =2.32.1 #2.34.1 needed due to error, not 2.34.0
-#    - bioconductor-rgraphviz =2.22.0
+    - bioconductor-rgraphviz =2.22.0
 #    - cairo =1.14.10

wrappers/topgo/wrapper.R

@@ -16,9 +16,7 @@ EG2GO <- read.table(snakemake@params[["gene_ID_to_GO"]], sep="\t", header=TRUE)
 colnames(EG2GO) <- c("ensembl_gene_id","go_id")
 }
 
-geneID2GO <- by(EG2GO$go_id,
-                EG2GO$ensembl_gene_id,
-                function(x) as.character(x))
+geneID2GO <- by(EG2GO$go_id, EG2GO$ensembl_gene_id, function(x) as.character(x))
 geneUniverse <- names(geneID2GO)
 
 #Read input gene list
@@ -31,29 +29,24 @@ names(geneList) <- geneUniverse
 ontologies <- if (snakemake@params[["ontology"]] == "all") c("BP", "CC", "MF") else snakemake@params[["ontology"]]
 
 #Combine data
-GOdata <- new("topGOdata",
-              ontology=ontologies,
-              allGenes = geneList,
-              annot = annFUN.gene2GO,
-              gene2GO = geneID2GO)
+GOdata <- new("topGOdata", ontology=ontologies, allGenes = geneList, annot = annFUN.gene2GO, gene2GO = geneID2GO)
 
 #Perform enrichment tests
 test.stat <- new("classicCount", testStatistic = GOFisherTest, name = "Fisher test")
 resultFisher <- getSigGroups(GOdata, test.stat)
-test.stat <- new("classicScore", testStatistic = GOKSTest, name = "KS tests")
-resultKS <- getSigGroups(GOdata, test.stat)
+#test.stat <- new("classicScore", testStatistic = GOKSTest, name = "KS tests")
+#resultKS <- getSigGroups(GOdata, test.stat)
 test.stat <- new("weightCount", testStatistic = GOFisherTest, name = "Fisher test", sigRatio = "ratio")
 resultWeight <- getSigGroups(GOdata, test.stat)
 
 #Make table with results
 allRes <- GenTable(GOdata, classic = resultFisher,
-                   KS = resultKS, weight = resultWeight,
+                   #KS = resultKS,
+                   weight = resultWeight,
                    orderBy = "weight", ranksOf = "classic", topNodes = 50)
 print(allRes)
-write.table(allRes, file=snakemake@output[[1]], quote=FALSE)
+write.table(allRes, file=snakemake@output[[1]], quote=FALSE, sep="\t", eol = "\n")
 
 #Visualize GO term relationships
-#Rgraphviz nodig!
-
-#showSigOfNodes(myGOdata, score(resultFisher), firstSigNodes = 5, useInfo ='all')
-#printGraph(GOdata, resultFisher, firstSigNodes = 5, fn.prefix = "tGO", useInfo = "all", pdfSW = TRUE)
+showSigOfNodes(GOdata, score(resultFisher), firstSigNodes = 5, useInfo ='all')
+printGraph(GOdata, resultFisher, firstSigNodes = 5, fn.prefix = "tGO", useInfo = "all", pdfSW = TRUE)