cutadapt doesn't trim all primers
0
0
Entering edit mode
3.2 years ago

Hi together, I'm new to the analysis of amplicon sequencing reads with dada2 and cutadapt in R. I'm following this tutorial (worked with the test data set) integrating my own data: https://benjjneb.github.io/dada2/ITS_workflow.html After the step, in which cutadapt trims the primers from my sequences, it doesn't detect the forward primer in all sequences. I assume that I would have to change some settings so that the leftover primer will also be found. Here is what I do:

# Loading cutadapt in R
cutadapt <- "~pathto directory"
system2(cutadapt, args = "--version")
# Define output files for cutadapt
path.cut <- file.path(path16S, "cutadapt")
if(!dir.exists(path.cut)) dir.create(path.cut)
fnFs.cut <- file.path(path.cut, basename(fnFs))
fnRs.cut <- file.path(path.cut, basename(fnRs))
# Define primers to search for
FWD.RC <- dada2:::rc(FWDPro341F)
REV.RC <- dada2:::rc(REVPro805R)
# Trim FWD and the reverse-complement of REV off of R1 (forward reads)
R1.flags <- paste("-g", FWDPro341F, "-a", REV.RC)
# Trim REV and the reverse-complement of FWD off of R2 (reverse reads)
R2.flags <- paste("-G", REVPro805R, "-A", FWD.RC)
# Run Cutadapt
for(i in seq_along(fnFs)) {
  system2(cutadapt, args = c(R1.flags, R2.flags, "-n", 2, "-o", fnFs.cut[i], "-p", fnRs.cut[i], fnFs.filtN[i], fnRs.filtN[i]))}
Summarized output: 
=== Summary ===

Total read pairs processed:             18,890
  Read 1 with adapter:                  18,502 (97.9%)
  Read 2 with adapter:                  18,739 (99.2%)
Pairs written (passing filters):        18,890 (100.0%)

Total basepairs processed:    11,334,000 bp
  Read 1:     5,667,000 bp
  Read 2:     5,667,000 bp
Total written (filtered):     10,521,217 bp (92.8%)
  Read 1:     5,297,096 bp
  Read 2:     5,224,121 bp

# Check for primers in the first cutadapted sample
rbind(FWD.ForwardReads = sapply(FWD341.orients, primerHits, fn = fnFs.cut[[1]]), 
      FWD.ReverseReads = sapply(FWD341.orients, primerHits, fn = fnRs.cut[[1]]), 
      REV.ForwardReads = sapply(REV805.orients, primerHits, fn = fnFs.cut[[1]]), 
      REV.ReverseReads = sapply(REV805.orients, primerHits, fn = fnRs.cut[[1]]))
After running this, I get the following output:
                 Forward Complement Reverse RevComp
FWD.ForwardReads       1          0       0       0
FWD.ReverseReads       0          0       0       0
REV.ForwardReads       0          0       0       0
REV.ReverseReads       0          0       0       0

So obviously I have a problem, by googling I got the hint I must change the errors/mismatches within the primer sequence, but I can't find how to do this. Any suggestions? Best regards, Veronika Flad

R sequencing • 1.5k views
ADD COMMENT
0
Entering edit mode

I didn't find any reference to FWDPro341F and REVPro805R in your code. Check the cutadapt guide primer trimming guide for more options and tweaks.

https://cutadapt.readthedocs.io/en/stable/recipes.html#trimming-amplicon-primers-from-both-ends-of-paired-end-reads

ADD REPLY

Login before adding your answer.

Traffic: 2304 users visited in the last hour
Help About
FAQ
Access RSS
API
Stats

Use of this site constitutes acceptance of our User Agreement and Privacy Policy.

Powered by the version 2.3.6