How to extract only intronic reads or region from BAM file?

[ClemBuntu]

Hello all,

I used picard tool CollectRnaSeqMetrics on my RNA-seq fastq files and I found a lot (> 50%) of intronic region.

I wonder if these intronic regions are equally distributed among the genome or otherwise if my reads are aligned one only few genes.

How can I see it ? Is there a way to extract only intronic regions from a BAM file ?

Or maybe there is a way to extract only reads which mapped on intronic regions and see where they're aligned ?

[GenoMax]

Look at the samtools view -L option. Provide your intronic regions as a bed file and extract reads aligned in those areas into a new file.

[WhatsOEver]

First, depending on the organism you are working with, it's not totally surprising that a lot of reads map to introns.
Second, if you have a lot of genes with different isoforms in your sample, the number of "intronic reads" may be influenced by the annotation file you are using.

There is at least one easy way I can think of to get non-exon regions from a BAM file. bedtools intersect (http://bedtools.readthedocs.org/en/l...intersect.html) allows you to select for reads that don't overlap defined regions (-v parameter). But, you will eventually get the inter-gene reads as well... So you would actually need to define gene regions first (maybe from 5'UTR start to 3'UTR end) and then use this file together with an exon defining file in bedtools intersect -v. This will give you an intron-only file which you can use for detecting unequal mapping distributions.

[ClemBuntu]

Hello,

Thank you guys for your answers, all I have is a GTF file with the transcripts coordinates (3th column contains only "CDS" "exon" "start_codon" and "stop_codon") so bedtools intersect with -v option seems more suitable for me.

Do you think a grep "*_codon" on my transcripts file will be do the trick to have only the gene regions ?

So if I well understood I should do something like that :
bedtools intersect -abam -a myfile.bam -b onlyGeneCoordinate.gtf > bamVSgene.bed #first I want the intersection between my reads and genes

bedtools intersect -v -abam -a bamVSgene.bed -b genesAndTranscripts.gtf > intronOnlyFile.bed #filter the exonic regions

bedtools intersect -abam -a myfile.bam -b intronOnlyFile.bed > onlyIntrons.bed #finally get only the intronic reads

I'm not sure the 2nd step will work because there is not intronic regions in the "genesAndTranscripts.gtf" file.

[minie]

Hi how did you get the intron.bed file? I have gff/gtf file with transcripts/exon/gene/cds information snd trying to parse intronic region. any help ?