Adding the feature of shifting and splitting the reads

Hi,

I was wondering if adding the feature of shifting and splitting the reads after the alignment would be helpful.

A recent review on the approaches to analyze ATAC-seq data (Yan et al., 2020) suggests reads should be shifted + 4 bp and − 5 bp for positive and negative strand respectively, to account for the 9-bp duplication created by DNA repair of the nick by Tn5 transposase and achieve base-pair resolution of TF footprint and motif related analyses.

Also, in many studies researchers have looked into the nucleosomal-free regions and nucleosome associated regions separately. For example Yoshida et al., 2019 .

Here's what I am proposing.

Let's have three modes of post-alignment analysis for the ATAC-seq pipeline.

1. No shifting and splitting the aligned reads
2. Shifting the reads and splitting them into two - NFR (Nucleosomal Free Regions) and NBR (Nucleosomal Bound regions)
3. Shifting the reads and splitting them into four - NFR (Nucleosomal Free Regions), Mono-, di- and tri-nucleosomal bound regions.

The first approach (which is the current atacseq pipeline) is useful in almost all cases where identifying the open chromatin regions is the objective. The latter would be helpful if the exact cut sites of transposase is important like motif analyses, footprinting etc., or when the analysis requires looking into the NFR/NBR regions separately.

The downstream analysis and QC could be done depending on which of these options the user chooses. For example, if the user chooses option 2, we could generate mergedLibrary and mergedReplicate BAM, bigWigs etc separately for NFR and NBR and perform peak calling and all the QC for them separately.

I have used deeptools in the past to perform shifting and splitting. Here's an example code for the same.

bam_loc="/path/to/aligned/bamfiles/directory"
out_loc="/path/to/output/directory"

## For NFR (Frag length of 120 bp and less)

alignmentSieve --bam $bam_loc/${SAMPLE}.mLb.clN.sorted.bam \
               --outFile $out_loc/${SAMPLE}.mLb.ss.NFR.bam \
               --ATACshift \
               -p 8 \
               --smartLabels \
               --minFragmentLength 0 \
               --maxFragmentLength 120 \
               --filterMetrics $out_loc/${SAMPLE}.mLb.ss.NFR.log

## For NBR (Frag length of 180 bp and above)

alignmentSieve --bam $bam_loc/${SAMPLE}.mLb.clN.sorted.bam \
               --outFile $out_loc/${SAMPLE}.mLb.ss.NBR.bam \
               --ATACshift \
               -p 8 \
               --smartLabels \
               --minFragmentLength 180 \
               --filterMetrics $out_loc/${SAMPLE}.mLb.ss.NBR.log

It would be great to know comments and suggestions from the nf-core members. I would be happy to discuss. I am a beginner in the nf-core community but would be happy to contribute to this pipeline enhancement.

Review the current atacseq pipeline and the proposed alignmentSieve commands first. Define how the three analysis modes should select shifted reads and fragment-size groups, then identify which merged BAMs, bigWigs, peak calls, and QC outputs must be produced for each selection. Done means the chosen mode consistently produces the corresponding split outputs and downstream analyses.