Phyloseq object contains less ASV than it should

[d4straub]

Description of the bug

The following was reported via slack:

I am clustering my data and the summary_report claims I should have 52,594 clusters after vsearch clustering, but when I open the phyloseq object in R, it contains only 10,000 OTUs.

The solution was increasing --dada_assign_chunksize to 60000.

That indicates that only the a chunk of 10k (which is default for --dada_assign_chunksize) was forwarded to phyloseq. The channel input might be the problem. It could be related to clustering, maybe not. TSE might be also affected, since it has a similar (same?) input.

Command used and terminal output

Relevant files

No response

System information

nf-core/ampliseq versions 2.13.0 and 2.14.0

[tom-brekke]

The only thing I can add here is that this bug happened with a custom database and without using ITSx.

When I used the in-built UNITE database and go down the ITSx route I was still getting fewer OTUs than expected in the phyloseq object (16,902 instead of 47,343) but I think that's a different issue as it's not a clean multiple of 10k (v2.13.0, I've not tried UNITE+ITSx with 2.14.0 yet). Hopefully that will help localise where the 10k issue is occurring.

In sum:
This params file results in 10k OTUs in phyloseq in 2.13.0 and 2.14.0: nf-params_PacBioFungalOTU_Eukaryome.json

This params file results in 16,902 OTUs in phyloseq in 2.13.0: nf-params_PacBioFungalOTU_UNITE.json

Here is the custom database: Eukaryome

[d4straub]

I cannot reproduce the problem. I tried locally with

nextflow run nf-core/ampliseq -r 2.14.0 -profile test_single,singularity --outdir results -resume #chunksize is set to 250
nextflow run nf-core/ampliseq -r 2.14.0 -profile test_single,singularity --outdir results  -resume --vsearch_cluster #chunksize is set to 250
nextflow run nf-core/ampliseq -r 2.14.0 -profile test_single,singularity --outdir results  -resume --dada_assign_chunksize 100
nextflow run nf-core/ampliseq -r 2.14.0 -profile test_single,singularity --outdir results  -resume --dada_assign_chunksize 1000

number of taxa in the phyloseq object was always 456 (not clustered) or 408 (clustered) taxa. This corresponds to the numbers in the summary report.

I also tested (on cluster cfc)

nextflow run nf-core/ampliseq -r 2.14.0 -profile test_pacbio_its,cfc --outdir results -resume --sintax_ref_taxonomy false --dada_ref_taxonomy "unite-alleuk" --dada_assign_chunksize 10 --skip_qiime --skip_phyloseq false --skip_dada_taxonomy false

The numbers in the summary report was equal to the number of taxa in the phyloseq object.

I produced the number of taxa with (typically dont use phyloseq, so please correct me if I am wrong):

singularity shell work/singularity/depot.galaxyproject.org-singularity-bioconductor-phyloseq-1.46.0--r43hdfd78af_0.img
R
library(phyloseq)
a <- readRDS("dada2_phyloseq.rds")
dim(tax_table(a))

[tom-brekke]

ugh, that's frustrating. Are any of those tests with custom databases?
The error only shows up for me if I use --dada_ref_tax_custom instead of one of the packaged databases with --dada_ref_taxonomy.

I'll have another look and see if I can figure what's going on. Otherwise I'll try to make a minimal reproducible example.
Thanks for looking into it.

[JeremieCote]

Hello,
I also got this problem using ampliseq 2.13.0

I did use --dada_ref_tax_custom with this database: silva_nr99_v138.1_wSpecies_train_set.fa.gz
And I also used this: --dada_ref_tax_custom_sp silva_species_assignment_v138.1.fa.gz

No problem using --dada_ref_taxonomy.

[d4straub]

Ah thanks, I'll test with the custom options next week to find a reproducible example, That will help immensely with finding the faulty code and verify the fix ;)

[d4straub]

I can reproduce the issue with nextflow run nf-core/ampliseq -r 2.14.0 -profile test_reftaxcustom,singularity --outdir results -resume --dada_assign_chunksize 100:

• dada2_phyloseq.rds contains 100 taxa
• kraken2_phyloseq.rds contains 347 taxa
• qiime2_phyloseq.rds contains 347 taxa
• the summary report reports 347 ASVs
• DADA2_TAXONOMY is run only once instead of multiple times.

[d4straub]

Ok, I opened a PR to dev. It took considerably longer to find a proper test than to fix the code.

[d4straub]

Fix is in the dev branch, so close it here.