Dada2: filter and trim per sample based on quality profile?

[fandemonium]

Hello,

I have been using dada2 to process my Illumina miseq amplicon sequences. I really like the tool and the flow. I do notice that in the workflow, we need to sample a few of our samples to evaluate where the quality starts to drop off by eyeballing the few plots. This is not very precise and very difficult to automate in a pipeline fashion because human interactions are required. Plus, some samples have very different drop off points (even within the read 2's).
I'm wondering if you guys have considered pull the stats part of the plotQualityProfile to a separate function and use the numbers to create trunclen arguments for each sample? Would you have any particular concerns about the approach?
Thanks! any thoughts are appreciated!
Fan

[benjjneb]

We have thought about it, but haven't put something specific together on that yet, although see this post on the filtering intuition issue: https://github.com/benjjneb/dada2/issues/236#issuecomment-422865307

We're certainly interested if anyone finds a good solution. Probably the first concern that comes to mind is ensuring that an automated method doesn't truncate paired reads too much so that they don't overlap anymore, although a solution to that might be to require the user to provide a minimum total length parameter or something similar.

[fandemonium]

Thanks for referencing the filtering intuition post. I played around with the functionplotQualityProfile by pulling out the plotdf and statdf. If I do a cut off of at the first occurrence of statdf$Mean < 25, it captures the drop really well. I did try a rolling average (window of 2) to overcome some outlier low quality base due to an air bubble or alike, but the criteria became too loose.

I had some samples did end up to be too short to merge. I looked at their error profiles and the error rates are indeed high at the truncation point. My thought is if sequences are full of errors to the point that no reliable reads can be merged, then the sample should be resequenced. Therefore, it may be more practical to set up a tryCatch statement to capture the samples that could not be merged instead of forcing the erroneous reads to merge.

[benjjneb]

My thought is if sequences are full of errors to the point that no reliable reads can be merged, then the sample should be resequenced. Therefore, it may be more practical to set up a tryCatch statement to capture the samples that could not be merged instead of forcing the erroneous reads to merge.

This would be a nice thing to include as some sort of an option in a functionalized version of what you are describing.

[fandemonium]

here is the loop I used to get pass the merge error. The below step comes after errors learned for each run (errF and errR). I passed filtered read 1 and read 2 for each sample to the loop. The loop reassigns the variable an.error.occurred based on the condition. Then based on the value of an.error.occurred, it saves to RDS or pass. There are probably better ways to do this... but I was under a time crunch and settled with this method...

an.error.occurred <- FALSE
tryCatch({
derepF <- derepFastq(A_SAMPLE_F_NAME.FQ)
ddsF <- dada(derepF, err = errF, multithread=T)
derepR <- derepFastq(gsub("_F_", "_R_", A_SAMPLE_F_NAME.FQ))
ddsR <- dada(derepR, err = errR, multithread=T)
mergePE <- mergePairs(ddsF, derepF, ddsR, derepR, verbose=T)
seqtab <- makeSequenceTable(mergePE)
seqtab.nochim <- removeBimeraDenovo(seqtab, method="consensus", multithread=TRUE, verbose=TRUE)
row.names(seqtab.nochim) <- sample.name
}, error=function(e){an.error.occurred <<- TRUE})

if (an.error.occurred == FALSE) {
saveRDS(seqtab.nochim, paste0(MAIN_PATH, "/", sample.name, "_seqtab_nochim.RDS"))
}