Hi everyone I have a cross-pipeline question about *mapping-quality filtering in nf-core pipelines. W*hile reviewing the source code, I found that different chromatin pipelines appear to use different default MAPQ thresholds: • nf-core/atacseq: MAPQ ≥1 atacseq/conf/modules.config at 2.1.2 · nf-core/atacseq • nf-core/chipseq: MAPQ ≥1 chipseq/conf/modules.config at master · nf-core/chipseq • nf-core/cutandrun: MAPQ ≥20 The filtering command uses:

samtools view -b -q ${params.minimum_alignment_q_score}

cutandrun/conf/modules.config at master · nf-core/cutandrun The default value of

minimum_alignment_q_score

is defined as 20 here: cutandrun/nextflow_schema.json at master · nf-core/cutandrun • nf-core/hic: MAPQ ≥10 hic/nextflow.config at master · nf-core/hic For nf-core/rnaseq, I could not find one pipeline-wide MAPQ cutoff. I may have missed it I understand that one universal MAPQ threshold may not be appropriate across different aligners and assays. However, ATAC-seq, ChIP-seq and CUT&RUN can all use Bowtie2, yet ATAC-seq and ChIP-seq use MAPQ ≥1 while CUT&RUN uses MAPQ ≥20. Could someone clarify: 1. Are these differences based on assay-specific benchmarking, downstream analysis requirements or historical pipeline decisions? 2. When the same aligner, such as Bowtie2, is used, what is the rationale for MAPQ ≥1 in ATAC-seq and ChIP-seq versus MAPQ ≥20 in CUT&RUN? 3. In ATAC-seq and ChIP-seq, when

keep_multi_map

is false, the pipelines apply

-q 1

. Does MAPQ ≥1 truly restrict the BAM to uniquely mapped reads? My understanding is that

-q 1

removes MAPQ 0 alignments, but does not necessarily prove that every retained read has only one possible genomic alignment. Would it therefore be more accurate to describe this as removing MAPQ 0 or highly ambiguous alignments rather than removing all multimapping reads? 4. When a pipeline has no explicit MAPQ-filtering parameter, does that mean filtering is left to the aligner or downstream tool? 5. Do the final BAM files, bigWigs, peak calls and count matrices all use the same MAPQ-filtered read set within each pipeline? 6. Would it be useful to document this consistently for every alignment-based nf-core pipeline, including: • the default MAPQ cutoff;where it is applied; • whether it is configurable; and • which downstream outputs are affected? Please correct any thresholds or interpretations above that are inaccurate. Thank you @Phil Ewels @Maxime Garcia

Hi @wani sajad that's impressive investigation. I would say these differences are mainly due to historical pipelines decisions. Which could have due to the fact that the underlying technique are different and requires different cutoff. I think @Jose Espinosa-Carrasco would have clearer answers for #CE5EL6326 and #CFP55NZ5G. @Jonathan Manning would be the best to reply for #CE8SSJV3N, and he might have some ideas about #C01FKFZ57SA. For #CJ3UVAZFT @nservant. I would say having a consistent cutoff might not be a solution for all pipelines using the same tools, but a consistent way to where it can be applied and configurable would definitively help out. I think this is actually a lot of effort to try to maintain a global cohesion on such a variety of pipelines taken care of by different teams and set of people, and documentation is often left behind. @Chris Hakkaart might have good ideas on how to deal with this. (I apologize for pinging a lot of people here)

IMO - Standardizing on a specific threshold would be extremely difficult, as pipelines are largely developed independently. Additionally, as pointed out already, each threshold may impact other tools downstream, so one rule is unlikely to fit all pipelines. However, @James Fellows Yates and I have been discussing pipeline-specific documentation to address issues like this. TLDR; we are looking at creating space in the pipeline documentation to explore nuanced technical details for pipeline configuration and customization. This is still a work in progress, so we will use this as a reference point for shaping what we are trying to do.

That's the kind of smart answer I knew you would have. Thanks Chris

Amazing job @wani sajad I agree that using the same threshold might not be the best solution but as @Chris Hakkaart mentioned above a clear documentation of why a given option is chosen in a given pipeline and how to tune it for different specific scenarios would be something very nice to have 🚀

Thank you all for the helpful responses. I agree that the most useful improvement would be pipeline-specific documentation explaining why the default threshold was selected, where it is applied, whether it is configurable, and which downstream outputs use the filtered read set. One technical point that may still be worth clarifying in the ATAC-seq and ChIP-seq documentation is the description of

-q 1

when

keep_multi_map

is false. My understanding is that this removes MAPQ 0 alignments but does not strictly demonstrate that every retained alignment is uniquely mapped. Therefore, wording such as “removes MAPQ 0 or highly ambiguous alignments” may be more precise than “retains only uniquely mapped reads,” unless uniqueness is also established through another Bowtie2 or SAM flag criterion. Thank you again. I am glad this discussion may be useful for the planned pipeline-specific technical documentation. @Maxime Garcia @Chris Hakkaart @Jose Espinosa-Carrasco

My two cents, dosen't different aligners give different MAPQ values? (e.g in contrast to Phred scores) -> Since the cut&run pipeline only uses bowtie2, I think it is easier to use a reasonable default value. I'm a bowtie2 guy and there I am pretty sure that "multimappers" qualified e.g by bowtie2 are given a MAPQ of 0. SO

keep_multi_map

false, does exactly what it says. Not sure for the other aligner but one would need to look at that. I fully agree, we should add a proper section on MAPQ usage with different aligners.... but this will become verbose 😉