Fix atacseq#107

CHANGELOG.md

@@ -15,12 +15,15 @@ and this project adheres to [Semantic Versioning](http://semver.org/spec/v2.0.0.
 * [nf-core/atacseq#76](https://github.com/nf-core/atacseq/issues/76) - featureCounts coupled to DESeq2
 * [nf-core/atacseq#79](https://github.com/nf-core/atacseq/issues/79) - Parallelize DESeq2
 * [nf-core/atacseq#97](https://github.com/nf-core/atacseq/issues/97) - PBC1, PBC2 from pipeline?
+* [nf-core/atacseq#107](https://github.com/nf-core/atacseq/issues/107) - Add options to change MACS2 parameters
 * [nf-core/atacseq#109](https://github.com/nf-core/atacseq/issues/109) - Specify custom gtf but gene bed is not generated from that gtf?
 * Regenerated screenshots and added collapsible sections for output files in `docs/output.md`
 * Update template to tools `1.9`
 * Replace `set` with `tuple` and `file()` with `path()` in all processes
 * Capitalise process names
 * Parameters:
+    * `--macs_fdr` to provide FDR threshold for MACS2 peak calling
+    * `--macs_pvalue` to provide p-value threshold for MACS2 peak calling
     * `--skip_peak_qc` to skip MACS2 peak QC plot generation
     * `--skip_peak_annotation` to skip annotation of MACS2 and consensus peaks with HOMER
     * `--skip_consensus_peaks` to skip consensus peak generation

docs/usage.md

@@ -35,6 +35,8 @@
 * [Peaks](#peaks)
     * [`--narrow_peak`](#--narrow_peak)
     * [`--broad_cutoff`](#--broad_cutoff)
+    * [`--macs_fdr`](#--macs_fdr)
+    * [`--macs_pvalue`](#--macs_pvalue)
     * [`--min_reps_consensus`](#--min_reps_consensus)
     * [`--save_macs_pileup`](#--save_macs_pileup)
     * [`--skip_peak_qc`](#--skip_peak_qc)
@@ -389,6 +391,14 @@ MACS2 is run by default with the [`--broad`](https://github.com/taoliu/MACS#--br
 
 Specifies broad cut-off value for MACS2. Only used when `--narrow_peak` isnt specified (Default: `0.1`).
 
+### `--macs_fdr`
+
+Minimum FDR (q-value) cutoff for peak detection, `--macs_fdr` and `--macs_pvalue` are mutually exclusive (Default: false).
+
+### `--macs_pvalue`
+
+p-value cutoff for peak detection, `--macs_fdr` and `--macs_pvalue` are mutually exclusive (Default: false). If `--macs_pvalue` cutoff is set, q-value will not be calculated and reported as -1 in the final .xls file.
+
 ### `--min_reps_consensus`
 
 Number of biological replicates required from a given condition for a peak to contribute to a consensus peak . If you are confident you have good reproducibility amongst your replicates then you can increase the value of this parameter to create a "reproducible" set of consensus of peaks. For example, a value of 2 will mean peaks that have been called in at least 2 replicates will contribute to the consensus set of peaks, and as such peaks that are unique to a given replicate will be discarded.

main.nf

@@ -56,6 +56,8 @@ def helpMessage() {
     Peaks
       --narrow_peak [bool]            Run MACS2 in narrowPeak mode (Default: false)
       --broad_cutoff [float]          Specifies broad cutoff value for MACS2. Only used when --narrow_peak isnt specified (Default: 0.1)
+      --macs_fdr [float]              Minimum FDR (q-value) cutoff for peak detection, --macs_fdr and --macs_pvalue are mutually exclusive (Default: false)
+      --macs_pvalue [float]           p-value cutoff for peak detection, --macs_fdr and --macs_pvalue are mutually exclusive (Default: false)
       --min_reps_consensus [int]      Number of biological replicates required from a given condition for a peak to contribute to a consensus peak (Default: 1)
       --save_macs_pileup [bool]       Instruct MACS2 to create bedGraph files normalised to signal per million reads (Default: false)
       --skip_peak_qc [bool]           Skip MACS2 peak QC plot generation (Default: false)
@@ -230,6 +232,8 @@ summary['MACS2 Genome Size']      = params.macs_gsize ?: 'Not supplied'
 summary['Min Consensus Reps']     = params.min_reps_consensus
 if (params.macs_gsize)            summary['MACS2 Narrow Peaks'] = params.narrow_peak ? 'Yes' : 'No'
 if (!params.narrow_peak)          summary['MACS2 Broad Cutoff'] = params.broad_cutoff
+if (params.macs_fdr)              summary['MACS2 FDR'] = params.macs_fdr
+if (params.macs_pvalue)           summary['MACS2 P-value'] = params.macs_pvalue
 if (params.skip_trimming) {
     summary['Trimming Step']      = 'Skipped'
 } else {
@@ -1121,6 +1125,8 @@ process MACS2 {
     broad = params.narrow_peak ? '' : "--broad --broad-cutoff ${params.broad_cutoff}"
     format = params.single_end ? 'BAM' : 'BAMPE'
     pileup = params.save_macs_pileup ? '-B --SPMR' : ''
+    fdr = params.macs_fdr ? "--qvalue ${params.macs_fdr}" : ''
+    pvalue = params.macs_pvalue ? "--pvalue ${params.macs_pvalue}" : ''
     """
     macs2 callpeak \\
         -t ${ipbam[0]} \\
@@ -1130,6 +1136,8 @@ process MACS2 {
         -g $params.macs_gsize \\
         -n $ip \\
         $pileup \\
+        $fdr \\
+        $pvalue \\
         --keep-dup all
 
     cat ${ip}_peaks.${PEAK_TYPE} | wc -l | awk -v OFS='\t' '{ print "${ip}", \$1 }' | cat $peak_count_header - > ${ip}_peaks.count_mqc.tsv

nextflow.config

@@ -36,6 +36,8 @@ params {
   // Options: Peaks
   narrow_peak = false
   broad_cutoff = 0.1
+  macs_fdr = false
+  macs_pvalue = false
   min_reps_consensus = 1
   save_macs_pileup = false
   skip_peak_qc = false