extracting data from bam file

Question

[Sara]

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

[poojasethiya]

You can also provide bed file and get the bam output

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

~Pooja

[cpad0112]

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

Answer 1

[Philipp Bayer]

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.