Dada2 MergePairs : consensus between merging & concatenating reads

[weber8thomas]

Authors

• @weber8thomas - Thomas Weber (EMBL Heidelberg)- Technical implementation, method parameter configuration accessible from the pipeline config, testing, and ensuring compliance with nf-core requirements.
• @nhenry50 - Nicolas Henry (CNRS - ABiMS) - Concept & Methods development, algorithm design, benchmarking
• @lplanat - Laurine Planat (EMBL Heidelberg)- Testing, feedback provision, benchmarking
• @FloraVincent - Flora Vincent (EMBL Heidelberg) - Conceptualisation, project supervision

Context & description

This pull request introduces an enhancement to the DADA2::mergePairs() process within the pipeline, enabling conditional merge or concatenation of sequences based on the overlap between forward and reverse reads. Previously, the pipeline allowed only to use either one of the two methods (merging or concatenating) of paired-end reads, via the --concatenate_reads parameter. This enhancement introduces the "consensus" method, allowing the pipeline to dynamically determine the appropriate method—merging or concatenation—thereby improving sequence assembly accuracy and downstream analysis outcomes.

The core enhancement revolves around incorporating conditional logic to assess the overlap between paired-end reads and decide whether to merge or concatenate them. This decision is based on a specified overlap threshold, ensuring that only reads with adequate overlap are merged, while others are concatenated with a defined spacer.

Enhancement Highlights:

• Dual Invocation of mergePairs:
  • Merging: Invoked with justConcatenate = FALSE to attempt merging where possible.
  • Concatenation: Invoked with justConcatenate = TRUE to concatenate reads where merging isn't feasible.
• Overlap Threshold Calculation:
  • Calculates a minimum overlap threshold (min_overlap_obs) based on accepted mergers.
  • (By default) Utilizes the 0.1th percentile (quantile(min_overlap_obs, 0.001)) to determine a stringent cutoff (can be customised through asv_percentile_cutoff). This ensures that only read pairs with exceptionally high overlap are merged into consensus sequences, while those with insufficient overlap are concatenated, thereby maintaining sequence accuracy and integrity.
• Conditional Replacement:
  • Iterates through each sample's mergers.
  • Replaces non-accepted mergers with concatenated sequences if the overlap falls below the threshold.
  • Filters out any non-accepted, non-concatenated sequences to maintain data integrity.

Parameters changed or introduced

• mergepairs_strategy (former concatenate_reads ; expect now a string as value)
  • Options:
  • "merge" (former TRUE): Enable concatenation (already existing)
  • "concatenate" (former FALSE): Disable concatenation (already existing)
  • "consensus": Enables conditional merging or concatenation based on the overlap between reads.

By default, when setting mergepairs_strategy to merge or concatenate, fixed generic values will still be used regarding the paired-reads alignment (see below):

match = 1, mismatch = -64, gap = -64, minOverlap = 12, maxMismatch = 0

However, by setting mergepairs_strategy to consensus, we offer now the possibility to tweak those values in order to adjust the alignment done between the reads of the same pair. This can be done my modifying the following newly introduced parameters:

• mergepairs_consensus_minoverlap
  • Description: Sets the minimum required overlap length for merging paired-end reads.
  • Default Value: 12
  • Usage: Determines the threshold below which reads will not be merged and may be concatenated instead.
• mergepairs_consensus_maxmismatch
  • Description: Defines the maximum allowed mismatches during the merging process.
  • Default Value: 0
  • Usage: Controls the stringency of the merging criteria by limiting the number of mismatches permitted between overlapping regions.
• mergepairs_consensus_gap
  • Description: Specifies the gap penalty used during the alignment process in merging.
  • Default Value: -64
  • Usage: Influences the alignment algorithm's handling of gaps, affecting the quality of merged sequences.
• mergepairs_consensus_match
  • Description: Sets the match score for the alignment algorithm.
  • Default Value: 5
  • Usage: Determines the scoring for matching bases during the alignment, impacting the alignment sensitivity.
• mergepairs_consensus_mismatch
  • Description: Sets the mismatch penalty for the alignment algorithm.
  • Default Value: -6
  • Usage: Determines the penalty for mismatched bases during alignment, affecting the alignment specificity.
• mergepairs_consensus_percentile_cutoff
  • Description: Sets the percentile cutoff determining the minimum observed overlap in the dataset.
  • Default Value: 0.001
  • Usage: This ensures that only read pairs with high overlap are merged into consensus sequences. Those with insufficient overlap are concatenated.

This approach is particularly beneficial for datasets containing both prokaryotic and eukaryotic sequences, where lengths vary, leading to differing overlap extents.

PR checklist

• [x] This comment contains a description of changes (with reason).
• [ ] If you've fixed a bug or added code that should be tested, add tests!
• [ ] If you've added a new tool - have you followed the pipeline conventions in the contribution docs
• [ ] If necessary, also make a PR on the nf-core/ampliseq branch on the nf-core/test-datasets repository.
• [x] Make sure your code lints (nf-core pipelines lint).
• [x] Ensure the test suite passes (nextflow run . -profile test,docker --outdir <OUTDIR>).
• [x] Check for unexpected warnings in debug mode (nextflow run . -profile debug,test,docker --outdir <OUTDIR>).
• [ ] Usage Documentation in docs/usage.md is updated.
• [ ] Output Documentation in docs/output.md is updated.
• [x] CHANGELOG.md is updated.
• [ ] README.md is updated (including new tool citations and authors/contributors).

[weber8thomas]

Thanks, that looks much cleaner!

(1) When I run

nextflow run weber8thomas/ampliseq -r consensus-pr -profile test,singularity --outdir results_consensus-pr -resume --skip_qiime

the warning

WARN: Access to undefined parameter `concatenate_reads` -- Initialise it to a default value eg. `params.concatenate_reads = some_value

(2) when appending --asv_concatenate_reads I get the warning

ERROR ~ Validation of pipeline parameters failed!

 -- Check '.nextflow.log' file for details
The following invalid input values have been detected:

* --asv_concatenate_reads (true): Expected any of [[false, true, consensus]]

(3) appending --concatenate_reads "consensus" results into

ERROR ~ Unknown config attribute `params.match` -- check config file: /home/daniel/.nextflow/assets/weber8thomas/ampliseq/nextflow.config

I think all problems should be solved by my comments. Edit: thats not true, all occurrences of concatenate_reads have to be replaced by asv_concatenate_reads to solve (1)

Thanks for the feedback @d4straub ! I implemented your suggestions and the execution worked on my side with and without --asv_concatenate_reads consensus

EDIT : just noticed I have an issue with how match & mismatch values should be defined, here is the original code from @nhenry50 below:

        ext.args2 = [
            'minOverlap = 12, maxMismatch = 0, propagateCol = character(0), gap = -64, homo_gap = NULL, endsfree = TRUE, vec = FALSE',
            params.concatenate_reads == "consensus" ? "returnRejects = TRUE, match = 5, mismatch = -6" :
                params.concatenate_reads == "concatenate" ? "justConcatenate = TRUE, returnRejects = FALSE, match = 1, mismatch = -64" :
                "justConcatenate = FALSE, returnRejects = FALSE, match = 1, mismatch = -64"
        ].join(',').replaceAll('(,)*$', "")

How could I define default values to match = 1, mismatch = -64 if consensus is not used, and match = 5, mismatch = -6 if used ?

[d4straub]

How could I define default values to match = 1, mismatch = -64 if consensus is not used, and match = 5, mismatch = -6 if used ?

That might be an issue. I can currently only think of using two params (introducing one more), which will complicate things, or remove the params and fix the values in the config. That way the values are still modifiable using a config, but less accessible for the standard user because working with configs requires more knowledge. Therefore, if those values do not need to be changed all the time and the values you propose here are usually fine, I'd say remove the params and put the values in the config (e.g. mismatch = -64 instead of mismatch = ${params.asv_mismatch} and remove --asv_mismatch altogether).

[weber8thomas]

How could I define default values to match = 1, mismatch = -64 if consensus is not used, and match = 5, mismatch = -6 if used ?

That might be an issue. I can currently only think of using two params (introducing one more), which will complicate things, or remove the params and fix the values in the config. That way the values are still modifiable using a config, but less accessible for the standard user because working with configs requires more knowledge. Therefore, if those values do not need to be changed all the time and the values you propose here are usually fine, I'd say remove the params and put the values in the config (e.g. mismatch = -64 instead of mismatch = ${params.asv_mismatch}).

Actually, match = 1, mismatch = -64 should be fixed when consensus is not used but we would like to give the possibility to adjust values when consensus is enabled. Would the last commit be okay in that regard? The user would only be able to play with parameters like asv_match asv_mismatch ... when --asv_concatenate_reads consensus is set.

[d4straub]

Would the last commit be okay in that regard?

Yes, but the documentation should be adjusted accordingly. And the parameter names are also not really fitting any more.

[weber8thomas]

Would the last commit be okay in that regard?

Yes, but the documentation should be adjusted accordingly. And the parameter names are also not really fitting any more.

Will work on this! Which prefix would you suggest regarding the parameter names : asvconsensus, asv_mergepairsconsensus, mergepairsconsensus, denoisingconsensus ?

[d4straub]

mergepairs_consensus_ sounds good I think. Thats precise. Could you also rename the asv_ prefix? mergepairs_ seems better suited, since it only affects PE reads!? And asv_concatenate_reads seems now that it has 3 choices also imperfect. Maybe just mergepairs, or mergepairs_type or such? Not sure, I'm sure you'll find a intuitive name.

[weber8thomas]

mergepairs_consensus_ sounds good I think. Thats precise. Could you also rename the asv_ prefix? mergepairs_ seems better suited, since it only affects PE reads!? And asv_concatenate_reads seems now that it has 3 choices also imperfect. Maybe just mergepairs, or mergepairs_type or such? Not sure, I'm sure you'll find a intuitive name.

Thanks for the feedback! :)

I updated all the different parameters introduced with prefix mergepairs_consensus_ and changed asv_concatenate_reads to mergepairs_strategy, I think it capture properly the idea of that one.