As I understand it, bowtie2 can easily be used to split reads into one of two groups:
- reads for which both of a pair align well to a reference (using e.g.
- reads for which one or both of a pair do not align a reference (using e.g
But I really want to split this second group into reads for which neither of a pair align to the reference.
bowtie2's report give a nice breakdown, but it's not clear to me how to extract specific members of these groups:
1000000 reads; of these: 1000000 (100.00%) were paired; of these: 958118 (95.81%) aligned concordantly 0 times 1329 (0.13%) aligned concordantly exactly 1 time 40553 (4.06%) aligned concordantly >1 times ---- 958118 pairs aligned concordantly 0 times; of these: 179 (0.02%) aligned discordantly 1 time ---- 957939 pairs aligned 0 times concordantly or discordantly; of these: 1915878 mates make up the pairs; of these: 711518 (37.14%) aligned 0 times 32134 (1.68%) aligned exactly 1 time 1172226 (61.18%) aligned >1 times 64.42% overall alignment rate
In this case, I really want those that didn't align at all (i.e. the 35.58% not included in the
64.42% overall alignment rate which is the same as the
37.14% of the pairs that did not align concordantly).
Since in my case the reference is rRNA and many of those sequences are shorter than my reads, I'm not surprised that so there are so few concordant reads. I'm really interested in reads that aren't concordant and aren't discordant discordant either.
So, if I'm sticking to bowtie2, it looks like my only option is to take the resulting SAM/BAM file and filter what I want based on flags and then convert that back into FASTQ files.
Is there another option that I'm missing? It seems that the remaining output options relate to unpaired reads or SAM files.
I'm looking into other aligners now...
My goal is to go from FASTQ files to FASTQ files without extra steps, but I can go through a SAM/BAM file stage if I need to and even build a solution around that if needed.