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: -4
  • 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: 1
  • 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: -2
  • 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

• 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.
• Make sure your code lints (nf-core pipelines lint).
• Ensure the test suite passes (nextflow run . -profile test,docker --outdir <OUTDIR>).
• 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.
• CHANGELOG.md is updated.
• README.md is updated (including new tool citations and authors/contributors).

[d4straub]

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)

[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 : asv_consensus_, asv_mergepairs_consensus_, mergepairs_consensus_, denoising_consensus_ ?

[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.

[d4straub]

I like --mergepairs_strategy !

Ok, almost there, I still found a few points. If you dont mind those can be also added by going into the "Files changed" tab and click on "Add suggestion to batch" and then committing those.

After that additions it should work as expected and I can test again.

The last piece would be to fix the formatting here with @nf-core-bot fix linting

[nf-core-bot]

Warning

Newer version of the nf-core template is available.

Your pipeline is using an old version of the nf-core template: 3.0.2.
Please update your pipeline to the latest version.

For more documentation on how to update your pipeline, please see the nf-core documentation and Synchronisation documentation.

[d4straub]

@nf-core-bot fix linting

[d4straub]

Dear @weber8thomas , that PR seems ready to be merged now. Only missing piece to perfection is testing. You mention that you did benchmarking, is there any data that you can share? The test data here would be only for testing technically if everything is fine. Do you think any paired-end data would be sufficient? If yes, then I would simply use --mergepairs_strategy consensus in any already existing test run.

[weber8thomas]

Dear @weber8thomas , that PR seems ready to be merged now. Only missing piece to perfection is testing. You mention that you did benchmarking, is there any data that you can share? The test data here would be only for testing technically if everything is fine. Do you think any paired-end data would be sufficient? If yes, then I would simply use --mergepairs_strategy consensus in any already existing test run.

Hi @d4straub , sorry for the late reply (and happy new year 🥳 ) We generated additional benchmarks in order to validate and justify custom values we are using for Dada2 alignemnt (gap, mismatch, match), especially using a more robust grid search. I'll share this soon.
Regarding the test dataset, any PE data should be sufficient to technically validate the function is working fine.
Just out of curiosity, would there be any light datasets containing both eukaryotes and prokaryotes we could leverage in ampliseq test datasets repo?
Thanks!

[d4straub]

I'll share this soon.

Great, looking forward!

Regarding the test dataset, any PE data should be sufficient to technically validate the function is working fine.

Ok, I'll do that then in a separate PR.

Just out of curiosity, would there be any light datasets containing both eukaryotes and prokaryotes we could leverage in ampliseq test datasets repo?

I am actually not sure, several people added test datasets and I never looked at them in that way.

[d4straub]

Regarding the test dataset, any PE data should be sufficient to technically validate the function is working fine.

Alright, I'm going to activate --mergepairs_strategy consensus in at least one test profile in a separate PR

Just out of curiosity, would there be any light datasets containing both eukaryotes and prokaryotes we could leverage in ampliseq test datasets repo?

I am not sure, I dont think so.

I would prefer to merge that PR now so that I can work towards a release (whatever little time I have for that). Any objections?

[weber8thomas]

Regarding the test dataset, any PE data should be sufficient to technically validate the function is working fine.

Alright, I'm going to activate --mergepairs_strategy consensus in at least one test profile in a separate PR

Just out of curiosity, would there be any light datasets containing both eukaryotes and prokaryotes we could leverage in ampliseq test datasets repo?

I am not sure, I dont think so.

I would prefer to merge that PR now so that I can work towards a release (whatever little time I have for that). Any objections?

Hi @d4straub , I'm just waiting for validation from @nhenry50 and @FloraVincent, will update you asap!

[weber8thomas]

I just updated the values accordingly to the best performance in our benchmark, Flora & Nicolas validated the merge. You're good to go! :)

[d4straub]

Suggested change

	saveRDS(mergers, "${meta.run}.mergers.rds")
	saveRDS(mergers, "${prefix}.mergers.rds")

better late than never, this seems to me as a breaking change under certain circumstances, was this intentional?