correct readme, add frip score

README.md

@@ -433,7 +433,7 @@ done
 - yes, but must to be the same as in config.yaml (`analysis_dir`) or use relative path
 
 **The pipeline fails because the IDR doesn't select enough reads?**
-- If none or very few peaks pass the IDR, this means that there is to much variability between your replicates aka they are probably not good replicates. If you want to proceed anyway with the analysis, you can use the intersection approach (less stringent) instead of the IDR by setting `intersectionApproach` to 'yes'.
+- If none or very few peaks pass the IDR, this means that there is to much variability between your replicates aka they are probably not good replicates. If you want to proceed anyway with the analysis, you can use the intersection approach (less stringent) instead of the IDR by setting `intersectionApproach:do` to 'yes'.
 
 **When should I use PhantompeakCalling's fragment size estimation instead of MACS2's one?**
 - If MACS2 cannot compute the fragment size estimation (or if you want), set `no_model` to yes, and the fragment length use for MACS2 will be the one computed by PhantompeakQualTools for each sample.

Snakefile

@@ -538,6 +538,7 @@ def IPandINPUT(wildcards):
 
 plotFingerprint_input = IPandINPUT
 plotFingerprint_output = os.path.join(analysis_dir, "05-QC/Fingerprint/{{IP}}_{}_fingerprint.pdf".format(ref))
+plotFingerprint_output_raw =  os.path.join(analysis_dir, "05-QC/Fingerprint/{{IP}}_{}_fingerprint_rawcounts.txt".format(ref))
 plotFingerprint_options = ""
 plotFingerprint_logs = os.path.join(analysis_dir, "05-QC/Fingerprint/logs/{{IP}}_{}_Fingerprint.out".format(ref))
 final_output.extend(expand(plotFingerprint_output,  IP=IP_ALL))    
@@ -569,8 +570,8 @@ if config['geneBody']['do']:
 #------------------------------------------------
 if paired:
     insert_size_metrics_input = "{}/{{SAMPLE}}_{{REF}}_sort{}.bam".format(biasedRegions_dir, biasedRegions)
-    insert_size_metrics_output_txt = os.path.join(analysis_dir, "05-QC/FragmentSizeDistribution/{SAMPLE}_{REF}_sort_dedup_fragmentSizeDistribution.txt")
-    insert_size_metrics_histo = os.path.join(analysis_dir, "05-QC/FragmentSizeDistribution/{SAMPLE}_{REF}_sort_dedup_fragmentSizeDistribution.pdf")
+    insert_size_metrics_output_txt = os.path.join(analysis_dir, "05-QC/FragmentSizeDistribution/{SAMPLE}_{REF}_fragmentSizeDistribution.txt")
+    insert_size_metrics_histo = os.path.join(analysis_dir, "05-QC/FragmentSizeDistribution/{SAMPLE}_{REF}_fragmentSizeDistribution.pdf")
     insert_size_metrics_log_out = os.path.join(analysis_dir, "05-QC/FragmentSizeDistribution/logs/{SAMPLE}_{REF}_fragment_size.out")
     insert_size_metrics_log_err = os.path.join(analysis_dir, "05-QC/FragmentSizeDistribution/logs/{SAMPLE}_{REF}_fragment_size.err")
     final_output.extend(expand(insert_size_metrics_output_txt, SAMPLE=samples, REF=ref))
@@ -738,7 +739,7 @@ def stats_pc_input(wildcards):
         return expand(os.path.join(analysis_dir, "06-PeakCalling/{{CALLER}}/{IP_REP}.{{MOD}}.bed"), IP_REP=IP_ALL)
 
 stats_peakCalling_input = stats_pc_input
-stats_peakCalling_csv = os.path.join(analysis_dir, "{CALLER}_{MOD}Peaks_metrics_mqc.out")
+stats_peakCalling_csv = os.path.join(analysis_dir, "{CALLER}_{MOD}_Peaks_metrics_mqc.out")
 stats_peakCalling_marks = marks
 stats_peakCalling_conds = conds
 stats_peakCalling_rep = rep_flag

config/multiqc_config.yaml

@@ -162,6 +162,7 @@ top_modules:
     - 'cutadapt'
     - 'bowtie2'
     - 'picard'
+    - 'deeptools'
     - 'phantompeakqualtools'
     - 'peaks_metrics'
     - 'idr_metrics'
@@ -200,6 +201,8 @@ extra_fn_clean_exts:
     - _sort
     - _sort_dedup_NoBlacklist
     - _sort_dedup
+    - _R1
+    - _R2
     - type: remove
       pattern: '.sorted'
     - type: regex
@@ -246,12 +249,19 @@ sp:
         fn: '*_spp.out'
     picard/markdups:
         fn: '*_dedup.txt'
+    picard/insertsize:
+        fn: '*_fragmentSizeDistribution.txt'
+    deeptools/plotFingerprintOutRawCounts:
+        fn: '*_fingerprint_rawcounts.txt'
     idr_metrics:
         fn: 'IDR_metrics_mqc.out'
     peaks_metrics:
         fn: 'Peaks_metrics_mqc.out'
     spikes_metrics:
         fn: 'Spikes_metrics_mqc.out'
+    frip_scores:
+        fn: 'frip_metrics_mqc.out'
+
 
 
 # Overwrite the defaults of which table columns are visible by default

workflow/rules/compute_frip.rules

+#########################################################################
+# ePeak: Standardize and reproducible ChIP-seq analysis from raw        #
+#           data to differential analysis                               #
+# Authors: Rachel Legendre, Maelle Daunesse                             #
+# Copyright (c) 2019-2020  Institut Pasteur (Paris) and CNRS.           #
+#                                                                       #
+# This file is part of ePeak workflow.                                  #
+#                                                                       #
+# ePeak is free software: you can redistribute it and/or modify         #
+# it under the terms of the GNU General Public License as published by  #
+# the Free Software Foundation, either version 3 of the License, or     #
+# (at your option) any later version.                                   #
+#                                                                       #
+# ePeak is distributed in the hope that it will be useful,              #
+# but WITHOUT ANY WARRANTY; without even the implied warranty of        #
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the          #
+# GNU General Public License for more details .                         #
+#                                                                       #
+# You should have received a copy of the GNU General Public License     #
+# along with ePeak (LICENSE).                                           #
+# If not, see <https://www.gnu.org/licenses/>.                          #
+#########################################################################
+
+
+rule compute_frip:
+    input:
+        bam = compute_frip_bam,
+        bed = compute_frip_peaks
+    output:
+        compute_frip_csv
+    params:
+        paired = compute_frip_paired
+    singularity:
+        "epeak.img"
+    log:
+        out = compute_frip_log
+    shell:
+        """
+        temp_file=$(mktemp)
+
+        if [[ {params.paired} == "True" ]]
+        then
+            #write reads as bedpe bed file to get fragments instead of reads
+            bedtools bamtobed -bedpe -i {input.bam} > $temp_file 
+        else
+            bedtools bamtobed -i {input.bam} > $temp_file
+        fi
+
+
+        #compute number of fragments overlap peaks
+        mapReads=$(bedtools intersect -a $temp_file -b {input.bed} -wa -u | wc -l)
+
+        #compute number of mapped reads
+        totReads=$(samtools view -c -F 260 {input.bam})
+
+        #compute FRiP
+        totReads/mapReads
+        """

workflow/rules/plotFingerprint.rules

@@ -26,7 +26,8 @@ rule plotFingerprint:
     input:
         plotFingerprint_input
     output:
-        plotFingerprint_output
+        plot = plotFingerprint_output,
+        raw  = plotFingerprint_output_raw
     params:
         options = plotFingerprint_options
     log:
@@ -44,5 +45,5 @@ rule plotFingerprint:
         IP=${{infiles[0]}} ; IP=$(basename ${{IP%_sort*}})
         INPUT=${{infiles[1]}} ; INPUT=$(basename ${{INPUT%_sort*}})
 
-        plotFingerprint --bam {input} --plotFile {output} --labels ${{IP}} ${{INPUT}} {params.options} --numberOfProcessors {threads} 2> {log.out}
+        plotFingerprint --bam {input} --plotFile {output.plot} --outRawCounts {output.raw} --labels ${{IP}} ${{INPUT}} {params.options} --numberOfProcessors {threads} 2> {log.out}
         """