Concatenating paired end reads when there are missing reads

[LizBent]

Hi all,

I am looking for a tool that can take Illumina fastq paired end reads (already trimmed and quality filtered, so that not all the read1 sequences are paired with read2 sequences, and vice versa), and concatenate them (by taking the reverse complement of read 2 and attaching it to the end of read1). These reads do not overlap. I'm concatenating them so that when I do database searches (e.g. BLAST), I have more information to use to determine what organism my amplicons came from (this is metagenomics work).

Does anyone have such a tool they would be willing to share? I have zero programming experience and our bioinformatician left months ago.

I've looked at ill2fastq.pl, but it seems to be designed for working with only pairs of reads, and can't handle unpaired reads.

[JackieBadger]

So what you need to do is:
1. Interlace your fastq files, grouping your paired ends by seq ID...i.e. cluster coordinates for read /1 and /2
2.De-interlace the paired reads
3. stitch together these two sequences.

This can easily be done on a local instance of Galaxy (I dont think the web portal has interlacer tool installed).

However, seeing as your paired ends come from opposite ends of sequence I don't see how stitching them together will help you in a BLAST search. You are creating an artificial sequence, and Genbank sequences are individual "real" fragments.
If I were you I wouldn't concatenate the PEs. Use FASTX sequence collapser in Galaxy and batch BLAST your unique reads individually.

[LizBent]

Hi Jackie- actually, if you BLAST the two ends of a sequence (with or without an artificial gap in the middle), you get better matches than if you BLAST just one end at a time. It is possible you'd get a nonsense window where the two ends meet, but the best overall matches would be for the longer ends that match real sequences, so those are the hits that will come out on top.

As for the solution you describe, I was rather hoping to find a script that would allow me to keep track of unpaired read1 and read2 sequences so I can use them as well.

[jbrwn]

do you know how to use the command line at all? post the first read name from each fastq and i'll try to help you out. my solution will require python.

[JackieBadger]

The solution I posted tracks unpaired reads.
otherwise check out here http://sfg.stanford.edu/quality.html
Their PECombiner.sh has a bug in it...they may have updated this on the site?
If they have not ask the authors to send you the working script.

[jbrwn]

uses a good amount of memory since it's storing one fq in a dict, but seems to work:

https://gist.github.com/3861828

edit: you didn't say anything about preserving the quals, so this prints a fasta.

[LizBent]

Thanks so much, I will try it