Adding peak-calling with Genrich (fix #108) and multi-mapping read analysis

assets/multiqc/merged_library_frip_score_header.txt

@@ -1,5 +1,5 @@
 #id: 'mlib_frip_score'
-#section_name: 'MERGED LIB: MACS2 peak FRiP score'
+#section_name: 'MERGED LIB: MACS2: Peak FRiP score'
 #description: "is generated by calculating the fraction of all mapped reads that fall
 #              into the MACS2 called peak regions. A read must overlap a peak by at least 20% to be counted.
 #              See <a href='https://www.encodeproject.org/data-standards/terms/' target='_blank'>FRiP score</a>."

assets/multiqc/merged_library_joint_genrich_frip_score_header.txt

@@ -0,0 +1,14 @@
+#id: 'mlib_joint_genrich_frip_score'
+#section_name: 'MERGED LIB: Genrich (joint): Peak FRiP score'
+#description: "is generated by calculating the fraction of all mapped reads that fall
+#              into the Genrich called peak regions. A read must overlap a peak by at least 20% to be counted.
+#              See <a href='https://www.encodeproject.org/data-standards/terms/' target='_blank'>FRiP score</a>."
+#plot_type: 'bargraph'
+#anchor: 'mlib_frip_score'
+#pconfig:
+#    title: 'FRiP score'
+#    ylab: 'FRiP score'
+#    ymax: 1
+#    ymin: 0
+#    tt_decimals: 2
+#    cpswitch: False

assets/multiqc/merged_library_joint_genrich_peak_annotation_header.txt

@@ -0,0 +1,9 @@
+#id: 'mlib_joint_genrich_peak_annotation'
+#section_name: 'MERGED LIB: Genrich (joint): HOMER peak annotation'
+#description: "is generated by calculating the proportion of peaks assigned to genomic features by
+#              <a href='http://homer.ucsd.edu/homer/ngs/annotation.html' target='_blank'>HOMER annotatePeaks.pl</a>."
+#plot_type: 'bargraph'
+#anchor: 'mlib_peak_annotation'
+#pconfig:
+#    title: 'Peak to feature proportion'
+#    ylab: 'Peak count'

assets/multiqc/merged_library_joint_genrich_peak_count_header.txt

@@ -0,0 +1,10 @@
+#id: 'mlib_joint_genrich_peak_count'
+#section_name: 'MERGED LIB: Genrich (joint): Peak count'
+#description: "is calculated from total number of peaks called by
+#	       <a href='https://github.com/jsh58/Genrich' target='_blank'>Genrich</a>"
+#plot_type: 'bargraph'
+#anchor: 'mlib_peak_count'
+#pconfig:
+#    title: 'Total peak count'
+#    ylab: 'Peak count'
+#    cpswitch: False

assets/multiqc/merged_library_peak_annotation_header.txt

@@ -1,5 +1,5 @@
 #id: 'mlib_peak_annotation'
-#section_name: 'MERGED LIB: HOMER peak annotation'
+#section_name: 'MERGED LIB: MACS2: HOMER peak annotation'
 #description: "is generated by calculating the proportion of peaks assigned to genomic features by
 #              <a href='http://homer.ucsd.edu/homer/ngs/annotation.html' target='_blank'>HOMER annotatePeaks.pl</a>."
 #plot_type: 'bargraph'

assets/multiqc/merged_library_peak_count_header.txt

@@ -1,5 +1,5 @@
 #id: 'mlib_peak_count'
-#section_name: 'MERGED LIB: MACS2 peak count'
+#section_name: 'MERGED LIB: MACS2: Peak count'
 #description: "is calculated from total number of peaks called by
 #	       <a href='https://github.com/taoliu/MACS' target='_blank'>MACS2</a>"
 #plot_type: 'bargraph'

assets/multiqc/merged_library_sep_genrich_frip_score_header.txt

@@ -0,0 +1,14 @@
+#id: 'mlib_sep_genrich_frip_score'
+#section_name: 'MERGED LIB: Genrich (sep): Peak FRiP score'
+#description: "is generated by calculating the fraction of all mapped reads that fall
+#              into the Genrich called peak regions. A read must overlap a peak by at least 20% to be counted.
+#              See <a href='https://www.encodeproject.org/data-standards/terms/' target='_blank'>FRiP score</a>."
+#plot_type: 'bargraph'
+#anchor: 'mlib_frip_score'
+#pconfig:
+#    title: 'FRiP score'
+#    ylab: 'FRiP score'
+#    ymax: 1
+#    ymin: 0
+#    tt_decimals: 2
+#    cpswitch: False

assets/multiqc/merged_library_sep_genrich_peak_annotation_header.txt

@@ -0,0 +1,9 @@
+#id: 'mlib_sep_genrich_peak_annotation'
+#section_name: 'MERGED LIB: Genrich (sep): HOMER peak annotation'
+#description: "is generated by calculating the proportion of peaks assigned to genomic features by
+#              <a href='http://homer.ucsd.edu/homer/ngs/annotation.html' target='_blank'>HOMER annotatePeaks.pl</a>."
+#plot_type: 'bargraph'
+#anchor: 'mlib_peak_annotation'
+#pconfig:
+#    title: 'Peak to feature proportion'
+#    ylab: 'Peak count'

assets/multiqc/merged_library_sep_genrich_peak_count_header.txt

@@ -0,0 +1,10 @@
+#id: 'mlib_sep_genrich_peak_count'
+#section_name: 'MERGED LIB: Genrich (sep): Peak count'
+#description: "is calculated from total number of peaks called by
+#	       <a href='https://github.com/jsh58/Genrich' target='_blank'>Genrich</a>"
+#plot_type: 'bargraph'
+#anchor: 'mlib_peak_count'
+#pconfig:
+#    title: 'Total peak count'
+#    ylab: 'Peak count'
+#    cpswitch: False

assets/multiqc/merged_replicate_frip_score_header.txt

@@ -1,5 +1,5 @@
 #id: 'mrep_frip_score'
-#section_name: 'MERGED REP: MACS2 peak FRiP score'
+#section_name: 'MERGED REP: MACS: Peak FRiP score'
 #description: "is generated by calculating the fraction of all mapped reads that fall
 #              into the MACS2 called peak regions. A read must overlap a peak by at least 20% to be counted.
 #              See <a href='https://www.encodeproject.org/data-standards/terms/' target='_blank'>FRiP score</a>."

assets/multiqc/merged_replicate_peak_annotation_header.txt

@@ -1,5 +1,5 @@
 #id: 'mrep_peak_annotation'
-#section_name: 'MERGED REP: HOMER peak annotation'
+#section_name: 'MERGED REP: MACS2: HOMER peak annotation'
 #description: "is generated by calculating the proportion of peaks assigned to genomic features by
 #              <a href='http://homer.ucsd.edu/homer/ngs/annotation.html' target='_blank'>HOMER annotatePeaks.pl</a>."
 #plot_type: 'bargraph'

assets/multiqc/merged_replicate_peak_count_header.txt

@@ -1,5 +1,5 @@
 #id: 'mrep_peak_count'
-#section_name: 'MERGED REP: MACS2 Peak count'
+#section_name: 'MERGED REP: MACS2: Peak count'
 #description: "is calculated from total number of peaks called by
 #	       <a href='https://github.com/taoliu/MACS' target='_blank'>MACS2</a>"
 #plot_type: 'bargraph'

assets/multiqc_config.yml

@@ -91,8 +91,14 @@ module_order:
         - "./macs2/merged_replicate/featurecounts/*.summary"
 
 report_section_order:
-  mlib_peak_count:
+  mlib_sep_genrich_peak_count:
     before: mlib_deeptools
+  mlib_sep_genrich_frip_score:
+    before: mlib_sep_genrich_peak_count
+  mlib_sep_genrich_peak_annotation:
+    before: mlib_sep_genrich_frip_score
+  mlib_peak_count:
+    before: mlib_sep_genrich_peak_annotation
   mlib_frip_score:
     before: mlib_peak_count
   mlib_peak_annotation:
@@ -115,6 +121,12 @@ report_section_order:
     before: mrep_featurecounts
   mrep_deseq2_clustering_1:
     before: mrep_deseq2_pca_1
+  mlib_joint_genrich_peak_count:
+    before: mrep_deseq2_clustering_1
+  mlib_joint_genrich_frip_score:
+    before: mlib_joint_genrich_peak_count
+  mlib_joint_genrich_peak_annotation:
+    before: mlib_joint_genrich_frip_score
   "nf-core-atacseq-methods-description":
     order: -1000
   software_versions:

bin/plot_homer_annotatepeaks.r

@@ -90,6 +90,7 @@ for (idx in 1:length(HomerFiles)) {
     dist.freq[which(unique.gene.dat$Distance.to.TSS < 10000)] <- "< 10kb"
     dist.freq[which(unique.gene.dat$Distance.to.TSS < 5000)] <- "< 5kb"
     dist.freq[which(unique.gene.dat$Distance.to.TSS < 2000)] <- "< 2kb"
+    dist.freq[which(unique.gene.dat$Distance.to.TSS < 1000)] <- "< 1kb"
     dist.freq <- as.data.frame(table(dist.freq))
     colnames(dist.freq) <- c("distance",sampleid)
     dist.melt <- melt(dist.freq)

bin/plot_macs2_qc.r → bin/plot_peaks_qc.r

@@ -20,7 +20,7 @@ library(scales)
 option_list <- list(make_option(c("-i", "--peak_files"), type="character", default=NULL, help="Comma-separated list of peak files.", metavar="path"),
                     make_option(c("-s", "--sample_ids"), type="character", default=NULL, help="Comma-separated list of sample ids associated with peak files. Must be unique and in same order as peaks files input.", metavar="string"),
                     make_option(c("-o", "--outdir"), type="character", default='./', help="Output directory", metavar="path"),
-                    make_option(c("-p", "--outprefix"), type="character", default='macs2_peakqc', help="Output prefix", metavar="string"))
+                    make_option(c("-p", "--outprefix"), type="character", default='_peakqc', help="Output prefix", metavar="string"))
 
 opt_parser <- OptionParser(option_list=option_list)
 opt <- parse_args(opt_parser)

conf/base.config

@@ -47,6 +47,11 @@ process {
     withLabel:process_long {
         time   = { check_max( 20.h  * task.attempt, 'time'    ) }
     }
+    withLabel:process_ultra {
+        cpus   = { check_max( 8    * task.attempt, 'cpus'    ) }
+        memory = { check_max( 150.GB * task.attempt, 'memory'  ) }
+        time   = { check_max( 120.h  * task.attempt, 'time'    ) }
+    }
     withLabel:process_high_memory {
         memory = { check_max( 200.GB * task.attempt, 'memory' ) }
     }