extracting reads from sam/bam file gene wise

chris_bioinfo · 07-24-2013, 05:58 AM

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 · 07-24-2013, 11:02 PM

Hope you got the answer..

kushald · 07-24-2013, 11:06 PM

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

dpryan · 07-25-2013, 12:10 AM

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 · 07-25-2013, 01:54 AM

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

Best,
Christopher

rajamoturi · 09-20-2013, 03:09 AM

Dear Chris

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

cheers
Raja

Similar Threads
Thread	Thread Starter	Forum	Replies	Last Post
Extracting unpaired reads from BAM file	JChase	Bioinformatics	4	08-13-2012 07:23 PM
Filter paired & mapped reads from SAM/BAM file	pirates.genome	Bioinformatics	2	06-19-2012 02:21 AM
extracting certain scaffolds from SAM file	lukas1848	Bioinformatics	1	01-18-2012 09:33 AM
Extract perfectly mapped reads from SAM/BAM file	Graham Etherington	Bioinformatics	2	07-21-2011 08:27 AM
Extracting unique reads from a .ma or .bam file?	JohnK	SOLiD	14	06-04-2010 01:32 AM

07-24-2013, 05:58 AM	#1
chris_bioinfo Junior Member Location: London Join Date: Aug 2012 Posts: 8	extracting reads from sam/bam file gene wise 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

07-24-2013, 11:02 PM	#2
sivasubramani Member Location: India Join Date: Apr 2011 Posts: 14	Hope you got the answer..

07-24-2013, 11:06 PM	#3
kushald Member Location: Bangalore, India Join Date: Mar 2013 Posts: 42	I had the same question Chris. Please let me know if you get know the answers. Cheers mate.

07-25-2013, 12:10 AM	#4
dpryan Devon Ryan Location: Freiburg, Germany Join Date: Jul 2011 Posts: 3,480	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).

07-25-2013, 01:54 AM	#5
chris_bioinfo Junior Member Location: London Join Date: Aug 2012 Posts: 8	Thanks guys, I wrote a perl program for this and it's working fine Best, Christopher

09-20-2013, 03:09 AM	#6
rajamoturi Junior Member Location: leeds Join Date: May 2011 Posts: 4	Dear Chris Please share the code if you can it will be very helpful to us cheers Raja