slice Bam

[Lilach1234]

I have a bam file and i wish to exclude all the reads that are part of known genes. i got the exon list + Coordinates from ucsc as a .gtf file, but i'm not sure what's next... should I use samtools? if so, what options?

Thank you very much in advance,
Lilach

[dpryan]

You might use htseq-count and use its option to write to SAM format. It then adds a flag indicating if there was a gene assignment or not. You can then simply use grep to extract the "no_feature" (or whatever that's called) reads. There are faster ways, but this is likely the simplest.

[dariober]

I have a bam file and i wish to exclude all the reads that are part of known genes. i got the exon list + Coordinates from ucsc as a .gtf file, but i'm not sure what's next... should I use samtools? if so, what options?

Thank you very much in advance,
Lilach

This might be a solution (not tested) using bedtools and samtools.

First, create a bed file of intergenic regions, then use samtools to extract reads in the bed intervals.

bedtools complement -i genes.gtf -g genome.txt > intergenic_regions.bed
samtools view aln.bam -L intergenic_regions.bed > intergenic_reads.sam

genome.txt is a tab-separated file of chromosome sizes which you can get from UCSC. See help in bedtools complement.

Hope this helps
Dario