Comparison of observed sequences to expected

[d4straub]

I think it would be a good idea to have a comparison of observed versus expected results for sequences, abundances, and taxonomic assignments for mock samples and the like.

This PR adds a first step, comparing observed ASV sequences against expected reference sequence(s) and report sequence mismatches/match rates.

Included steps:

• compare generated ASV sequences to --expected_sequences (fasta file) with VSEARCH --usearchglobal
• interpret results for each sample with tables and a barplot (y = number of sequences, x = number of mismatches)
• also supplying an abundance table (tsv) with --expected_abundances will filter for each sample the expected sequences (without it, it is assumed that all samples have all sequences in --expected_sequences)

Example barplot:

Also added a CI test, data via nf-core/test-datasets#2117

Further PRs might add in tentative priority
(1) FN/FP/TP, F1 score of (perfect) matching sequences and
(2) evaluate how abundances compare to expected
(3) compare taxonomic classification to expected

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

[erikrikarddaniel]

Looks fine, and I see the point.

It strikes me, however, that this might be better implemented as a separate pipeline that reads ampliseq output files. If a user always uses --outdir results the parameters could even be reused between calls, except maybe for the path to the mock/expected data. I believe this would be easier in terms of maintenance, particularly for new features of this type. (I can imagine building some sort of test on ecological distances between observed and expected mock structure.) Maybe this could one day be turned into a more general tool also applicable for e.g. metagenome and metatranscriptome mock data.

With metapipelines coming one day, calling the two pipelines would be easier too.

Happy to have in Ampliseq, and happy to discuss more.

[d4straub]

Thanks @erikrikarddaniel for the review and discussion! I found a few more points, I would be happy to hear your thoughts on them, if you have the time.

It strikes me, however, that this might be better implemented as a separate pipeline that reads ampliseq output files. If a user always uses --outdir results the parameters could even be reused between calls, except maybe for the path to the mock/expected data. I believe this would be easier in terms of maintenance, particularly for new features of this type.

Yes, that would be possibly a good idea.

[d4straub]

This is dangerous in case column names overlap. This can be if at least one sample name is identical to a column name of the VSEARCH output. The likelihood can be reduced by minimizing VSEARCH output columns and adding a test whether column names and samples overlap. Do you agree?

[erikrikarddaniel]

When samples are column names this is of course a risk one should look for. Personally, I only work with "long" tables, i.e. tables with, typically, three or four columns: sample, asv_id, count (and relab if I've calculated relative abundance). This avoids the problem you see here, and makes data smaller since one can remove zeroes.

[d4straub]

Thats a great point. I didnt consider that. I'll think about it.

[d4straub]

While a long format would have been a great solution, it would need some adaptions. So I omitted merging tables that could have identical column names. This way it works as well, no conflicts possible whatsoever.

[d4straub]

Thanks a lot for your advice!!