This is a test version of Biostars. For the public version, visit https://www.biostars.org.
extracting data from bam file

I have bam file (result of alignment of RNAseq to the genome). now I would like to get the coverage of a specific region for instance chr3: 123545-123965. do you know how I can do that?

here is the bam file structure:

@HD VN:1.0  SO:coordinate
@SQ SN:chr1 LN:249250621
@SQ SN:chr10    LN:135534747
@SQ SN:chr11    LN:135006516
@SQ SN:chr12    LN:133851895
@SQ SN:chr13    LN:115169878
@SQ SN:chr14    LN:107349540
@SQ SN:chr15    LN:102531392
@SQ SN:chr16    LN:90354753
@SQ SN:chr17    LN:81195210
@SQ SN:chr18    LN:78077248
@SQ SN:chr19    LN:59128983
@SQ SN:chr2 LN:243199373
@SQ SN:chr20    LN:63025520
@SQ SN:chr21    LN:48129895
@SQ SN:chr22    LN:51304566
@SQ SN:chr3 LN:198022430
@SQ SN:chr4 LN:191154276
@SQ SN:chr5 LN:180915260
@SQ SN:chr6 LN:171115067
@SQ SN:chr7 LN:159138663
@SQ SN:chr8 LN:146364022
@SQ SN:chr9 LN:141213431
@SQ SN:chrM LN:16569
@SQ SN:chrX LN:155270560
@SQ SN:chrY LN:59373566
rna-seq
• 3,131 views • link
updated by Philipp Bayer Philipp Bayer • written by Sara Sara

You can also provide bed file and get the bam output

samtools view -b -L myGenes.bed Input.bam > out.bam

~Pooja

• 0 views • link
written by poojasethiya poojasethiya

syntax:

java -jar igvtools.jar count --query <querycoordinates> -w <windowsize> <input bam> <output.wig> <genomeversion>

Example code for Overall coverage:

java -jar igvtools_2.3.95/igvtools.jar count --query chr3:123545-123965 -w 1  --bases input.bam output.wig hg19

Example code for Per base coverage:

java -jar igvtools_2.3.95/igvtools.jar count --query chr3:123545-123965 -w 1  --bases input.bam output.wig hg19
• 0 views • link
written by cpad0112 cpad0112

1 answer

You could use bedtools genomecov aka genomeCoverageBed, slightly adjusted from the manual:

samtools view -b <BAM> chr3:123545-123965 | genomeCoverageBed -ibam stdin

Add -bga to genomeCoverageBed if you want a different format which will also report bases with 0 coverage.

• 0 views • link
written by Philipp Bayer Philipp Bayer

Log in to answer this question.

More posts like this