Softwares For Calculation Of Genome-Wide Alignments Coverage

Question

[Bioscientist]

Any recommendation for packages to calculate alignment coverage from bam/sam files? 1. I need whole genome-wide calculation instead of certain genomic regions. 2. I mean available softwares instead of playing with bash code myself, :)

Bedtools, GATK, Picardtools, Samtools....I search a little bit but cannot find what I want.

thx

[Sean Davis]

Does samtools mpileup not give you what you need?

[cobalym7]

I think Qualimap can help you. It is easy to use and give you breadth of coverage in each deep of coverage.

Answer 1

[Aaronquinlan]

The bedtools command bedtools genomecov (aka genomeCoverageBed) will do this for you. (The commands below are for version 2.15.0 and above).

By default, you will get a histogram of coverage for each chromosome, and the last group of results will be a histogram of coverage for the genome:

bedtools genomecov -ibam aln.sorted.bam | grep genome

It will also report genome-wide coverage in BEDGRAPH format (only covered intervals):

bedtools genomecov -ibam aln.sorted.bam -bg

BEDGRAPH format with both covered and uncovered intervals:

bedtools genomecov -ibam aln.sorted.bam -bga

Only coverage on the positive strand:

bedtools genomecov -ibam aln.sorted.bam -s +

Scaling coverage by a constant for things like RPM:

bedtools genomecov -ibam aln.sorted.bam -scale 0.73

Taking into account spliced BAM alignments:

bedtools genomecov -ibam aln.sorted.bam -split

Only counting MAPQ >= 30:

samtools view -q 30 -u aln.bam | bedtools genomecov -ibam -

If you want per-base coverage at each base in the genome, you can use the -d option:

bedtools genomecov -ibam aln.sorted.bam -d

samtools will also do this much more quickly, but I believe it will only report positions with coverage > 0:

samtools depth aln.bam

[Madelaine Gogol]

Ooh, fancy new way of calling bedtools.

[Bioscientist]

Oh thanks so much! Is this a novel function? I checked the manual but didn't find that. (http://code.google.com/p/bedtools/wiki/Usage#coverageBed) Probably you need to upgrade the site. :)

thx a lot! Bedtools is super~

[Sean Davis]

I had missed genomecov. Thanks for the heads-up.

[Aaronquinlan]

It is also known as genomeCoverageBed and is in the PDF manual. That said, you are right that the docs are dated. Sorry. Manual follows http://code.google.com/p/bedtools/downloads/detail?name=BEDTools-User-Manual.v4.pdf

[Aaronquinlan]

Note the input param should be -ibam. I erroneously wrote -I at first.

[Bioscientist]

I find it requires a genome file, with chromsize. How can I get that? thx....

[moorem]

So can bedtools be used to determine the 'alignment length' of short reads? For example if I wanted to calculate percent identity of a region that has large deletions where there's no coverage and then get mismatches/small indels from GATK to get percent identity overall?

In a similar way to do so from BLAST