create igv session file

workflow/Snakefile

@@ -22,6 +22,7 @@ include: "rules/utils.smk"
 include: "rules/post-analysis.smk"
 include: "rules/peak_analysis.smk"
 include: "rules/consensus_peak_analysis.smk"
+include: "rules/igv_session.smk"
 
 rule all:
     input: all_input

workflow/rules/common.smk

@@ -161,6 +161,38 @@ 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_igv_input():
+    igv_input = []
+    igv_input.extend([
+        "results/IGV/big_wig/merged_library.bigWig.igv.txt"
+    ])
+    igv_input.extend(
+        expand(
+            [
+                "results/IGV/macs2_callpeak-{peak}/merged_library.{sample_control_peak}_peaks.igv.txt",
+
+            ],
+            sample_control_peak=get_sample_control_peak_combinations_list(),
+            peak=config["params"]["peak-analysis"]
+        )
+    )
+    unique_antibodies = set([antibody for antibody in samples["antibody"] if not pd.isnull(antibody)])
+    for antibody in unique_antibodies:
+        igv_input.extend(
+            expand(
+                [
+                    "results/IGV/consensus/merged_library.{antibody}.consensus_{peak}-peaks.igv.txt",
+                    "results/IGV/consensus/merged_library.{antibody}.consensus_{peak}-peaks.deseq2.FDR_0.05.igv.txt"
+                ],
+                antibody=antibody,
+                peak=config["params"]["peak-analysis"]
+            )
+        )
+    return igv_input
+
 def get_multiqc_input(wildcards):
     multiqc_input = []
     for (sample, unit) in units.index:
@@ -277,9 +309,10 @@ def all_input(wildcards):
 
     wanted_input = []
 
-    # QC with fastQC and multiQC
+    # QC with fastQC and multiQC, igv session
     wanted_input.extend([
-        "results/qc/multiqc/multiqc.html"
+        "results/qc/multiqc/multiqc.html",
+        "results/IGV/igv_session.xml"
     ])
 
     # trimming reads
@@ -312,7 +345,7 @@ def all_input(wildcards):
         wanted_input.extend(
             expand (
                 [
-                    "results/IGV/big_wig/merged_library.bigWig.igv.txt",
+                    # "results/IGV/big_wig/merged_library.bigWig.igv.txt",
                     "results/phantompeakqualtools/{sample}.phantompeak.pdf"
                 ],
                 sample = sample
@@ -364,7 +397,7 @@ def all_input(wildcards):
                                         [
                                             "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/IGV/consensus/merged_library.{antibody}.consensus_{peak}-peaks.igv.txt"
                                         ],
                                         peak = config["params"]["peak-analysis"],
                                         antibody = antibody
@@ -410,7 +443,7 @@ 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/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

+import os
+
+from snakemake import rules
+
 rule bedtools_merge_broad:
     input:
         get_macs2_peaks()
@@ -43,7 +47,7 @@ rule create_consensus_bed:
     input:
         "results/macs2_merged_expand/{antibody}.consensus_{peak}-peaks.boolean.txt"
     output:
-        "results/macs2_merged_expand/{antibody}.consensus_{peak}-peaks.boolean.bed"
+        report("results/macs2_merged_expand/{antibody}.consensus_{peak}-peaks.boolean.bed", category="accessory files")
     conda:
         "../envs/gawk.yaml"
     log:
@@ -89,7 +93,9 @@ rule create_consensus_igv:
         "logs/igv/consensus/merged_library.{antibody}.consensus_{peak}-peaks.igv.log"
     shell:
         "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}"
+        "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:
@@ -186,7 +192,11 @@ 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=directory("results/deseq2/FDR/bed_files/FDR_0.05_{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"),
+        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"),
             patterns=["{antibody}.{{group_1_vs_group_2}}.MA-plot_FDR_0.01.pdf"],
@@ -228,14 +238,3 @@ rule featurecounts_deseq2:
         "../envs/featurecounts_deseq2.yaml"
     script:
         "../scripts/featurecounts_deseq2.R"
-
-rule create_deseq2_igv:
-    input:
-        "results/deseq2/results/{antibody}.consensus_{peak}-peaks.deseq2.FDR_0.05.results.bed"
-    output:
-        "results/IGV/consensus/merged_library.{antibody}.consensus_{peak}-peaks.deseq2.FDR_0.05.igv.txt"
-    log:
-        "logs/igv/consensus/merged_library.{antibody}.consensus_{peak}-peaks.deseq2.FDR_0.05.igv.log"
-    shell:
-        "find {input} -type f -name '*.consensus_{wildcards.peak}-peaks.deseq2.FDR_0.05.results.bed' -exec "
-        "echo -e 'results/IGV/consensus/{wildcards.antibody}/deseq2/\"{{}}\"\t255,0,0' \; > {output} 2> {log}"

workflow/rules/igv_session.smk

+rule create_igv_input_file:
+    input:
+        get_igv_input()
+    output:
+        # "results/IGV/merged_igv_files.txt"
+        "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"
+
+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")
+    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}"

workflow/rules/peak_analysis.smk

@@ -52,9 +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"
@@ -82,9 +82,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:
@@ -182,7 +182,7 @@ rule create_igv_peaks:
         "logs/igv/create_igv_peaks/merged_library.{sample}-{control}.{peak}_peaks.log"
     shell:
         " find {input} -type f -name '*_peaks.{wildcards.peak}Peak' -exec "
-        "echo -e 'results/IGV/macs2_callpeak/{wildcards.peak}/\"{{}}\"\t0,0,178' \; > {output} 2> {log}"
+        "echo -e '{{}}\t0,0,178' \; > {output} 2> {log}"
 
 rule homer_annotatepeaks:
     input:

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:
-        "results/big_wig/{sample}.bigWig"
+        report("results/big_wig/{sample}.bigWig", category="accessory files")
     log:
         "logs/big_wig/{sample}.log"
     params:
@@ -120,7 +120,7 @@ rule create_igv_bigwig:
         "logs/igv/create_igv_bigwig.log"
     shell:
         "find {input} -type f -name '*.bigWig' -exec "
-        "echo -e 'results/IGV/big_wig/\"{{}}\"\t0,0,178' \;  > {output} 2> {log}"
+        "echo -e '{{}}\t0,0,178' \;  > {output} 2> {log}"
 
 rule compute_matrix:
     input:

workflow/rules/qc.smk

@@ -16,7 +16,7 @@ rule multiqc:
     input:
         get_multiqc_input
     output:
-        "results/qc/multiqc/multiqc.html"
+        report("results/qc/multiqc/multiqc.html", category="MultiQC report")
     log:
         "logs/multiqc.log"
     wrapper:

workflow/scripts/featurecounts_deseq2.R

@@ -305,6 +305,8 @@ if (file.exists(ResultsFile) == FALSE) {
                     # AVI: dynamically file name extensions added
                     CompResultsFile <- paste(snakemake@output[["FDR_5_perc_res"]], "/", snakemake@params[["antibody"]], ".", CompPrefix, ".FDR_0.05_results.txt", sep="")
                     CompBEDFile <- paste(snakemake@output[["FDR_5_perc_bed"]], "/", snakemake@params[["antibody"]], ".", CompPrefix, ".FDR_0.05_results.bed", sep="")
+                    # AVI: write paths of FDR 0.05 bed files to igv file
+                    cat(paste0(CompBEDFile, "\t255,0,0\n"),file=snakemake@output[["igv_FDR_5_bed"]],append=TRUE)
                     MAplotFile <- paste(snakemake@output[["plot_FDR_5_perc_MA"]], "/", snakemake@params[["antibody"]], ".", CompPrefix, ".MA-plot_FDR_0.05.pdf", sep="")   # AVI: added to create separate pdf files
                     VolcanoPlotFile <- paste(snakemake@output[["plot_FDR_5_perc_volcano"]], "/", snakemake@params[["antibody"]], ".", CompPrefix, ".volcano_FDR_0.05.pdf", sep="")  # AVI: added to create separate pdf files
                 }

workflow/scripts/igv_files_to_session.py

+#!/usr/bin/env python
+
+### source: https://github.com/nf-core/chipseq/blob/master/bin/igv_files_to_session.py
+### own changes and adjustments for snakemake-workflow chipseq are marked with "# AVI: " in comment
+
+#MIT License
+#
+#Copyright (c) Philip Ewels
+#
+#Permission is hereby granted, free of charge, to any person obtaining a copy
+#of this software and associated documentation files (the "Software"), to deal
+#in the Software without restriction, including without limitation the rights
+#to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
+#copies of the Software, and to permit persons to whom the Software is
+#furnished to do so, subject to the following conditions:
+#
+#The above copyright notice and this permission notice shall be included in all
+#copies or substantial portions of the Software.
+#
+#THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
+#IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+#FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
+#AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+#LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
+#OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
+#SOFTWARE.
+
+
+
+#######################################################################
+#######################################################################
+## Created on July 4th 2018 to create IGV session file from file list
+#######################################################################
+#######################################################################
+
+import os
+import errno
+import argparse
+
+############################################
+############################################
+## PARSE ARGUMENTS
+############################################
+############################################
+
+Description = 'Create IGV session file from a list of files and associated colours - ".bed", ".bw", ".bigwig", ".tdf", ".gtf" files currently supported.'
+Epilog = """Example usage: python igv_files_to_session.py <XML_OUT> <LIST_FILE> <GENOME>"""
+
+argParser = argparse.ArgumentParser(description=Description, epilog=Epilog)
+
+## REQUIRED PARAMETERS
+argParser.add_argument('XML_OUT', help="XML output file.")
+argParser.add_argument('LIST_FILE', help="Tab-delimited file containing two columns i.e. file_name\tcolour. Header isnt required.")
+argParser.add_argument('GENOME', help="Full path to genome fasta file or shorthand for genome available in IGV e.g. hg19.")
+
+## OPTIONAL PARAMETERS
+argParser.add_argument('-pp', '--path_prefix', type=str, dest="PATH_PREFIX", default='', help="Path prefix to be added at beginning of all files in input list file.")
+args = argParser.parse_args()
+
+############################################
+############################################
+## HELPER FUNCTIONS
+############################################
+############################################
+
+def makedir(path):
+
+    if not len(path) == 0:
+        try:
+            os.makedirs(path)
+        except OSError as exception:
+            if exception.errno != errno.EEXIST:
+                raise
+
+############################################
+############################################
+## MAIN FUNCTION
+############################################
+############################################
+
+def igv_files_to_session(XMLOut,ListFile,Genome,PathPrefix=''):
+
+    makedir(os.path.dirname(XMLOut))
+
+    fileList = []
+    fin = open(ListFile,'r')
+    while True:
+        line = fin.readline()
+        if line:
+            ifile,colour = line.strip().split('\t')
+            if len(colour.strip()) == 0:
+                colour = '0,0,178'
+            fileList.append((PathPrefix.strip()+ifile,colour))
+        else:
+            break
+            fout.close()
+
+    ## ADD RESOURCES SECTION
+    XMLStr =  '<?xml version="1.0" encoding="UTF-8" standalone="no"?>\n'
+    XMLStr += '<Session genome="%s" hasGeneTrack="true" hasSequenceTrack="true" locus="All" version="8">\n' % (Genome)
+    XMLStr += '\t<Resources>\n'
+    for ifile,colour in fileList:
+        XMLStr += '\t\t<Resource path="%s"/>\n' % (ifile)
+    XMLStr += '\t</Resources>\n'
+
+    ## ADD PANEL SECTION
+    XMLStr += '\t<Panel height="1160" name="DataPanel" width="1897">\n'
+    for ifile,colour in fileList:
+        extension = os.path.splitext(ifile)[1].lower()
+        if extension in ['.bed','.broadpeak','.narrowpeak']:
+            XMLStr += '\t\t<Track altColor="0,0,178" autoScale="false" clazz="org.broad.igv.track.FeatureTrack" color="%s" ' % (colour)
+            XMLStr += 'displayMode="SQUISHED" featureVisibilityWindow="-1" fontSize="10" height="20" '
+            XMLStr += 'id="%s" name="%s" renderer="BASIC_FEATURE" sortable="false" visible="true" windowFunction="count"/>\n' % (ifile,os.path.basename(ifile))
+        elif extension in ['.bw', '.bigwig', '.tdf']:
+            XMLStr += '\t\t<Track altColor="0,0,178" autoScale="true" clazz="org.broad.igv.track.DataSourceTrack" color="%s" ' % (colour)
+            XMLStr += 'displayMode="COLLAPSED" featureVisibilityWindow="-1" fontSize="10" height="30" '
+            XMLStr += 'id="%s" name="%s" normalize="false" renderer="BAR_CHART" sortable="true" visible="true" windowFunction="mean">\n' % (ifile,os.path.basename(ifile))
+            XMLStr += '\t\t\t<DataRange baseline="0.0" drawBaseline="true" flipAxis="false" maximum="10" minimum="0.0" type="LINEAR"/>\n'
+            XMLStr += '\t\t</Track>\n'
+        elif extension in ['.gtf']:
+            XMLStr += '\t\t<Track altColor="0,0,178" autoScale="false" clazz="org.broad.igv.track.FeatureTrack" color="%s" ' % (colour)
+            XMLStr += 'displayMode="COLLAPSED" featureVisibilityWindow="-1" fontSize="10" '
+            XMLStr += 'id="%s" name="%s" renderer="BASIC_FEATURE" sortable="false" visible="true" windowFunction="count"/>\n' % (ifile,os.path.basename(ifile))
+        elif extension in ['.bam']:
+            pass
+        else:
+            XMLStr += '\t\t<Track altColor="0,0,178" autoScale="false" clazz="org.broad.igv.track.FeatureTrack" color="%s" ' % (colour)
+            XMLStr += 'displayMode="SQUISHED" featureVisibilityWindow="-1" fontSize="10" height="20" '
+            XMLStr += 'id="%s" name="%s" renderer="BASIC_FEATURE" sortable="false" visible="true" windowFunction="count"/>\n' % (ifile,os.path.basename(ifile))
+
+    XMLStr += '\t</Panel>\n'
+    #XMLStr += '\t<HiddenAttributes>\n\t\t<Attribute name="DATA FILE"/>\n\t\t<Attribute name="DATA TYPE"/>\n\t\t<Attribute name="NAME"/>\n\t</HiddenAttributes>\n'
+    XMLStr += '</Session>'
+    XMLOut = open(XMLOut,'w')
+    XMLOut.write(XMLStr)
+    XMLOut.close()
+
+############################################
+############################################
+## RUN FUNCTION
+############################################
+############################################
+
+igv_files_to_session(XMLOut=args.XML_OUT,ListFile=args.LIST_FILE,Genome=args.GENOME,PathPrefix=args.PATH_PREFIX)
+
+############################################
+############################################
+############################################
+############################################

workflow/scripts/igv_report_cleanup.py

+log = snakemake.log_fmt_shell(stdout=False, stderr=True)
+
+merged_igv = snakemake.input[0]