igv session download from report

workflow/report/igv_session.rst

+The IGV session can be downloaded here. It contains all files necessary for the session.
+Unzip the directory 'report_igv_session.zip'.
+`Install IGV <http://software.broadinstitute.org/software/igv/download>`_ and open session by selecting
+'report_igv_session.xml'.
+
+Alternatively you can also start the IGV session directly from the results of the workflow.
+The session file is located in 'results/IGV/igv_session.xml'.

workflow/rules/common.smk

@@ -161,8 +161,8 @@ def get_read_group(wildcards):
         unit=wildcards.unit,
         platform=units.loc[(wildcards.sample, wildcards.unit), "platform"])
 
-# def get_deseq2_igv_input():
-#     return [file for file in os.listdir("results/deseq2/FDR/bed_files/FDR_0.05_{antibody}.consensus_{peak}-peaks")]
+def get_unique_antibodies():
+    return set([antibody for antibody in samples["antibody"] if not pd.isnull(antibody)])
 
 def get_igv_input():
     igv_input = []
@@ -179,8 +179,8 @@ def get_igv_input():
             peak=config["params"]["peak-analysis"]
         )
     )
-    unique_antibodies = set([antibody for antibody in samples["antibody"] if not pd.isnull(antibody)])
-    for antibody in unique_antibodies:
+
+    for antibody in get_unique_antibodies():
         igv_input.extend(
             expand(
                 [
@@ -193,6 +193,31 @@ def get_igv_input():
         )
     return igv_input
 
+def get_files_for_igv():
+    igv_files = []
+    for sample in samples.index:
+        igv_files.extend(
+            expand(
+                [
+                    "results/big_wig/{sample}.bigWig"
+                ],
+                sample=sample
+            )
+        )
+
+    igv_files.extend(
+        expand(
+            [
+                "results/macs2_callpeak/{sample_control_peak}_peaks.{peak}Peak",
+                "results/macs2_merged_expand/{antibody}.consensus_{peak}-peaks.boolean.bed"
+            ],
+            sample_control_peak=get_sample_control_peak_combinations_list(),
+            peak=config["params"]["peak-analysis"],
+            antibody=get_unique_antibodies()
+        )
+    )
+    return igv_files
+
 def get_multiqc_input(wildcards):
     multiqc_input = []
     for (sample, unit) in units.index:
@@ -312,7 +337,8 @@ def all_input(wildcards):
     # QC with fastQC and multiQC, igv session
     wanted_input.extend([
         "results/qc/multiqc/multiqc.html",
-        "results/IGV/igv_session.xml"
+        "results/IGV/igv_session.xml",
+        "results/IGV/report_igv_session.zip"
     ])
 
     # trimming reads
@@ -345,7 +371,6 @@ def all_input(wildcards):
         wanted_input.extend(
             expand (
                 [
-                    # "results/IGV/big_wig/merged_library.bigWig.igv.txt",
                     "results/phantompeakqualtools/{sample}.phantompeak.pdf"
                 ],
                 sample = sample
@@ -396,8 +421,7 @@ def all_input(wildcards):
                                     expand(
                                         [
                                             "results/macs2_merged_expand/{antibody}.consensus_{peak}-peaks.boolean.saf",
-                                            "results/macs2_merged_expand/plots/{antibody}.consensus_{peak}-peaks.boolean.intersect.plot.pdf",
-                                            # "results/IGV/consensus/merged_library.{antibody}.consensus_{peak}-peaks.igv.txt"
+                                            "results/macs2_merged_expand/plots/{antibody}.consensus_{peak}-peaks.boolean.intersect.plot.pdf"
                                         ],
                                         peak = config["params"]["peak-analysis"],
                                         antibody = antibody
@@ -443,7 +467,6 @@ def all_input(wildcards):
                         "results/macs2_callpeak/{sample}-{control}.{peak}_treat_pileup.bdg",
                         "results/macs2_callpeak/{sample}-{control}.{peak}_control_lambda.bdg",
                         "results/macs2_callpeak/{sample}-{control}.{peak}_peaks.{peak}Peak",
-                        # "results/IGV/macs2_callpeak-{peak}/merged_library.{sample}-{control}.{peak}_peaks.igv.txt",
                         "results/macs2_callpeak/plots/plot_{peak}_peaks_count.pdf",
                         "results/macs2_callpeak/plots/plot_{peak}_peaks_frip_score.pdf",
                         "results/macs2_callpeak/plots/plot_{peak}_peaks_macs2.pdf"

workflow/rules/consensus_peak_analysis.smk

@@ -47,7 +47,7 @@ rule create_consensus_bed:
     input:
         "results/macs2_merged_expand/{antibody}.consensus_{peak}-peaks.boolean.txt"
     output:
-        report("results/macs2_merged_expand/{antibody}.consensus_{peak}-peaks.boolean.bed", category="accessory files")
+        "results/macs2_merged_expand/{antibody}.consensus_{peak}-peaks.boolean.bed"
     conda:
         "../envs/gawk.yaml"
     log:
@@ -94,8 +94,6 @@ rule create_consensus_igv:
     shell:
         "find {input} -type f -name '*.consensus_{wildcards.peak}-peaks.boolean.bed' -exec "
         "echo -e '{{}}\t0,0,0' \; > {output} 2> {log}"
-        # "find {input} -type f -name '*.consensus_{wildcards.peak}-peaks.boolean.bed' -exec "
-        # "echo -e 'results/IGV/consensus/{wildcards.antibody}/\"{{}}\"\t0,0,0' \; > {output} 2> {log}"
 
 rule homer_consensus_annotatepeaks:
     input:
@@ -192,10 +190,7 @@ rule featurecounts_deseq2:
         FDR_1_perc_res=directory("results/deseq2/FDR/results/FDR_0.01_{antibody}.consensus_{peak}-peaks"),
         FDR_5_perc_res=directory("results/deseq2/FDR/results/FDR_0.05_{antibody}.consensus_{peak}-peaks"),
         FDR_1_perc_bed=directory("results/deseq2/FDR/bed_files/FDR_0.01_{antibody}.consensus_{peak}-peaks"),
-        FDR_5_perc_bed=report(
-            directory("results/deseq2/FDR/bed_files/FDR_0.05_{antibody}.consensus_{peak}-peaks"),
-            patterns=["{antibody}.{{group_1_vs_group_2}}.FDR_0.05_results.bed"],
-            category="accessory files"),
+        FDR_5_perc_bed=directory("results/deseq2/FDR/bed_files/FDR_0.05_{antibody}.consensus_{peak}-peaks"),
         igv_FDR_5_bed="results/IGV/consensus/merged_library.{antibody}.consensus_{peak}-peaks.deseq2.FDR_0.05.igv.txt",
         plot_FDR_1_perc_MA=report(
             directory("results/deseq2/comparison_plots/MA_plots/FDR_0.01_{antibody}consensus_{peak}-peaks"),

workflow/rules/igv_session.smk

@@ -2,34 +2,73 @@ rule create_igv_input_file:
     input:
         get_igv_input()
     output:
-        # "results/IGV/merged_igv_files.txt"
-        "results/IGV/igv_files.txt"
+        temp("results/IGV/igv_files.txt")
     log:
         "logs/igv/merged_igv_files.log"
     shell:
         "cat {input} > {output} 2> {log}"
 
-# remove paths for igv session download from report.zip
-# rule igv_report_cleanup:
-#     input:
-#         "results/IGV/merged_igv_files.txt"
-#     output:
-#         "results/IGV/igv_files.txt"
-#     log:
-#         "logs/igv/igv_files.log"
-#     script:
-#         " ../scripts/igv_report_cleanup.py"
-
+# igv session that can be started directly from the generated files of workflow
 rule igv_files_to_session:
     input:
-        # igv_data=get_igv_data,
         igv="results/IGV/igv_files.txt",
         fasta="resources/ref/genome.fasta"
     output:
-        report("results/IGV/igv_session.xml", category = "igv session")
+        "results/IGV/igv_session.xml"
     params:
         "--path_prefix '../../'"
     log:
         "logs/igv/igv_session_to_file.log"
     shell:
         " ../workflow/scripts/igv_files_to_session.py {output} {input.igv} ../../{input.fasta} {params} 2> {log}"
+
+# remove paths for igv session to download from report.zip
+rule igv_report_cleanup:
+    input:
+        "results/IGV/igv_files.txt"
+    output:
+        temp("results/IGV/report_igv_files.txt")
+    log:
+        "logs/igv/igv_files.log"
+    script:
+        "../scripts/igv_report_cleanup.py"
+
+rule igv_files_to_report:
+    input:
+        igv="results/IGV/report_igv_files.txt",
+        fasta="resources/ref/genome.fasta"
+    output:
+        temp("results/IGV/report_igv_session.xml")
+    params:
+        ""
+    log:
+        "logs/igv/igv_files_to_report.log"
+    shell:
+        " ../workflow/scripts/igv_files_to_session.py {output} {input.igv} $(basename {input.fasta}) {params} 2> {log}"
+
+rule collect_igv_report_session_files:
+    input:
+        igv_data=get_files_for_igv(),
+        deseq2_files=directory(expand("results/deseq2/FDR/bed_files/FDR_0.05_{antibody}.consensus_{peak}-peaks",
+        antibody=get_unique_antibodies(), peak=config["params"]["peak-analysis"])),
+        fasta="resources/ref/genome.fasta",
+        igv_session="results/IGV/report_igv_session.xml"
+    output:
+        temp(directory("results/IGV/report_igv_session"))
+    log:
+        "logs/igv/collect_igv_report_session_files.log"
+    shell:
+        "mkdir -p {output}; cp {input.igv_data} {output}/; "
+        "cp -R {input.deseq2_files}/* {output}/; cp {input.fasta} {output}/; "
+        "cp {input.igv_session} {output}/"
+
+# igv session that can be downloaded from generated report
+rule zip_igv_report_session:
+    input:
+        directory("results/IGV/report_igv_session")
+    output:
+        report("results/IGV/report_igv_session.zip", caption = "../report/igv_session.rst", category="IGV session")
+    log:
+        "logs/igv/collect_igv_report_session_files.log"
+    shell:
+        "cd $(dirname {input}); zip $(basename {output}) $(basename {input})/*"

workflow/rules/peak_analysis.smk

@@ -52,10 +52,9 @@ rule macs2_callpeak_broad:
                  "_treat_pileup.bdg",
                  "_control_lambda.bdg",
                  # these output extensions internally set the --broad option:
+                 "_peaks.broadPeak",
                  "_peaks.gappedPeak"
-                 ),
-        report("results/macs2_callpeak/{sample}-{control}.broad_peaks.broadPeak", category="accessory files")
-
+                 )
     log:
         "logs/macs2/callpeak.{sample}-{control}.broad.log"
     params:
@@ -82,9 +81,9 @@ rule macs2_callpeak_narrow:
                  "_treat_pileup.bdg",
                  "_control_lambda.bdg",
                  # these output extensions internally set the --broad option:
+                 "_peaks.narrowPeak",
                  "_summits.bed"
-                 ),
-        report("results/macs2_callpeak/{sample}-{control}.narrow_peaks.narrowPeak", category="accessory files")
+                 )
     log:
         "logs/macs2/callpeak.{sample}-{control}.narrow.log"
     params:

workflow/rules/post-analysis.smk

@@ -102,7 +102,7 @@ rule bedGraphToBigWig:
         bedGraph="results/bed_graph/{sample}.sorted.bedgraph",
         chromsizes="resources/ref/genome.chrom.sizes"
     output:
-        report("results/big_wig/{sample}.bigWig", category="accessory files")
+        "results/big_wig/{sample}.bigWig"
     log:
         "logs/big_wig/{sample}.log"
     params:

workflow/scripts/igv_report_cleanup.py

+from smart_open import open
+import os
+
 log = snakemake.log_fmt_shell(stdout=False, stderr=True)
 
-merged_igv = snakemake.input[0]
+with open(snakemake.input[0]) as fin:
+    with open(snakemake.output[0], 'w') as fout:
+        for line in fin:
+            items = line.split("\t")
+            items[0] = os.path.basename(items[0])
+            line = "{path}\t{col}".format(path=items[0], col=items[1])
+            fout.write(line)