cannot remove primers on reverse reads using cutadapt

[kdbchau]

Cross-posted from cutadapt github because the author hasn't responded to other posts in a while...my issue is that for my ITS sequences, I am wanting to use cutadapt to remove primers but I cannot remove them entirely from my reverse reads no matter what I do.

Here is my code and output (I do have to run cutadapt on the terminal because I have had issues with it on qiime2 on our cluster and on R, but the flags and primers would essentially be identical):

# Primers
FWD <- "CTTGGTCATTTAGAGGAAGTAA"
REV <- "GCTGCGTTCTTCATCGATGC"

> FWD.orients <- allOrients(FWD); FWD.orients
                 Forward               Complement                  Reverse                  RevComp 
"CTTGGTCATTTAGAGGAAGTAA" "GAACCAGTAAATCTCCTTCATT" "AATGAAGGAGATTTACTGGTTC" "TTACTTCCTCTAAATGACCAAG" 

> REV.orients <- allOrients(REV); REV.orients
               Forward             Complement                Reverse                RevComp 
"GCTGCGTTCTTCATCGATGC" "CGACGCAAGAAGTAGCTACG" "CGTAGCTACTTCTTGCGTCG" "GCATCGATGAAGAACGCAGC"

# Check primers in original uncut/unfiltered files; not filtering because I lose way too many reads
rbind(FWD.ForwardReads = sapply(FWD.orients, primerHits, fn = fnFs[[1]]),
      FWD.ReverseReads = sapply(FWD.orients,primerHits, fn = fnRs[[1]]),
      REV.ForwardReads = sapply(REV.orients, primerHits,fn = fnFs[[1]]),
      REV.ReverseReads = sapply(REV.orients, primerHits, fn = fnRs[[1]]))

# Output
                 Forward Complement Reverse RevComp
FWD.ForwardReads   27211          0       0       0
FWD.ReverseReads    1720       1715    1707   18415
REV.ForwardReads       0          0       0   26660
REV.ReverseReads   26766       1826    1728    1744

FWD <- FWD.orients[["Forward"]]; FWD  # keep the same as is
REV <- REV.orients[["RevComp"]]; REV  # because it matches to the revcomp

> REV.orients <- allOrients(REV); REV.orients
               Forward             Complement                Reverse                RevComp 
"GCTGCGTTCTTCATCGATGC" "CGACGCAAGAAGTAGCTACG" "CGTAGCTACTTCTTGCGTCG" "GCATCGATGAAGAACGCAGC"

I then run cutadapt in terminal and normally this has never caused an issue with my single-end data for ITS, but this time I am working with paired-end and have issues with removing primers for reverse reads.

for r1 in *_R1.fastq; do r2=${r1/_R1.fastq/_R2.fastq}; sample=${r1%%_R1.fastq}; /home/chauk/.local/bin/cutadapt \
-a CTTGGTCATTTAGAGGAAGTAA...GCATCGATGAAGAACGCAGC -A GCTGCGTTCTTCATCGATGC...TACTTCCTCTAAATGACCAAG \
-n 2 -o original/${sample}_R1_cut1.fastq -p original/${sample}_R2_cut1.fastq $r1 $r2; done

Then I go back into R and run my sanity check which gives me the following output:

                 Forward Complement Reverse RevComp
FWD.ForwardReads       0          0       0       0
FWD.ReverseReads    1720       1715    1707    1518
REV.ForwardReads       0          0       0       0
REV.ReverseReads    1834       1826    1728    1744

Now, I've tried running the -A flag with different orientations for the reverse primer, and even did multiple rounds of primer removal but it doesn't change the fact that these original primers are still showing up in every orientation only in the reverse reads (even though a good chunk were removed from the reverse reads in the RevComp orientation and in the reverse reads in the Forward orientation.

I have done multiple rounds after seeing that the sequencing company has added adaptors to the primers and did it 4 times both for the forward and reverse primers - HOWEVER, it still works completely fine for forward reads where just using the actual primer sequence I can remove it entirely. Below are the additional adaptor+primer sequences used:

Forward1 = ACACTCTTTCCCTACACGACGCTCTTCCGATCTCTTGGTCATTTAGAGGAAGTAA
Forward2 = ACACTCTTTCCCTACACGACGCTCTTCCGATCTTCTTGGTCATTTAGAGGAAGTAA
Forward3 = ACACTCTTTCCCTACACGACGCTCTTCCGATCTACCTTGGTCATTTAGAGGAAGTAA
Forward4 = ACACTCTTTCCCTACACGACGCTCTTCCGATCTCAACTTGGTCATTTAGAGGAAGTAA

Reverse1 = GTGACTGGAGTTCAGACGTGTGCTCTTCCGATCTGCTGCGTTCTTCATCGATGC
Reverse2 = GTGACTGGAGTTCAGACGTGTGCTCTTCCGATCTTGCTGCGTTCTTCATCGATGC
Reverse3 = GTGACTGGAGTTCAGACGTGTGCTCTTCCGATCTACGCTGCGTTCTTCATCGATGC
Reverse4 = GTGACTGGAGTTCAGACGTGTGCTCTTCCGATCTCATGCTGCGTTCTTCATCGATGC

I am at a loss at what to do with this. Any advice would be appreciated!

[kdbchau]

Update. I ran a second round of cutadapt on these trimmed reads using the following command to see if I can just remove these reads with adapters still present - it made zero difference. I don't understand where these adaptor sequences are being detected in my reverse reads.

cd original
for r1 in *_R1_cut1.fastq; do r2=${r1/_R1_cut1.fastq/_R2_cut1.fastq}; sample=${r1%%_R1_cut1.fastq}; /home/chauk/.local/bin/cutadapt -a CTTGGTCATTTAGAGGAAGTAA...GCATCGATGAAGAACGCAGC -A GCTGCGTTCTTCATCGATGC...TACTTCCTCTAAATGACCAAG --pair-filter=both --discard-trimmed -o round2/${sample}_R1_cut2.fastq -p round2/${sample}_R2_cut2.fastq $r1 $r2; done


# Sanity check in R shows:

# go to path where the first set of trimmed reads go
path.cut2 <- "~/../Desktop/.../original/round2"
myfiles.For.cut2 <- file.path(path.cut, paste0(strsplit(basename(fnFs),".fastq"),"_cut2.fastq"))
myfiles.Rev.cut2 <- file.path(path.cut, paste0(strsplit(basename(fnRs),".fastq"),"_cut2.fastq"))


# Sanity check; call the original primers and check if anymore are removed
rbind(FWD.ForwardReads = sapply(FWD.orients, primerHits, fn = myfiles.For.cut2[[1]]), 
      FWD.ReverseReads = sapply(FWD.orients, primerHits, fn = myfiles.Rev.cut2[[1]]), 
      REV.ForwardReads = sapply(REV.orients, primerHits, fn = myfiles.For.cut2[[1]]), 
      REV.ReverseReads = sapply(REV.orients, primerHits, fn = myfiles.Rev.cut2[[1]]))

# Output; NO DIFFERENCE
                 Forward Complement Reverse RevComp
FWD.ForwardReads       0          0       0       0
FWD.ReverseReads    1720       1715    1707    1518
REV.ForwardReads       0          0       0       0
REV.ReverseReads    1744       1728    1826    1834

I did try with --discard-untrimmed but I lost all my reads.

[benjjneb]

Do any of your original sequences contain ambiguous nucleotides in them? Like "N"?

[kdbchau]

Yes - while not all reads, there are a bunch in reverse reads that look like the following:

+
IIIIIIIIIIIIIIIIIIIIIIIIII9IIII#-##-I#-##I--99#I#99###########-###########--#############-9I9999-9---999----9----99###9#-###-#--#-999--#----##-##-9-##9#9----#----#--9#-9---#-#-9----9##9-99--9#####-#####99--#----####99--#9-#---#9#---9--99##--##-9##--#9-###9##-######-#-#-----#9##--###-#99--#-##-#-###9#
@VH01519:242:AAHMVKKM5:1:1304:56614:1435 2:N:0:GTGTATGCGT+AGTGCTCCAT
TGCTGCGTTCTTCATCGATGCCAGAACCAAGNCNNCCNCNNCCCCCNCNCTNNNNNNNNNNNGNNNNNNNNNNNCCNNNNNNNNNNNNNCTTCTTCTCCTCTCGTCCTTCTTCGTNNNTNTNNNTNTTNTCCCTCNTTCCNNTNNTTTNNTNCCTCCNCCGGNTTCNTTCCTNCNCTGCTTCNNTCTTTCGNNNNNCNNNNNCTCGNGCCTNNNNTTCGNGCNTTCNTNTCTTCGCTNNTTNNCGNNTTNCTNNNTNNCNNNNNNTNTNTCTATNTNNTTNNNCNTTCTNCNNTNCNNNTN
+
IIIIIIIIIIIIIIIIIIIIIIIIIIIIIII#9##--#-##9----#-#--###########-###########9-#############-99-9--9-9999------9--9--9###-#-###9#99#-9-99-#-9--##9##--9##9#-99--#9999#-9-#9--99#-#9--999-##9-----9#####9#####999-#9-9-####---9#-9#9--#-#-9--9999##--##--##-9#--###-##9######-#9#-9--9#-##--###-#-9-9#9##-#-###9#
@VH01519:242:AAHMVKKM5:1:1304:22908:1511 2:N:0:GTGTATGCGT+AGTGCTCCAT
ACGCTGCGTTCTTCATCGATGCCAGAACCAANANNTTNGNNGTTAANGNTTNNNNNNNNNNNANNNNNNNNNNNCANNNNNNNNNNNNNCGCAAGAGTTTGAAGTGGCCAGCCGGNNNANCNNNCNGANGGCAGGNTCACNNANNGGGNNANAGAGTNGGAANACANAGGGTNANGGACATCNNGGTGGTTNNNNNGNNNNNGATTNTAATNNNNCATGNGCNGGTNTNTTTAGGGANNAANNGTNNAGNATNNNANNTNNNNNNANGNAAAAANGNNGANNNGNAGAGNGNNGNGNNNGN
+
IIIIIIIIIIIIIIIIIIIIII-IIIIIIII#9##I9#9##II999#-#9I###########9###########9-#############9999-99I99II9-9I--9-9--999###-#-###-#-9#II9999#I9I-##9##99I##9#9999-#I999#999#I9I9I#-#9-9999-##999I99I#####9#####-9--#9-9-####9---#99#9--#-#--9-----##9-##9-##-9#--###9##-######-#-#9---9#9##--###-#9---#9##-#9###9#
@VH01519:242:AAHMVKKM5:1:1304:22719:1700 2:N:0:GTGTATGCGT+AGTGCTCCAT
ACGCTGCGTTCTTCATCGATGCCAGAACCAANANNTCNGNNGTTAGNGNTTNNNNNNNNNNNANNNNNNNNNNNCANNNNNNNNNNNNNCGCAAGAGTTTGAAGTGGCCAACCGANNNCNCNNNCNGANGGCAGGNTCGCNNANNGGGNNANAGAGTNGGAANACGNAGGGTNANGAACATCNNGGTGGTTNNNNNGNNNNNAATTNTAATNNNNCCTGNGCNGGTNCNAGTAAGGANNCTNNGTNNTGNCTNNNANNTNNNNNNGNTNAAAGANGNNGGNNNGNAGAGNGNNGNGNNNGN
+
II9IIIIII9IIIIIIIIIIIII-9IIIIII#9##9-#-##I999-#9#9-###########9###########--#############999-99-I99I-9-9-9--9--999-###-#-###-#-9#II9999#9999##-##--9##9#-9999#9999#99-#I9I9I#9#9-9-999##9-II999#####9#####99-9#9-9-####I-99#9-#9-9#-#-9-9-9-9##9-##99##--#99###-##9######-#9#-9-99#-##9-###-#9---#9##9#9###9#
@VH01519:242:AAHMVKKM5:1:1304:64869:1700 2:N:0:GTGTATGCGT+AGTGCTCCAT
TGCTGCGTTCTTCATCGATGCCAGAACCAAGNGNNCCNCNNCCACCNANGCNNNNNNNNNNNTNNNNNNNNNNNCCNNNNNNNNNNNNNTTTGCCCTCATGTACTCTTCAGCCGTNNNTNTNNNTNTTNCCCGCCNTTACNNGNNGTCNNTNCCCCCNCCGGNTCTNTCCCCNCNACTCTGGNNTCGTCTGNNNNNCNNNNNTTCGNGGTCNNNNTCCCNGGNTCGNTNCTGCCGCTNNTTNNTTNNCTNCTNNNTNNTNNNNNNCNCNTAGGGNGNNGANNNCNGCCGNCNNTNCNNNTN
+
IIIIII9IIIIIIIIII9IIIIIIIII9I9I#-##--#-##I----#-#--###########-###########-9#############---9-------9-999-99--9-999###-#-###-#-9#--9---#9---##-##--9##9#-----#--99#9--#---99#-#---9---##---9--9#####-#####--9-#--99####----#9-#---#-#-9--9--9##9-##99##-9#-9###9##-######9#9#-----#9##--###-#--9-#-##9#-###-#
@VH01519:242:AAHMVKKM5:1:1304:20163:1947 2:N:0:GTGTATGCGT+AGTGCTCCAT
TGCTGCGTTCTTCATCGATGCCAGAACCAAGNANNCCNTNNTCGAANGNTANNNNNNNNNNNTNNNNNNNNNNNAGNNNNNNNNNNNNNAAAGGTTTATGTTTGTAATAGTGATGNNNGNANNNANCANGGAAACNAGGANNANNGAANNGNCAAGGNTGAANAATNCAGCANGNTTGTAACNNAGAGGGGNNNNNGNNNNNATTCNTATTNNNNTAATNATNCATNCNGGGGTTCANNTANNGANNACNTGNNNANNCNNNNNNTNANTATGANTNNCTNNNANGTGGNANNANTNNNGN
+
IIIIIIIIIIIIIIIIIIIIIIIIIIIIIII#-##--#-##9-9-9#I#9-###########9###########--#############--I99999-9I9999I--I-99--9-###9#9###9#-9#II999-#I9I-##I##--I##-#--9-9#-999#9--#I99-I#9#--9999-##-9-999-#####9#####9-9-#--99####9-9-#II#--9#-#--99-999##-I##99##--#-9###9##-######-#-#--999#-##--###-#9---#9##-#-###9#
@VH01519:242:AAHMVKKM5:1:1304:48339:1984 2:N:0:GTGTATGCGT+AGTGCTCCAT
GCTGCGTTCTTCATCGATGCCAGAACCAAGANGNNCTNTNNCCTCGNCNAANNNNNNNNNNNGNNNNNNNNNNNCCNNNNNNNNNNNNNCTTCCTCTTCGCTTTGCTATCGTTGCNNNGNTNNNTNTTNTTCTTCNTTTTNNTNNGTTNNGNCCGTCNCCGCNCGCNCCCCTNCNCCCCAGCNNGGCGTCTNNNNNTNNNNNGGCGNCGAGNNNNTTTGNGCNCGCNCNTCTGTGTCNNTCNNTTNNTCNCCNNNTNNTNNNNNNTNCNTTCTTNTNNTTNNNTNTTCGNTNNTNCNNNGN
+
I9IIIIIIIIIIIIIIIIIIIIIIIIIIIII#-##9-#9##--9--#9#--###########-###########9-#############9----9--9--9--------9---9-###-#-###-#--#99---9#9999##-##9--##9#--99-#--99#9--#---9-#-#--9--99##--99--9#####-#####----#--99####---9#-9#99-#-#9-9--9--##9-##--##-9#-9###-##-######-#-#----9#-##9-###-#----#-##-#-###-#
@VH01519:242:AAHMVKKM5:1:1304:65210:2041 2:N:0:GTGTATGCGT+AGTGCTCCAT
CATGCTGCGTTCTTCATCGATGCCAGAACCANANNACNCNNCCACCNCNGCNNNNNNNNNNNGNNNNNNNNNNNCCNNNNNNNNNNNNNGTTGTCTTAAGGTCCTCTTCAGCAGTNNNTNTNNNTNTTNCCCGCCNTCGCNNGNNCTCNNTNCTCCCNCCGGNTCCNTCTCCNCNACTCTCTNNCCCTCTGNNNNNCNNNNNTCCGNGGTCNNNNTCTCNGGNTCTNTNCCTCCACTNNTTNNTTNNCTNTTNNNTNNTNNNNNNCNCNTTTGGNGNNGGNNNCNGATGNCNNTNCNNNCN

[kdbchau]

I will mention that after running cutadapt once and then running the filterAndTrim I did manage to remove such reads and I ended up with the following outcome:

out <- filterAndTrim(cut.files, filt.reads, maxN = 0,
                     minLen = 100, rm.phix = TRUE, compress = FALSE, multithread = FALSE)

I then ran the sanity check again on these cutadapt-ed and filtered files and now it looks like all the primers are removed (except a couple samples that have 1 primer showing up in a read):

Some examples of output for different samples:

# Sample 1
                 Forward Complement Reverse RevComp
FWD.ForwardReads       0          0       0       0
FWD.ReverseReads       0          0       0       0
REV.ForwardReads       0          0       0       0
REV.ReverseReads       0          0       0       0

# Sample 10
                 Forward Complement Reverse RevComp
FWD.ForwardReads       0          0       0       0
FWD.ReverseReads       0          0       0       0
REV.ForwardReads       0          0       0       0
REV.ReverseReads       1          0       0       0

So I think, considering that I will use the trimmed and filtered reads anyways that this should be fine that I don't see all the primers removed right after the initial cutadapt step right? It looks like those problematic reads get removed anyways after this step.