Bam file to junctions.bed

[adrian]

Hi:

I received aligned BAM file and do not have a raw sequence file.

My aim is to count how many reads skip exon 7 of a gene and how many reads do not skip, in addition to reads that span exons 6 and 7 ; 7 and 8.

Ex 6-------- Ex 7 -------- Ex 8
___________ __________ => Condition1: reads that span 6-7 and exs 7-8
____//////////////////////___ => Condition 2:reads skipping exon 7
___________________ => Condition 3: reads that span exon 6,7 and end in 8.



Is there a way to get numbers from my BAM file for above 3 conditions.

Yes, if I have raw reads, i would use TopHat and get junctions.bed to deduce these. However, the BAM file was not generated using TopHat and I don't have access to raw sequences.

Is there a way to get junctions.bed - perhaps convert bam to FASTA and then realign using Tophat. This potentially 'would' corrupt the paired end structure..leading to loss of read numbers..( I am not so sure about this though)

Or

Is there any other smart way to just count reads that jump exon 7.

Appreciate any response. Thanks a lot.

Adrian

[dpryan]

The simplest method would be to just script this in pysam (or whatever language you prefer).

BTW, you can convert the BAM file to fastq and realign that, but it's faster to just write a little python script.

[adrian]

Thank you.
Is there a particular function that I could use? If not what would be the logic to get those read stats.

thanks
Adrian

[dpryan]

The general idea is to:

1. Iterate over the reads
2. For each read, get its start and end position.
3. If at least one of the exons could be between those coordinates then get the CIGAR
4. Parse the CIGAR string into a sequence of aligned regions
5. For each region, note if it overlaps one of your exons. Add that to a vector or a data structure of your choice (you could even just use an integer as a bitmap).
6. Once you've iterated through the aligned regions for a read of interest, look at the structure from the previous step and proceed as desired.

That's the general idea. If your BAM file is coordinate sorted and indexed, then you can simply request the reads covering the regions of interest, which will make things a bit quicker.

[gringer]

Here's a rough idea of how to do bam2fastq:

samtools view file.bam | awk -F '\t' '{print ">"$1"\n"$10"\n+\n"$11}' > file.fastq

Unfortunately this will give you a fastq file with interleaved reads, which can be a little bit of a pain to use. You can use the filter function (-f / -F) of samtools view to get around that, reading through the BAM file twice:

samtools view -f 0x40 file.bam | awk -F '\t' '{print ">"$1"\n"$10"\n+\n"$11}' > file_R1.fastq
samtools view -f 0x80 file.bam | awk -F '\t' '{print ">"$1"\n"$10"\n+\n"$11}' > file_R2.fastq

The SAM File format specification is your friend, see section 1.4.

The process of BAM -> FASTQ -> Tophat is slower in terms of computer time, but from your description it sounds like it will be quicker in terms of bum-on-seat time.