Dada2: Multithread option for collapseNoMismatch()

[masumistadler]

Hello Benjamin,

I've merged sequence tables of several dada() runs (in summary around 680 samples) and I am running the collapseNoMismatch() right now and it takes quite a while (it's running for at least a day).

Now, I have noticed that this function does not have a multithread = TRUE option as most of your computationally more expensive functions do. I was wondering whether there is a way to enable multithreading. It would be great to speed up the process.

[benjjneb]

There currently is not, but this could be an enhancement request. We have not really run into cases where this is the rate limiting step before, so haven't made it a priority, but we can consider it in the future.

[masumistadler]

That would be great for the future. Thank you!

[leholman]

I am using a set of primers that are very 'slippy' so I generate many different lengths of merged reads for the same biological sequence. collapseNoMismatch seems like just the thing to sort this out but it takes a while!

A multithreaded version would be really useful.

[masumistadler]

For your information, the collapseNoMismatch() step took around 4 days for 630 16S amplicon sequencing samples with approx. 199000 ASVs.

[benjjneb]

Thanks for the reports. I've added this issue to the 1.14 milestone, which is our next release in ~6 months.

No promises, but it's at least on the formal agenda now!

[leholman]

🥺 4 days. Thanks for the info Masumi!

I'll try and benchmark the function while I run it on the various datasets and post the results up here.

[leholman]

I have finished running the function on some datasets. See benchmarks below.

Dataset 1: marine eDNA metabarcoding of COI of ~70 samples: 37674 OTUs - 7476 seconds
Dataset 2: marine eDNA metabarcoding of 18S of ~70 samples: 10260 OTUs - 12029 seconds
Dataset 3: marine eDNA metabarcoding of 16S of ~70 samples: 19228 OTUs - 6210 seconds

Roughly 1% of sequences for the COI dataset were merged, while only 2 were merged for the 18S dataset. I pulled the code for the function and added a counter to monitor progress. The code seems to run exponentially slower as it runs through the dataset.

I also ran the RStudio code profiler on the function (http://rpubs.com/lholman/codeprofile1). Looks like grepl, substr and an if statement use most of the resources. I run into problems like this with my own code so I'm interested in different approaches to speed up these functions.

[ptpell77]

Hi,
Curious if there is now a beta version of the parallelized collapseNoMismatch command? This seems to be bottleneck in collapsing dozens of sequencing runs

[benjjneb]

Curious if there is now a beta version of the parallelized collapseNoMismatch command? This seems to be bottleneck in collapsing dozens of sequencing runs

Nope there is not. But more reports that this is a relevant bottleneck will help push this further up the queue, so thank you.

[ARBramucci]

Hi Benjamin,
I agree that this step is really slowing up the processing of multiple plates, however along with that I am curious, if processing say 60 plates wouldn't it be better for the collapse no mismatch step to come at the end of the process, so that when you combine multiple plates you are not introducing potential mismatches say if seq A and B are different by 1 nucleotide at the end of the seq, and in plate 1 A > B but in plate 2 B> A then when combined you will once again have both A and B even though they are 1 nucleotide different in length without mismatches...??

[masumistadler]

Hi,
I run the collapseNoMismatch() step always on the final sequence table that includes all plates. I do not really see the point of doing it on all the individual sequence tables separately. To my knowledge, if one would collapse ASVs on a plate basis, if you have a collapsed ASV from plate A and collapsed ASV from plate B, they will only be merged together in the merged sequence table if they have the exact same sequence. And they may be different, if you check the help page of the function, it says, the sequence that is being kept for the collapsed ASV is that of the most abundant. This means, if the sequence that is most abundant for a collapsed ASV is the same across all plates, they should be merged together when merging sequence tables, but if not, they will appear as two separate ASVs. You could run the collapse on a plate to plate basis, but I would suggest running another collapse on the final sequence table. Hope this answers your concern.

[ARBramucci]

So masumistadler, in order to do as you suggest above, you export your ASV tables as csv files and then you compile all the plates together into one file and feed it back into dada2 for the collapse no mismatch step? and you don't run into formatting issues?

[masumistadler]

No, I save all seqtab files from the individual plates as an .rds file. I read them in and join them with mergeSequenceTables(), it's a function to merge several seqtabs provided by dada2.

Here is some example code:

# read in all sequence tables from plates (there are 24)
merger <- list() # create empty list
# read in data with loop
for(i in 1:24){
merger[[i]] <- readRDS(paste0("./Objects/",i,"_seqtab.rds"))
}

# merge all tables
st.all <- mergeSequenceTables(merger[[1]],merger[[2]],merger[[3]],merger[[4]],merger[[5]],merger[[6]],merger[[7]],merger[[8]],merger[[9]],
merger[[10]],merger[[11]],merger[[12]],merger[[13]],merger[[14]],merger[[15]],merger[[16]],merger[[17]],merger[[18]],merger[[19]],
merger[[20]],merger[[21]],merger[[22]],merger[[23]],merger[[24]])
saveRDS(st.all, "./Objects/all_seqtab.rds") # save

# collapse ASVs with identical sequence
col.st <- collapseNoMismatch(st.all, minOverlap = 20, verbose = T)
saveRDS(col.st, "./Objects/collapsed_seqtab.rds")

Hope this helps!

[ARBramucci]

Thank you so much!!! I will try it now but that looks like it will work!!

[madkinsz]

Hi! Just a +1 that this is something we are noticing may be a bottleneck in our pipeline. We can report actual compute times if that would be helpful too.

[ARBramucci]

Quick question about the collapse no mismatch and running multiple plates with dada2. I have 2 plates and they have slightly different optimal truncLengths (optimal defined as what % make it all the way through the pipeline), should I optimise them both to use the same TruncLengths or does that not matter if I am doing collapse no mismatch on them after the fact because any different lengths will be collapse together?
thanks heaps!
Anna

[benjjneb]

@ARBramucci Is this paired-end data? If so, the merged reads should cover the same genetic region (i.e spanning from primer to primer) even if the truncated forward/reverse read lengths were slightly different.

If this is single-end data, I think the best option is to use the same truncation length across the runs.

[ARBramucci]

It is paired-end, 300 bp reads. So yes I agree with you that slightly different lengths should not matter. Thanks a lot for the confirmation!
Anna