extracting reads from sam/bam file gene wise

[chris_bioinfo]

Dear All,

I have aligned Illumina transcriptome data ( paired-end 100bp) onto genes using bowtie2, which consists of ~10,000 genes.

Now I want to know for EACH GENE,

1) how many reads are aligned concordantly
2) how many reads are aligned discordantly
3) how many reads are aligned as single-end reads

Any help is appreciated.

Thank you,
Christopher

[sivasubramani]

Hope you got the answer..

[kushald]

I had the same question Chris. Please let me know if you get know the answers. Cheers mate.

[dpryan]

There are a number of ways to do this, depending on your comfort in coding. One way would be to write a shell script to read in the SAM/BAM header with samtools to parse out the gene list and then use "samtools view -f XX aligned.bam gene" on each of those genes with an appropriate flag option for each of your cases.

A better method would be to just write a script in python/C/R/whatever to read through the SAM/BAM file once and tally accordingly as it goes. That's pretty straight-forward to do if you can program (if not, you'll need to learn anyway).

[chris_bioinfo]

Thanks guys, I wrote a perl program for this and it's working fine

Best,
Christopher

[rajamoturi]

Dear Chris

Please share the code if you can it will be very helpful to us


cheers
Raja