Filtering out reads from a reference (e.g. rRNA) using bowtie2

Question

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. --al-conc-gz)
• reads for which one or both of a pair do not align a reference (using e.g --un-conc-gz)

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.

Answer 1

If you want to remove rRNA, use tools that are made for this. We usually go for SortMeRNA, which can do exactly what you describe your goal is, i.e. go from FASTQ to FASTQ, splitting your reads into rRNA and non-rRNA.

In our experience, it is very accurate, if a bit slow at times.

Answer 2

1. Bam flags will indicate which reads have both pairs not mapping.

2. Nothing mapped. You need to figure out why, and you can look at just about anything out figure that out.