Mispaired reads following trimming and demultiplexing

[Tony Kess]

Hello,

I'm using cutadapt to trim off the leading 2 basepairs from reads from several bestRAD libraries to enable finding RAD loci as all reads start with NA/G/T/C, as well as remove over-represented sequence identified using fastQC. When I attempt to demultiplex these libraries I find that individual reads are now out of sync (different read names).  Alignment in bwa mem gives the error:

[mem_sam_pe] paired reads have different names: "489_7_1208_10581_48597", "489_7_1208_8430_48597"

Checking of read pairing for demultiplexed samples in cutadapt gives this output: 

reads/minutecutadapt: error: Error in FASTQ file at line 47956: Length of sequence and qualities differ

Has anyone encountered this problem before? 

I'm using stacks 2.3, and cutadapt 2.1.

Here is the command used for trimming in cut adapt:

cutadapt \

  -u 2 \

  -U 2 \

  --length 90 \

  -a TCGGAAGAGCACACGTCTGAACTCCAGTCACCGATGTATCTCGTATGCCG \

  -A TCGGAAGAGCGTCGTGTAGGGAAAGAGTGTAGATCTCGGTGGTCGCCGTA \

  -o HI.5003.006.NEBNext_Index_14.HhiRAD14_R1_trim.fastq \

  -p HI.5003.006.NEBNext_Index_14.HhiRAD14_R2_trim.fastq \

  HI.5003.006.NEBNext_Index_14.HhiRAD14_R1.fastq.gz \

  HI.5003.006.NEBNext_Index_14.HhiRAD14_R2.fastq.gz \

  -j 12

And here is the command used for demultiplexing

   process_radtags -1 HI.5003.006.NEBNext_Index_14.HhiRAD14_R1_trim.fastq \

  -2 HI.5003.006.NEBNext_Index_14.HhiRAD14_R2_trim.fastq \

  -o Lib14 \

  -b Halibut_Barcodes_Lib14.txt \

  -r  --bestrad  \

  -c -q \

  -e 'sbfI' \

  --barcode_dist_1 2 \

  --adapter_1 GATCGGAAGAGCACACGTCTGAACTCCAGTC \

  --adapter_2 CACTCTTTCCCTACACGACGCTCTTCCGATCT 

[Amanda Stahlke]

I'm interested in this because we also are using the bestrad flag and have some issues when the usual GG overhang has errors at the start of a read. I hadn't tried just trimming the reads - could be a good solution!

You could always try another trimmer (trimmomatic?) and/or alignment software.  It sounds like cutadapt could be the root of the issue. 

Following!

On Thursday, April 18, 2019 at 7:50:39 AM UTC-7, Tony Kess wrote

[Tony Kess]

The most recent alignments in bwa mem worked without a mis-pairing error. It turns out that cutadapt will keep reads of 0 length, but it appears process RADtags wasn't incorporating them. Importantly, the --length option in cutadapt doesn't override the retention of 0 basepair reads. Here is our command for cut adapt that worked:

cutadapt \

  -u 2 \

  -U 2 \

  --length 90 \

  --minimum-length 90 \

  -a TCGGAAGAGCACACGTCTGAACTCCAGTCACCGATGTATCTCGTATGCCG \

  -a GATCGGAAGAGCACACGTCTGAACTCCAGTC \

  -A CACTCTTTCCCTACACGACGCTCTTCCGATCT \

  -A TCGGAAGAGCGTCGTGTAGGGAAAGAGTGTAGATCTCGGTGGTCGCCGTA \

  -o {}_R1_trim.fastq \

  -p {}_trim.fastq \

  HI.5003.006.NEBNext_Index_14.HhiRAD14_R1.fastq.gz \

  HI.5003.006.NEBNext_Index_14.HhiRAD14_R2.fastq.gz \

SO we just needed to specify a minimum length option as well. Hope this is helpful in the future!