aligners filtering out reads with excessive read depth?

[efoss]

I have RNA seq data for total yeast RNA (not poly-A purified). 90% of my reads should be ribosomal RNA. When aligning with bwa and then looking at read depth at the rDNA using bedtools' genomeCoverageBed or IGV, I see little to no coverage of the rDNA - certainly far far less than I should be seeing. If I go into the original fastq file and just pull out sequences at random and blast them to the yeast genome, almost all of them map to the rDNA, as they should. I imagine that when bwa sees greater than a certain read depth that it just tosses those reads. (I can't think of another way to explain this.) Looking at the help pages for bwa and gsnap (I'd be happy to use either), I don't see an explanation of this nor an option to turn it off. Can anyone suggest what is going on and what I can do to fix it? I want to see all of the reads.

Thanks.

Eric

[swbarnes2]

bwa doesn't throw out reads. I bet you can confirm that by counting the # of lines in your fastq with the # of lines in your .bam file,.

Well, it looks like genomeCoverageBed has a max setting, and IGV probably does too.

I bet if you used samtools view to filter down to the reads that align in the rRNA region of the genome, and then counted how many reads were in that .bam, you'd see what you expect.

[efoss]

Hi swbarnes2,

Thanks very much. This is really helpful.

Best wishes,

Eric