Hello,
I am going to use cutadapt for a 16S rRNA marine vertebrate marker (~245 bp). Since I have NovaSeq data with a different quality scoring system, I know I cannot inspect the .qzv file I generated to determine the cutoff point where read quality drops below a certain threshold. Further, since my marker is of variable length (different length amplicons produced for different species groups), is it recommended to instead set p-minimum-length in cutadapt?
Specifically, I am using 16S rRNA marine vertebrate marker (~245 bp): MarVer3F: 5′-AGACGAGAAGACCCTRTG-3′ and MarVer3R: 5′- GGATTGCGCTGTTATCCC-3’ (Valsecchi et al., 2020), and the approximate amplicon size (bp) ranges from 232-265. My current code is as follows:
qiime cutadapt trim-paired \
--i-demultiplexed-sequences [pathway] \
--p-front-f AGACGAGAAGACCCTRTG \
--p-adapter-f GGGATAACAGCGCAATCC \
--p-front-r GGATTGCGCTGTTATCCC \
--p-adapter-r CAAGGGTCTTCTCGTCT \
--p-error-rate 0.2 \
--p-overlap 5 \
--p-minimum-length 150 \
--p-discard-untrimmed \
--o-trimmed-sequences [pathway] \
--verbose
So, my question is 1) Is it recommended to set p-minimum-length in cutadapt for a marker of variable length, and 2) If so, which length should I set it to, given the smallest approximate amplicon generated (as reported in the paper) is 232 bp?
Thanks!
Hello @bioinfo_rookie,
I think it depends on whether you want partial amplicons in your dataset or not. If you're comfortable with that, then no need for minimum length filtering. If not, then you could add a minimum length filter, just note that this step is applied after primer removal (as explained here). So take that into account when calculating lengths, i.e. is the 232bp quoted length inclusive of the primer sequences?. To further complicate things, you potentially have your amplicon split across read pairs depending on the read length. I think it will be very difficult to use a length-based filter using trim-paired effectively.
Another option to enforce the no partial amplicons restriction that could be more accurate is to used linked adapters and --p-discard-untrimmed. This is explained in the docs here. Though I don't know off the top of my head if cutadapt lets you expand the linked adapter concept across read pairs. We're working on exposing the --pair-filter option (documented here), which could be useful to enforce that a pair is removed unless adapters are found in both of the read directions. This should be available next release.
Finally, remember that you can just ignore this problem right now, and remove amplicons less than your length requirement after denoising. This is probably the simplest and easiest option.
Hello again @bioinfo_rookie,
Actually, I think you might already have all the correct parts in your command. The last paragraph of this section of the docs states that the way --discard-untrimmed interacts with paired end reads is such that the pair is discarded if one or the other read direction is untrimmed. So you're basically enforcing that in a pair of reads both the forward and reverse primer must be found, which makes the length requirement unnecessary--you are now guaranteeing that each read pair covers the entire amplicon. Remember that --p-minimum-length should be removed because this is enforcing a per-read (not per-amplicon) minimum length.
Hi @colinvwood, thanks for your quick and detailed response! I'll leave out the --p-minimum-length 150 argument in this case.
Hi again @colinvwood, I have one more question about cutadapt 
For NovaSeq metabarcoding data using markers for taxonomic groups (e.g., marine vertebrates in this example, or COI universal barcode marker) is it recommended to use the default parameters for --p-error-rate and --p-overlap, or are the ones specified in the above command (0.2 and 5) better? I noticed in User guide — Cutadapt 5.2 documentation under Adapter-search parameters, the table shows the default parameters, then the "global option." It's unclear to me what the global option means, and what the advantages of using this option are. Any clarification would be much appreciated!
Hello @bioinfo_rookie,
The "global option" just means that the value will applied "globally" to all adapters, when more than one adapter is searched for. I don't think the docs are making a recommendation for --p-error-rate 0.2 and --p-overlap 5, rather these are just examples being used.
I see, thanks. Just out of curiosity, I tried comparing the output between a) --p-error-rate 0.2 and --p-overlap 5 with b) the default parameters: --p-error-rate 0.1 and --p-overlap 3. The results are as follows:
Before trimming with cutadapt:
a) --p-error-rate 0.2 and --p-overlap 5:
b) Default parameters: --p-error-rate 0.1 and --p-overlap 3
Here, I can see that the default results in slightly fewer read output, however the difference does not seem to me to be very stark. In this case, is it recommended to stick with the default, that allows for fewer errors to increase the accuracy of the output? Or, should I aim for whichever option gives the higher read output?
Hello @bioinfo_rookie,
Without a real reason to use --p-error-rate 0.2 and --p-overlap 5 I would stick with the defaults. Again, these values were arbitrarily chosen by the documentation. One is more specific than its corresponding default and the other is more sensitive than its.
