Merging paired reads with SeqPrep

[jagodoy]

I'm trying to use SeqPrep to merge MiSeq paired reads and I observed many of the alignments show the five terminal bases of the second read as "-" and these are excluded from the merged sequences. These do not seem to be really bad quality calls (34-36 in the example below), so I don't quite understand why this happens. We are dealing with authentically short fragments so losing this five bases will have a noticable impact on the global output. Any clue will be most appreciated ...

-Jose.

An example alignment along with the corresponding raw reads is shown below.

Read Alignment Score:107, Suboptimal Score:-82
ID:M00476:52:000000000-A4KH3:1:1101:18272:1674 1:N:0:1
READ1: TCAGATTCTGTGTCTCCCTCTCTCTGACCCTCCCTTATTCATGCTCTGTCTCTCTCTGTTTCAAAAATACATAAGTGTTAAAAG
|||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
READ2: -----TTCTGTGTCTCCCTCTCTCTGACCCTCCCTTATTCATGCTCTGTCTCTCTCTGTTTCAAAAATACATAAGTGTTAAAAG
MERGD: TTCTGTGTCTCCCTCTCTCTGACCCTCCCTTATTCATGCTCTGTCTCTCTCTGTTTCAAAAATACATAAGTGTTAAAAG

a1_S1_L001_R1_001.fastq.gz:@M00476:52:000000000-A4KH3:1:1101:18272:1674 1:N:0:1
a1_S1_L001_R1_001.fastq.gz:TCAGATTCTGTGTCTCCCTCTCTCTGACCCTCCCTTATTCATGCTCTGTCTCTCTCTGTTTCAAAAATACATAAGTGTTAAAAG
a1_S1_L001_R1_001.fastq.gz:+
a1_S1_L001_R1_001.fastq.gz:CCDCDFFFFFFFGGGGGGGGGGHHHHHHHHHHGGHHHHHHHHHHHHHHHHHHHHHHHIHIHIIIIHHHHHHHHHHHHHHHHHHH

a1_S1_L001_R2_001.fastq.gz:@M00476:52:000000000-A4KH3:1:1101:18272:1674 2:N:0:1
a1_S1_L001_R2_001.fastq.gz:CTTTTAACACTTATGTATTTTTGAAACAGAGAGAGACAGAGCATGAATAAGGGAGGGTCAGAGAGAGGGAGACACAGAATCTGAGGAAGAGCGTCGTGTAGGGAAAGAGTGTAGATCTCGGTGGTCGCCGTATCATTAAAAAAAAAAAAA
a1_S1_L001_R2_001.fastq.gz:+
a1_S1_L001_R2_001.fastq.gz:CEECDFFFFFFFGGGGGGGGGGGGHHHHHHHHGGGGHHHHHHHHHHHHHHHHHHGGGGHHHHHHGGGGGGGGGHHHHHHHGHHHHGHHHHHGGGGGGGHHHHGHHHHHHGGGHHHHEHHHFFGGGHCGGGG?FGHHHHGHGH@>AF@<-<

[GenoMax]

Have you tried FLASH to see if that suits your needs better?

[mastal]

Merging paired reads with SeqPrep

Hi,

I notice that your raw reads R1 and R2 are different lengths, so they have presumably already been trimmed to remove poor quality ends or adapters.

Is SeqPrep trying to also do some adapter trimming?

What version of Illumina sample prep kit was used for library prep (i.e. what adapter sequences are possibly present in your reads)?

[jagodoy]

Have you tried FLASH to see if that suits your needs better?

Thanks, GenoMax, for the advice. I'll give FLASH a try if I cannot solve this with SeqPrep.

[jagodoy]

Mastal, because two slightly different libraries were run simultaneously, R1 read was already trimmed from adaptors, but R2 reads still show them. I asked SeqPrep to remove adaptors and this was done efficiently, leaving R1 and R2 of exactly the same length, as in the example.

[mastal]

Merging paired reads with SeqPrep

Mastal, because two slightly different libraries were run simultaneously, R1 read was already trimmed from adaptors, but R2 reads still show them.

with respect to Illumina paired reads, this doesn't make any sense, as R1 and R2 come from the same fragment.

I asked SeqPrep to remove adaptors and this was done efficiently, leaving R1 and R2 of exactly the same length, as in the example.

the R2 that you showed in the example is 5 bases shorter than R1, I think that is because the 5 bases are part of the Illumina adapter.

[jagodoy]

The reason why R2 still has the adapter is because it is not the standard adapter for which the run was trimmed, while the other R1 adapter was standard and thus properly trimmed.

Note that once R2 is trimmed, R1 and R2 are the same size and completely complementary over their full length.

[jagodoy]

OK. I think I just realized what maybe going on. Our modified adapter is 5 bases shorter, but other than that identical to the standard. I'm providing the truncated real adapter sequence using the option -B, but the program seems to be applying the default standard adapter and deleting the five bases upstream the matched adapter sequence, even though they do not match at all.

If this is the case, I might not be using the custom adapters option correctly, or it might not be functioning properly. This an example of how I'm running SeqPrep:

SeqPrep -f a1_S1_L001_R1_001.fastq.gz -r a1_S1_L001_R2_001.fastq.gz -1 a1_S1_merged_R1.fastq.gz -2 a1_S1_merged_R2.fastq.gz -3 a1_S1_merged_discard_R1.fastq.gz -4 a1_S1_merged_discard_R2.fastq.gz –B GGAAGAGCGTCGTGTAGGGA –g -s a1_S1_merged.fastq.gz -E a1_S1_merged.align.gz –o 11


Any clue?