This is a test version of Biostars. For the public version, visit https://www.biostars.org.
Extract the sequence of a gene from an aligned BAM file

I would like to calculate the percentage of mapping of a particular gene to the genome. Say for eg, am looking at gene A

I extracted the regions of that gene from the genome file and created a bed file

 Homo_sapiens.GRCh38_geneA.bed
 6       52186375        52190638

Then I tried to get the corresponding regions of geneA from the BAM file

samtools view -h -b -L  Homo_sapiens.GRCh38_geneA.bed  SampleAligned.sortedByCoord.out.bam >geneA.bam

I tried 

 samtools flagstat geneA.bam

I do not see much information from the output.

How do I go about getting the percentage of reads from gene A that mapped to the genome? Is this the right way of doing this?

bam samtools bedtools
• 195 views • link
updated by swbarnes2 swbarnes2 • written by rykerklie7 • rykerklie7

Could you please show the first several reads of your original bam file and the results of samtools flagstat?

• 0 views • link
written by Jianyu Jianyu

1 answer

You have the # of reads in the region, what more do you want? Can't you do the division yourself?

• 0 views • link
written by swbarnes2 swbarnes2

Did not know to interpret the % of reads for geneA from my sample mapping to the genome

• 0 views • link
written by rykerklie7 • rykerklie7

Log in to answer this question.

More posts like this