Problem trimming primers

[Gian77]

Hello,
I am trying to trim the primer off of 16S Miseq reads ( I am using just 4 samples, 8 fastq PE files as a toy dataset ) and I am having trouble understanding how the removePrimers function works.

Perimers are present as you can see, when I search for them

> CountPrimer <- function(primer, file_names){
+   require(ShortRead)
+   
+   num_sequences <- vector("integer", length(file_names))
+   
+   for (i in 1:length(file_names)) {
+     file_path <- file_names[i]
+     print(file_path)
+     num_hits <- 
+       vcountPattern(primer, ShortRead::sread(readFastq(file_path)), 
+                     with.indels = TRUE,
+                     max.mismatch=0,
+                     fixed = FALSE)
+     num_sequences[i] <- sum(num_hits > 0)
+   }
+   return(num_sequences)
+ }
> CountPrimer(FWD_primer, fastq_names)
[1] "/home/gian/Dropbox/2_BIOINFORMATICS/bioinfo_pipelines/pipe_amplicon_DADA2/rawdata/col-10nbd_S115_L001_R1_001.fastq"
[1] "/home/gian/Dropbox/2_BIOINFORMATICS/bioinfo_pipelines/pipe_amplicon_DADA2/rawdata/col-10nbd_S115_L001_R2_001.fastq"
[1] "/home/gian/Dropbox/2_BIOINFORMATICS/bioinfo_pipelines/pipe_amplicon_DADA2/rawdata/col-10scd_S294_L001_R1_001.fastq"
[1] "/home/gian/Dropbox/2_BIOINFORMATICS/bioinfo_pipelines/pipe_amplicon_DADA2/rawdata/col-10scd_S294_L001_R2_001.fastq"
[1] "/home/gian/Dropbox/2_BIOINFORMATICS/bioinfo_pipelines/pipe_amplicon_DADA2/rawdata/col-10scw_S206_L001_R1_001.fastq"
[1] "/home/gian/Dropbox/2_BIOINFORMATICS/bioinfo_pipelines/pipe_amplicon_DADA2/rawdata/col-10scw_S206_L001_R2_001.fastq"
[1] "/home/gian/Dropbox/2_BIOINFORMATICS/bioinfo_pipelines/pipe_amplicon_DADA2/rawdata/col-11nbw_S55_L001_R1_001.fastq"
[1] "/home/gian/Dropbox/2_BIOINFORMATICS/bioinfo_pipelines/pipe_amplicon_DADA2/rawdata/col-11nbw_S55_L001_R2_001.fastq"
[1] 22511     0  8552     0  1074     0 14210     0
> CountPrimer(REV_primer, fastq_names)
[1] "/home/gian/Dropbox/2_BIOINFORMATICS/bioinfo_pipelines/pipe_amplicon_DADA2/rawdata/col-10nbd_S115_L001_R1_001.fastq"
[1] "/home/gian/Dropbox/2_BIOINFORMATICS/bioinfo_pipelines/pipe_amplicon_DADA2/rawdata/col-10nbd_S115_L001_R2_001.fastq"
[1] "/home/gian/Dropbox/2_BIOINFORMATICS/bioinfo_pipelines/pipe_amplicon_DADA2/rawdata/col-10scd_S294_L001_R1_001.fastq"
[1] "/home/gian/Dropbox/2_BIOINFORMATICS/bioinfo_pipelines/pipe_amplicon_DADA2/rawdata/col-10scd_S294_L001_R2_001.fastq"
[1] "/home/gian/Dropbox/2_BIOINFORMATICS/bioinfo_pipelines/pipe_amplicon_DADA2/rawdata/col-10scw_S206_L001_R1_001.fastq"
[1] "/home/gian/Dropbox/2_BIOINFORMATICS/bioinfo_pipelines/pipe_amplicon_DADA2/rawdata/col-10scw_S206_L001_R2_001.fastq"
[1] "/home/gian/Dropbox/2_BIOINFORMATICS/bioinfo_pipelines/pipe_amplicon_DADA2/rawdata/col-11nbw_S55_L001_R1_001.fastq"
[1] "/home/gian/Dropbox/2_BIOINFORMATICS/bioinfo_pipelines/pipe_amplicon_DADA2/rawdata/col-11nbw_S55_L001_R2_001.fastq"
[1]     0 16873     0  8249     0  1074     0 13498

When I try to trime them off the reads, it seem no primers is detected

> fastq_filt <- dada2::removePrimers(fastq_names, 
+                                    fastq_no_primer, 
+                              primer.fwd = FWD_primer, 
+                              primer.rev = REV_COMP_primer, 
+                              orient = TRUE,
+                              compress = TRUE,
+                              verbose = TRUE,
+                              max.mismatch = 0,
+                              allow.indels = TRUE)
Primer matching with indels allowed is currently significantly (~4x) slower.
Read in 23980, output 0 (0%) filtered sequences.
Read in 23980, output 0 (0%) filtered sequences.
Read in 9340, output 0 (0%) filtered sequences.
Read in 9340, output 0 (0%) filtered sequences.
Read in 1516, output 0 (0%) filtered sequences.
Read in 1516, output 0 (0%) filtered sequences.
Error in sapply(match.fwd, end) + 1 : 
  non-numeric argument to binary operator
> fastq_filt
                                 reads.in reads.out
col-10nbd_S115_L001_R1_001.fastq    23980         0
col-10nbd_S115_L001_R2_001.fastq    23980         0
col-10scd_S294_L001_R1_001.fastq     9340         0
col-10scd_S294_L001_R2_001.fastq     9340         0
col-10scw_S206_L001_R1_001.fastq     1516         0
col-10scw_S206_L001_R2_001.fastq     1516         0
col-11nbw_S55_L001_R1_001.fastq     15565         0
col-11nbw_S55_L001_R2_001.fastq     15565         0

Thanks much,

Gian

[benjjneb]

A couple things.

We don't recommend using removePrimers on Illumina data. It was developed to solve some specific issues with long-read sequencing data, but it is not very performant with read numbers, and there are other more fully featured methods like cutadapt or trimmomatic for custom primer detection and removal.

That said, in most cases the best solution is to use filterAnd Trim(..., trimLeft=c(FWD_PRIM_LEN, REV_PRIM_LEN) to remove primers. This works as long as the primers are all a constant length, and always at the start of the forward and reverse reads respectively (which is usually the case for the common amplicon library strategies).

As to what you are seeing, it seems that you may have replaced REV_primer the REV_COMP_primer in your removePrimers call. In the code above it was REV_primer that was being detected in many of the R2 reads. Probably REV_COMP_primer is not appearing in the R2 reads at all, and thus no reads are passing removePrimers.

[Gian77]

Thank you so much for the answer, I will use cutadapt.
Gian