Updated circos gene plots.

rules/Circos.smk

-from Circos import InputCircos, bed_to_circos
+from Circos import InputCircos, bed_to_circos, make_CIRCOS_legend
+from PIL import Image
 
 rule make_CIRCOS_input:
-    """Make input files for making a circos diagram."""
+    """Make input files for making a CIRCOS plot."""
     input:
         seg="Input/Segments_tumor.txt",
         gistic="GISTIC/all_lesions.conf_" + config["gistic_precision"] + ".txt",
@@ -15,7 +16,7 @@ rule make_CIRCOS_input:
         InputCircos(input.seg, input.gistic, input.rubic_gains, input.rubic_losses, output.seg, output.gistic, output.rubic)
 
 rule make_CIRCOS_plot:
-    """Make circos diagram of recurrent regions in RUBIC and GISTIC2.0"""
+    """Make CIRCOS plot of recurrent regions in RUBIC and GISTIC2.0"""
     input:
         seg="Reports/Circos/Segments.txt",
         gistic="Reports/Circos/GISTIC_results.txt",
@@ -29,39 +30,17 @@ rule make_CIRCOS_plot:
     shell:
         "circos -conf {params.conf} -outputfile {output[0]} -param cnv_file={input.seg} -param gistic_file={input.gistic} -param rubic_file={input.rubic}"
 
-rule make_legend_CIRCOS:
-    output:
-        "Reports/Circos/legend.png"
-    run:
-        import matplotlib.pyplot as plt
-        import matplotlib.patches as mpatches
-        gains = mpatches.Patch(color='#74C476', label='Gains')
-        losses = mpatches.Patch(color='#FB6A4A', label='Losses')
-        gistic = mpatches.Patch(color='#5975A4', label='GISTIC2.0 regions')
-        rubic = mpatches.Patch(color='#5F9E6E', label='RUBIC regions')
-        legend = plt.legend(handles=[gains, losses, gistic, rubic], loc=3, framealpha=1, frameon=False)
-        fig = legend.figure
-        fig.canvas.draw()
-        bbox = legend.get_window_extent().transformed(fig.dpi_scale_trans.inverted())
-        fig.savefig(output[0], dpi=400, bbox_inches=bbox)
-
 rule add_legend_CIRCOS:
+    """Add a custom legend to the CIRCOS plot."""
     input:
-        "Reports/Circos/RecurrentRegions.png",
-        "Reports/Circos/legend.png"
+        circos="Reports/Circos/RecurrentRegions.png",
     output:
-        "Reports/Circos/RecurrentRegions_legend.png"
+        legend="Reports/Circos/legend.png",
+        circos="Reports/Circos/RecurrentRegions_legend.png"
     run:
-        from PIL import Image
-        circos = Image.open(input[0], 'r')
-        c_w, c_h = circos.size
-        legend = Image.open(input[1], 'r')
-        l_w, l_h = legend.size
-        offset = ((c_w - l_w), 0)
-        circos.paste(legend, offset)
-        circos.save(output[0])
+        make_CIRCOS_legend(input.circos, output.legend, output.circos)
 
-rule circos_input_zoom: #necessary?
+rule make_CIRCOS_zoom_input: #necessary?
     """Make input files for making a circos diagram."""
     input:
         bed="Reports/Overlap_known_genes.bed"
@@ -72,24 +51,54 @@ rule circos_input_zoom: #necessary?
     run:
         bed_to_circos(input.bed, output.rubic, output.gistic, output.genes)
 
-def get_list_genes(gene_file):
-    with open(gene_file, 'r') as genes:
-        for line in genes:
-            print(line)
+def get_list_genes(overlapping_genes, locations_known_genes):
+    plot_list = []
+    list_overlapping_genes = []
+    with open(overlapping_genes, 'r') as plot_genes:
+        for line in plot_genes:
+            chrom, start, end = line.strip().split("\t")
+            chrom = "chr" + chrom.strip("hs")
+            list_overlapping_genes.append([chrom, start, end])
+    with open(locations_known_genes, 'r') as known_genes:
+        known_genes.readline()
+        for line in known_genes:
+            chrom, start, end, gene_name = line.strip().split("\t")
+            if [chrom, start, end] in list_overlapping_genes:
+                plot_list.append(gene_name)
+    return(plot_list)
+
+def get_plot_region(gene_name, locations_known_genes):
+    with open(locations_known_genes, 'r') as known_genes:
+        known_genes.readline()
+        for line in known_genes:
+            if gene_name == line.strip().split("\t")[3]:
+                chrom, start, end, name = line.strip().split("\t")
+                chrom_region = "hs" + chrom.strip("chr")
+                if int(start) < 1000000:
+                    start_region = 0
+                else:
+                    start_region = int(start) - 1000000
+                end_region = int(end) + 1000000
+                plot_region = chrom_region + ":" + str(start_region) + "-" + str(end_region)
+    return plot_region
 
-rule make_zoomed_CIRCOS_plot:
+rule make_CIRCOS_zoom_plots:
     """Compare locations of known genes, recurrent regions from RUBIC and recurrent regions from GISTIC2."""
     input:
-        gistic="Reports/Circos/Zoom/GISTIC.txt",
-        rubic="Reports/Circos/Zoom/RUBIC.txt",
-        genes="Reports/Circos/Zoom/Genes.txt"
-    output:
-        plots=expand("Reports/Overlap_plots/{gene}.png", gene="NR6") #get_list_genes(input.genes))
+        gistic="Reports/Circos/Zoom/GISTIC.txt",    #all gistic regions
+        rubic="Reports/Circos/Zoom/RUBIC.txt",      #all rubic regions
+        genes="Reports/Circos/Zoom/Genes.txt",
+        known="Reports/Locations_known_genes.bed",   #use locations known genes to extract list of known genes and run rule for each gene    maybe also not overlapping, but region around gene
+        list_genes = lambda
     params:
-        workflow.basedir + "/scripts/circos/circos_zoom.conf",
-        chrom='hs12:67-71'
+        known="Reports/Locations_known_genes.bed",
+        genes="Reports/Circos/Zoom/Genes.txt",
+        conf=workflow.basedir + "/scripts/circos/circos_zoom.conf",
+        chrom=get_plot_region(wildcards.gene, input.known)
+    output:
+        plot=expand("Reports/Overlap_plots/{gene}.png", gene=get_list_genes(params.genes, params.known))
     conda:
         workflow.basedir + "/envs/circos.yaml"
     shell:
-        "circos -conf {params[0]} -outputfile {output[0]} -param gistic_file={input.gistic} -param rubic_file={input.rubic} \
+        "circos -conf {params.conf} -outputfile {output.plot} -param gistic_file={input.gistic} -param rubic_file={input.rubic} \
          -param gene_file={input.genes} -param chrom={params.chrom}"

rules/GenePrioritization.smk

@@ -3,7 +3,7 @@ rule do_GO_analysis:
     input:
         gene_list="Reports/Genes_{tool}.txt"
     output:
-        table="GO/Enriched_GOs_{tool}.txt",
+        table="GO/Enriched_GOs_{tool}.txt",                 #Add gene names to table
         plot="GO/Enriched_GOs_{tool}.jpg"
     params:
         organism="hsapiens",    #default is human
@@ -30,7 +30,7 @@ def compare_go(GO_files, output_file):
     go_terms = get_go(GO_files)
     overlap = set(go_terms[0]).intersection(set(go_terms[1]))
     with open(output_file, 'w') as out:
-        out.write("Percentage overlapping: " + str(int((len(overlap) / 50 * 100))) + "\n")
+        out.write("Percentage overlapping: " + str(int((len(overlap) / 50 * 100))) + "\n\n")
         out.write("\n".join(list(overlap)))
 
 def get_go(GO_files):
@@ -41,8 +41,7 @@ def get_go(GO_files):
         with open(tool_file, 'r') as tool:
             tool.readline()
             for line in tool:
-                go_term = line.split(' ')[1].strip('"')
+                go_term = line.split("\t")[1]
                 go_terms_tool += [go_term]
         go_terms.append(go_terms_tool)
-        go_terms.append(go_terms_tool)
     return go_terms

scripts/Circos.py

 from math import log10
 import os.path
+from PIL import Image
+import matplotlib.pyplot as plt
+import matplotlib.patches as mpatches
 
 class InputCircos:
     """Make input files for a Circos plot."""
@@ -46,7 +49,33 @@ class InputCircos:
                         out_line = ["hs" + chrom, start, end, qval]
                         out.write(" ".join(out_line) + "\n")
 
+
+def make_CIRCOS_legend(CIRCOS_png, legend_png, concatted_png):
+    """Make a legend to add to CIRCOS plot."""
+    gains = mpatches.Patch(color='#74C476', label='Gains')
+    losses = mpatches.Patch(color='#FB6A4A', label='Losses')
+    gistic = mpatches.Patch(color='#5975A4', label='GISTIC2.0 regions')
+    rubic = mpatches.Patch(color='#5F9E6E', label='RUBIC regions')
+    legend = plt.legend(handles=[gains, losses, gistic, rubic], loc=3, framealpha=1, frameon=False)
+    fig = legend.figure
+    fig.canvas.draw()
+    bbox = legend.get_window_extent().transformed(fig.dpi_scale_trans.inverted())
+    fig.savefig(legend_png, dpi=400, bbox_inches=bbox)
+    concat_CIRCOS_legend(CIRCOS_png, legend_png, concatted_png)
+
+def concat_CIRCOS_legend(CIRCOS_png, legend_png, concatted_png):
+    """Combine legend png file and CIRCOS plot."""
+    circos = Image.open(CIRCOS_png, 'r')
+    c_w, c_h = circos.size
+    legend = Image.open(legend_png, 'r')
+    l_w, l_h = legend.size
+    offset = ((c_w - l_w), 0)
+    circos.paste(legend, offset)
+    circos.save(concatted_png)
+
+
 def bed_to_circos(bed_file, rubic_file, gistic_file, gene_file):
+    """Convert bed file to CIRCOS input files."""
     rubic = open(rubic_file, 'w')
     gistic = open(gistic_file, 'w')
     genes = open(gene_file, 'w')

wrappers/topgo/wrapper.R

@@ -30,6 +30,7 @@ ontologies <- if (snakemake@params[["ontology"]] == "all") c("BP", "CC", "MF") e
 
 #Combine data
 GOdata <- new("topGOdata", ontology=ontologies, allGenes = geneList, annot = annFUN.gene2GO, gene2GO = geneID2GO)
+pint(GOdata)
 
 #Perform enrichment tests
 test.stat <- new("classicCount", testStatistic = GOFisherTest, name = "Fisher test")