Bedtools Genomecoveragebed Usage : How To do per base coverage analysis for only chr21?

Question

[gskbioinfo143]

Hey any one pls explain how to do per base coverage analysis for only particular chr from given bam file using bedtools or GATK.

Answer 1

[Irsan]

As a start, you need a file with the start and end coordinates of your chromosomes. Then prepare a bed-file for your given chr, e.g. chr21:

bedtools makewindows -g hg19-chrLengths.txt -w 1 | grep '^chr21\t' > chr21-bases.bed

Calculate coverage for each base

bedtools coverage -abam yourFile.bam -b chr21-bases.bed

You can circumvent the grep-part if you make sure that in the hg19-chrLengths.txt file only chr21 is there

[gskbioinfo143]

Hey thanks for the reply but the first command is taking time..i have given genome file which bedtools provided. need help.

Can you please explain how to do the same with GATK.

[Irsan]

Actually the only thing that this genome file needs to contain is:

chr1    249250621
chr10   135534747
chr11   135006516
chr12   133851895
chr13   115169878
chr14   107349540
chr15   102531392
chr16   90354753
chr17   81195210
chr18   78077248
chr19   59128983
chr2    243199373
chr20   63025520
chr21   48129895
chr22   51304566
chr3    198022430
chr4    191154276
chr5    180915260
chr6    171115067
chr7    159138663
chr8    146364022
chr9    141213431
chrX    155270560
chrY    59373566

Optionally you can remove the chrs for which you have no interest, can you try that?

[gskbioinfo143]

pls mention how much time it will take for first command..

[Irsan]

I did it only for chr21 produced an output file of 1.1 Gb and it took 3 minutes on 1 cpu at 2.67 Ghz. The more chromosomes you include and the longer they are, the longer it will take. Note: when you do this for the complete human genome, (1 row per base in the output) you get a file of 3 billion rows!

[gskbioinfo143]

hey i have extracted chr21 from bam using samtools it took 1min then did per base coverage with bedtools coverage

[Irsan]

So it worked? Nice job!

[gskbioinfo143]

yes it worked thanks

[gskbioinfo143]

how to find the no of reads matched in chr21 region? any visualization blocks?