How to properly 'filter' with bowtie via --un?

[rxzlmn]

Hi,

I would like to filter my small RNA-seq data vs. tRNA/rRNAs. However, after I did the filtering, I still end up with reads there when I count.

What I did is the following:

1. Download the tRNA/rRNA subset from UCSC by selecting repClass = 'tRNA' and 'rRNA'

2. Create a FASTA file from that using bedtools getfasta

3. Map my FASTQ files to a bowtie index built from those FASTA files using less strict criteria than I use for my actual mapping (-v 3, ie. 3 mismatches allowed), write unmapped reads to a new FASTQ via --un

4. Map those unmapped reads to the genome using more strict criteria (-v 0)

5. Count reads on the tRNA/rRNA bed files that I originally downloaded in R using summarizeOverlaps with mode='IntersectionStrict'

Albeit few, I still end up with reads that are being counted there. I wonder how is that possible? What's a possible reason for that?

Any help's greatly appreciated

[Cytosine]

Hard to say... could be that the UCSC sets are not complete, or maybe a high sequence redundancy between ncRNAs?

I'd recommend SortMeRNA for getting rid of rRNA reads in samples though.

[rxzlmn]

I figured it out myself, it was a problem with 0-based coordinates from the UCSC data.

After reading the bed files through read.table with makeGRangesFromDataFrame and starts.in.df.are.0based=TRUE, no more hits are being reported.

[dpryan]

Out of curiousity, do the reads that map to tRNA/rRNA after filtering have more than 3 mismatches? Regardless, one thing to keep in mind is that there is often some randomness in the results of various aligners (e.g., if you reorder the reads, you'll often get slightly different alignments). That may also be the cause of what you're seeing.

Edit: I guess you figured out the problem, nevermind!