RNA expression - reads mapping outside of the annotated regions of the genome

[feralBiologist]

I have data from an experiment which is looking into differential expression of small RNA. It seems that a fraction of the reads are mapped outside of the annotated regions of the mouse genome (as listed in the Ensemble GTF file). To be honest, I am not expecting to find new genes in such a popular model organism but I am still curious where these reads map. I would like to produce statistics for each chromosome listing the number of reads that are mapped outside of the annotated regions. What's the easiest way to produce such statistics? I am using bowtie1, samtools and htseq-count in my current pipeline but so far I have always limited my analysis to the annotated regions so I am not sure how to best approach this - any advice is most welcome!

[syfo]

If you have read alignments in bam format on the one hand and a gene annotation in gtf/gff/bed format on the other hand it should be simple to compare both with bedtools.

For instance something like "bedtools intersect -v -abam mapped_reads.bam -b exons.bed" should give you the reads that do not intersect any annotated exon (you can use entire gene regions if you prefer to discard intronic reads). Pipe this to "samtools view -c -" to get the number of reads.

More precisely:

bam1=mapped_reads.bam
bam2=intergenic_reads.bam
annotation=annotated_genes.bed

bedtools intersect -v -abam $bam1 -b $annotation > $bam2
samtools index $bam2

And then either something like

for chr in `cat list_of_chromosomes.txt` ; do
 N=`samtools view -c $bam1 $chr` ;
 n=` samtools view -c $bam2 $chr` ;
 echo $chr $n $N ;
done

Or simply

samtools idxstats $bam2

You might also be interested in how close your reads are from annotated exons. Have a look at "bedtools closest".

[kokonech]

Qualimap (http://qualimap.org) outputs percentage of reads in intronic and intergenic regions as one of the quality metrics for RNA-seq data.
Example report can be found here.