make duplicate handling configurable

#47

CITATION.cff

@@ -24,6 +24,9 @@ authors:
     family-names: Folkman
     orcid: 'https://orcid.org/0000-0002-5811-8875'
     affiliation: CeMM Research Center for Molecular Medicine
+  - given-names: Fangwen
+    family-names: Zhao
+    affiliation: CeMM Research Center for Molecular Medicine
   - given-names: Lina
     family-names: Dobnikar
     affiliation: CeMM Research Center for Molecular Medicine

README.md

@@ -24,6 +24,7 @@ A [Snakemake 8](https://snakemake.readthedocs.io/en/stable/) workflow implementa
 - [Bekir Ergüner](https://github.com/berguner)
 - [Daniele Barreca](https://github.com/DanieleBarreca)
 - [Lukas Folkman](https://github.com/lukas-folkman)
+- [Fangwen Zhao](https://github.com/fwzhao)
 - [Lina Dobnikar](https://github.com/ld401)
 - [Christoph Bock](https://github.com/chrbock)
 
@@ -51,7 +52,7 @@ This project wouldn't be possible without the following software and their depen
 This is a template for the Methods section of a scientific publication and is intended to serve as a starting point. Only retain paragraphs relevant to your analysis. References [ref] to the respective publications are curated in the software table above. Versions (ver) have to be read out from the respective conda environment specifications (`workflow/envs/*.yaml file`) or post-execution in the result directory (`atacseq_pipeline/envs/*.yaml`). Parameters that have to be adapted depending on the data or workflow configurations are denoted in squared brackets e.g., [X].
 
 **Processing.**
-Sequencing adapters were removed using the software fastp (ver) [ref]. Bowtie2 (ver) [ref] was used for the alignment of the short reads (representing locations of transposition events) to the [GRCh38 (hg38)/GRCm38 (mm10)] assembly of the [human/mouse] genome using the “--very-sensitive” parameter. PCR duplicates were marked using samblaster (ver) [ref]. Aligned BAM files were then sorted, filtered using ENCODE blacklisted regions [ref], and indexed using samtools (ver) [ref]. To detect the open chromatin regions, peak calling was performed using MACS2 (ver) [ref] using the “--nomodel”, “--keep-dup auto” and “--extsize 147” options on each sample. HOMER (ver) [ref] function findMotifs was used for motif enrichment analysis of the detected open chromatin regions. Quality control metrics were aggregated and reported using MultiQC (ver) [ref], [X] sample(s) needed to be removed.
+Sequencing adapters were removed using the software fastp (ver) [ref]. Bowtie2 (ver) [ref] was used for the alignment of the short reads (representing locations of transposition events) to the [GRCh38 (hg38)/GRCm38 (mm10)] assembly of the [human/mouse] genome using the “--very-sensitive” parameter. PCR duplicates were marked using samblaster (ver) [ref]. Aligned BAM files were then sorted, filtered using ENCODE blacklisted regions [ref], samtools view flags [SAM_flag], and indexed using samtools (ver) [ref]. To detect the open chromatin regions, peak calling was performed using MACS2 (ver) [ref] using the “--nomodel”, “--keep-dup [macs2_keep_dup]” and “--extsize 147” options on each sample. HOMER (ver) [ref] function findMotifs was used for motif enrichment analysis of the detected open chromatin regions. Quality control metrics were aggregated and reported using MultiQC (ver) [ref], [X] sample(s) needed to be removed.
 
 **Quantification.**
 A consensus region set, comprising of [X] genomic regions, was generated, by merging the identified peak summits, extended by [slop_extension]bp on both sides using the slop function from bedtools (ver) [ref] and pybedtools (ver) [ref], across all samples while again discarding peaks overlapping blacklisted features as defined by the ENCODE project [ref]. The consensus region set was used to quantify the chromatin accessibility in each sample by summing the number of reads overlapping each consensus region. The consensus region set, and sample-wise quantification of accessibility was performed using bedtools (ver) [ref] and pybedtools (ver) [ref]. Furthermore, the consensus region set was used to quantify the peak support per sample and each region was mapped to its closest TSS according to the HOMER annotation within proximal TSS up and downstream distances [proximal_size_up/down] yielding a gene/TSS-based quantification. Complementary, all promoter regions, defined by the same parameters, were quantified for each sample and aggregated to yield a gene/promoter-based quantification. Finally, all sample-wise enriched known motifs according to HOMER were aggregated.
@@ -62,30 +63,43 @@ Consensus regions were annotated using annotatePeaks function from HOMER (ver) [
 The processing and quantification described here was performed using a publicly available Snakemake [ver] (ref) workflow [[10.5281/zenodo.6323634](https://doi.org/10.5281/zenodo.6323634)].
 
 # 🚀 Features
-- Processing (results/)
+- Processing (`results/`)
     - Alignment of both single-end and paired-end reads in raw/unaligned BAM format with Bowtie2.
-    - Peak calling with MAC2.
-      - Note: even though the peak support of a region in a certain sample is 0, does not mean that there are no reads counted in the count matrix, it just states that there was no peak called.
-      - Note: the peak support can be >1 for certain samples in case of a consensus region spanning more than one peak within the respective sample.
+      - Filtering using `samtools view` can be configured using [SAM Flags](https://broadinstitute.github.io/picard/explain-flags.html) (`SAM_flag`).
+    - Peak calling with `MACS2`.
+      - Duplicate handling can be configured using `macs2_keep_dup`.
+      - Even though the peak support of a region in a certain sample is 0, does not mean that there are no reads counted in the count matrix, it just states that there was no peak called.
+      - The peak support can be >1 for certain samples in case of a consensus region spanning more than one peak within the respective sample.
     - Peak annotation and motif analysis HOMER.
-    - Quantification of TSS coverage (results/).
-- Reporting (report/)
+    - Quantification of TSS coverage.
+- Reporting (`report/`)
     - MultiQC report generation using MultiQC, extended with an in-house developed plugin [atacseq_report](./workflow/scripts/multiqc_atacseq).
-- Quantification (counts/)
-    - Consensus region set generation across all called peaks (consensus_regions.bed).
-    - Read count quantification of the consensus regions across samples, yielding a count matrix with dimensions consensus regions X samples (consensus_regions.bed, consensus_counts.csv).
-    - Peak support quantification of the consensus regions across samples, yielding a count matrix with dimensions consensus regions X samples (support_counts.csv).
-    - Consensus regions mapped to closest gene TSS according to HOMER (Distance to TSS) within proximal TSS up and downstream distances (TSS_regions.bed, TSS_counts.csv, TSS_annotation.csv).
-    - Read count quantification of promoter regions based on provided proximal TSS up and downstream distances (promoter_regions.bed, promoter_counts.csv, promoter_annotation.csv).
-      - [Pseudoautosomal regions in human](https://www.ensembl.org/info/genome/genebuild/human_PARS.html) chromosome Y are skipped.
-    - Aggregation of all sample-wise HOMER known motif enrichment results into one CSV in long-format (HOMER_knownMotifs.csv).
-- Annotation (counts/)
+- Quantification (`counts/`)
+    - Consensus region set generation across all called peaks (`consensus_regions.bed`).
+    - Read count quantification of the consensus regions across samples, yielding a count matrix with dimensions consensus regions X samples (`consensus_counts.csv`).
+    - Peak support quantification of the consensus regions across samples, yielding a count matrix with dimensions consensus regions X samples (`support_counts.csv`).
+    - Consensus regions mapped to closest gene TSS according to HOMER (Distance to TSS) within proximal TSS up and downstream distances (`TSS_regions.bed`, `TSS_counts.csv`, `TSS_annotation.csv`).
+    - Read count quantification of promoter regions based on provided proximal TSS up and downstream distances (`promoter_regions.bed`, `promoter_counts.csv`, `promoter_annotation.csv`).
+      - [Pseudoautosomal regions in human](https://www.ensembl.org/info/genome/genebuild/human_PARS.html) chromosome `Y` are skipped.
+    - Aggregation of all sample-wise HOMER known motif enrichment results into one CSV in long-format (`HOMER_knownMotifs.csv`).
+- Annotation (`counts/`)
     - Sample annotation file based on MultiQC general stats and provided annotations for downstream analysis (`sample_annotation.csv`).
     - Consensus region set annotation using (`consensus_annotation.csv`)
       - UROPA with regulatory build and gencode as references, configurable here: `workflow/resources/UROPA/*.txt`.
       - HOMER with `annotatePeaks.pl`
       - bedtools for nucleotide counts/content (e.g., % of GC)
 
+> [!IMPORTANT]  
+> **Duplciate reads** can be filtered during the alignment step by `samtools` and/or ignored during peak calling by `MACS2`.
+> 
+> The inclusion of duplicates should be intentional, and may lead to a large number of consensus regions.
+> 
+> The removal of duplicates should be intentional, might remove real biological signal.
+> 
+> The decision depends on your downstream analysis steps e.g., rigorous filtering (e.g., using `edgeR::filterByExpr`) and/or accounting for PCR bias by normalization conditional on genomic region length and GC content (e.g., [CQN](https://academic.oup.com/biostatistics/article/13/2/204/1746212)) and goals (e.g., differential accessibility analysis).
+> 
+> We recommend reading this ChIP-seq tutorial's section on ["Removing redundancy"](https://hbctraining.github.io/Intro-to-ChIPseq/lessons/05_peak_calling_macs.html).
+
 # 🛠️ Usage
 These steps are the recommended usage for this workflow:
 
@@ -170,6 +184,7 @@ Finally, a previous PhD student in our lab, [André Rendeiro](https://orcid.org/
   - [Unsupervised Analysis](https://github.com/epigen/unsupervised_analysis) to understand and visualize similarities and variations between cells/samples, including dimensionality reduction and cluster analysis. Useful for all tabular data including single-cell and bulk sequencing data.
   - [Differential Analysis with limma](https://github.com/epigen/dea_limma) to identify and visualize statistically significantly different features (e.g., genes or genomic regions) between sample groups.
   - [Enrichment Analysis](https://github.com/epigen/enrichment_analysis) for biomedical interpretation of (differential) analysis results using prior knowledge.
+- [Introduction to ChIP-seq using high performance computing](https://hbctraining.github.io/Intro-to-ChIPseq/)
 
 # 📑 Publications
 The following publications successfully used this module for their analyses.

config/config.yaml

@@ -23,19 +23,19 @@ sequencing_center: CeMM_BSF
 annot_columns: ['pass_qc','read_type','organism'] # (optional) can be empty [""]. must be columns in the annotation sheet
 
 ##### QUANTIFICATION #####
+
+# samtools view filtering SAM flags
+# https://broadinstitute.github.io/picard/explain-flags.html
+SAM_flag: 2316  #  int; keep duplictaes 2316 or filter duplicates 3340
+
 # coverage calculation with bedtools for QC report
 tss_slop: 2000
 noise_lower: 100
 
-# specify if duplicates should be kept during filtering bam files to define samtools view filtering flags
-# filtered reads used as input for macs2 peak calling and counts
-# warning: inclusion of duplicates should be intentional, and may lead to a large number of consensus regions
-remove_dup: True  #  bool: True by default
-
-# specify how duplicate reads should be handled by macs2
-# warning: inclusion of duplicates should be intentional, and may lead to a large number of consensus regions
-# see documentation for parameter: https://manpages.ubuntu.com/manpages/mantic/man1/macs2_callpeak.1.html
-keep_dup: 1  # int: default = 1, int for kept reads at given genomic coordinates. or "auto"
+# MACS2 --keep-dup paramater for duplicate read handling during peak calling
+# https://manpages.ubuntu.com/manpages/mantic/man1/macs2_callpeak.1.html
+# https://hbctraining.github.io/Intro-to-ChIPseq/lessons/05_peak_calling_macs.html
+macs2_keep_dup: "auto"  # int: default = 1, int for kept reads at given genomic coordinates; or "auto"
 
 # determination of consensus regions using (py)bedtools
 slop_extension: 250

test/hg38test/config/hg38test_atacseq_pipeline_config.yaml

@@ -2,7 +2,6 @@
 ##### RESOURCES #####
 mem: '32000'
 threads: 2
-partition: 'shortq'
 
 ##### GENERAL #####
 project_name: hg38test #MyATACproject # name of the project/dataset
@@ -24,10 +23,20 @@ sequencing_center: CeMM_BSF
 annot_columns: ['pass_qc','read_type','organism'] # (optional) can be empty [""]. must be columns in the annotation sheet
 
 ##### QUANTIFICATION #####
+
+# samtools view filtering SAM flags
+# https://broadinstitute.github.io/picard/explain-flags.html
+SAM_flag: 2316  #  int; keep duplictaes 2316 or filter duplicates 3340
+
 # coverage calculation with bedtools for QC report
 tss_slop: 2000
 noise_lower: 100
 
+# MACS2 --keep-dup paramater for duplicate read handling during peak calling
+# https://manpages.ubuntu.com/manpages/mantic/man1/macs2_callpeak.1.html
+# https://hbctraining.github.io/Intro-to-ChIPseq/lessons/05_peak_calling_macs.html
+macs2_keep_dup: "auto"  # int: default = 1, int for kept reads at given genomic coordinates; or "auto"
+
 # determination of consensus regions using (py)bedtools
 slop_extension: 250
 

test/mm10test/config/mm10test_atacseq_pipeline_config.yaml

@@ -2,7 +2,6 @@
 ##### RESOURCES #####
 mem: '32000'
 threads: 2
-partition: 'shortq'
 
 ##### GENERAL #####
 project_name: mm10test #MyATACproject # name of the project/dataset
@@ -24,10 +23,20 @@ sequencing_center: CeMM_BSF
 annot_columns: ['pass_qc','read_type','organism'] # (optional) can be empty [""]. must be columns in the annotation sheet
 
 ##### QUANTIFICATION #####
+
+# samtools view filtering SAM flags
+# https://broadinstitute.github.io/picard/explain-flags.html
+SAM_flag: 2316  #  int; keep duplictaes 2316 or filter duplicates 3340
+
 # coverage calculation with bedtools for QC report
 tss_slop: 2000
 noise_lower: 100
 
+# MACS2 --keep-dup paramater for duplicate read handling during peak calling
+# https://manpages.ubuntu.com/manpages/mantic/man1/macs2_callpeak.1.html
+# https://hbctraining.github.io/Intro-to-ChIPseq/lessons/05_peak_calling_macs.html
+macs2_keep_dup: "auto"  # int: default = 1, int for kept reads at given genomic coordinates; or "auto"
+
 # determination of consensus regions using (py)bedtools
 slop_extension: 250
 

workflow/rules/processing.smk

@@ -17,14 +17,12 @@ rule align:
         samtools_flagstat_log = os.path.join(result_path, 'results', "{sample}", 'mapped', '{sample}.samtools_flagstat.log'),
         stats = os.path.join(result_path, 'results', "{sample}", '{sample}.align.stats.tsv'),
     params:
-        # alignment parameters
         interleaved_in = lambda w: "--interleaved_in" if samples["{}".format(w.sample)]["read_type"] == "paired"  else " ",
         interleaved = lambda w: "--interleaved" if samples["{}".format(w.sample)]["read_type"] == "paired" else " ",
-        filtering = lambda w: "-q 30 -F {flag} -f 2 -L {whitelist}".format(flag=3340 if config['remove_dup'] else 2316, whitelist=config["whitelisted_regions"]) if samples["{}".format(w.sample)]["read_type"] == "paired" else "-q 30 -F {flag} -L {whitelist}".format(flag=3340 if config['remove_dup'] else 2316, whitelist=config["whitelisted_regions"]),
+        filtering = lambda w: "-q 30 -F {flag} -f 2 -L {whitelist}".format(flag=config['SAM_flag'], whitelist=config["whitelisted_regions"]) if samples["{}".format(w.sample)]["read_type"] == "paired" else "-q 30 -F {flag} -L {whitelist}".format(flag=config['SAM_flag'], whitelist=config["whitelisted_regions"]),
         add_mate_tags = lambda w: "--addMateTags" if samples["{}".format(w.sample)]["read_type"] == "paired" else " ",
         adapter_sequence = "-a " + config["adapter_sequence"] if config["adapter_sequence"] !="" else " ",
         adapter_fasta = "--adapter_fasta " + config["adapter_fasta"] if config["adapter_fasta"] !="" else " ",
-        # configs
         sequencing_platform = config["sequencing_platform"],
         sequencing_center = config["sequencing_center"],
         mitochondria_name = config["mitochondria_name"],
@@ -64,10 +62,8 @@ rule tss_coverage:
     output:
         tss_hist = os.path.join(result_path,"results","{sample}","{sample}.tss_histogram.csv"),
     params:
-        # parameters for coverage
         noise_upper = ( config["tss_slop"] * 2 ) - config["noise_lower"],
         double_slop = ( config["tss_slop"] * 2 ),
-        # configs
         genome_size = config["genome_size"],
         tss_slop = config["tss_slop"],
         unique_tss = config["unique_tss"],
@@ -106,17 +102,14 @@ rule peak_calling:
         macs2_log = os.path.join(result_path, 'results', "{sample}", 'peaks', '{sample}.macs2.log'),
         stats = os.path.join(result_path, 'results', "{sample}", '{sample}.peak.stats.tsv'),
     params:
-        # paths
         peaks_dir = os.path.join(result_path,"results","{sample}","peaks"),
         homer_dir = os.path.join(result_path,"results","{sample}","homer"),
         homer_bin = os.path.join(HOMER_path,"bin"),
-        # peak calling parameters
         formating = lambda w: '--format BAMPE' if samples["{}".format(w.sample)]["read_type"] == "paired" else '--format BAM',
-        # configs
         genome_size = config["genome_size"],
         genome = config["genome"],
         regulatory_regions = config["regulatory_regions"],
-        keep_dup = config['keep_dup'],
+        keep_dup = config['macs2_keep_dup'],
     resources:
         mem_mb=config.get("mem", "16000"),
     threads: 4*config.get("threads", 2)