I would like to get the average coverage of all the captured bases in a bam file. What would be the best way to do this? What I am looking is a simple one number like 40X. Given that there may be millions of bases sequenced in a next gen study I would like to get the overall average coverage for all these bases as a single number.
Thanks
if you have samtools and awk, you can keep it simple with this:
samtools depth *bamfile* | awk '{sum+=$3} END { print "Average = ",sum/NR}'
and with stdev:
samtools depth *bamfile* | awk '{sum+=$3; sumsq+=$3*$3} END { print "Average = ",sum/NR; print "Stdev = ",sqrt(sumsq/NR - (sum/NR)**2)}'
Not sure this works. It calculates coverage over regions but does not account for regions that were not covered with a read. Then it divides all the sum with the number of rows that were reported... but the positions without a read hit are not included. Basically - it gives a too high result.
That's true. To be fair, the denominator needs to be the size of genome where the BAM file was generated against, rather than the number of bases covered at least once (that's what NR is in this context). To get the genome size from the BAM file:
samtools view -H *bamfile* | grep -P '^@SQ' | cut -f 3 -d ':' | awk '{sum+=$1} END {print sum}'
Try the genomeCoverageBed tool in the BEDtools package, which takes a BAM or BED file as input and:
computes a histogram of feature coverage (e.g., aligned sequences) for a given genome. Optionally, by using the –d option, it will report the depth of coverage at each base on each chromosome in the genome file (-g).
The way we run this is to first make a tab-delimited genome.txt file with "contig_name contig_size" information for all contigs in your genome and then run genomeCoverageBed on a sorted BAM file:
$ genomeCoverageBed -ibam sortedBamFile.bam -g genome.txt > coverage.txt
This will give you a summary histogram of coverage across each contig and for the entire genome, from which you can obtain the modal or mean value (the single fold-coverage value you are looking for).
That's because it reports a histogram. From the documentation: The values in the output are:
[all values taken from your first row]
If you want to get the average coverage: add up the product of bases per coverage [2* 4+3* 28+4* 10+...] and divide by the total number of bases [137928].
And the mandatory R solution. The Rsamtools package can be used to read in the BAM file.
Here is a code example:
library(Rsamtools)
bamcoverage <- function (bamfile) {
# read in the bam file
bam <- scanBam(bamfile)[[1]] # the result comes in nested lists
# filter reads without match position
ind <- ! is.na(bam$pos)
## remove non-matches, they are not relevant to us
bam <- lapply(bam, function(x) x[ind])
ranges <- IRanges(start=bam$pos, width=bam$qwidth, names=make.names(bam$qname, unique=TRUE))
## names of the bam data frame:
## "qname" "flag" "rname" "strand" "pos" "qwidth"
## "mapq" "cigar" "mrnm" "mpos" "isize" "seq" "qual"
## construc: genomic ranges object containing all reads
ranges <- GRanges(seqnames=Rle(bam$rname), ranges=ranges, strand=Rle(bam$strand), flag=bam$flag, readid=bam$rname )
## returns a coverage for each reference sequence (aka. chromosome) in the bam file
return (mean(coverage(ranges)))
}
And the output:
> bamcoverage("~/test.bam")
gi|225184640|emb|AL009126.3|
15.35720
The good thing about this is you can do much more with the genomic ranges object than computing the coverage, the downside is, that the code is not that efficient.
did you try the GATK "DepthOfCoverage" ? http://www.broadinstitute.org/gsa/wiki/index.php/Depth_of_Coverage_v3.0
or you can run 'samtools pileup' and calculate the mean value for the coverage column.
I was interested in this issue too, but I wanted to know the coverage for certain genes. after looking to this tool's wiki I've realized that I can feed it with a list of genes intervals (which of course I'll have to retrieve from UCSC or Ensembl) and I'll surely sort my problem out! look what an unexpected great piece of information I've just found! once again thanks a lot, Pierre.
Here's a very simple perl script I wrote that does exactly this that you're free to use:
($num,$den)=(0,0);
while ($cov=<STDIN>) {
$num=$num+$cov;
$den++;
}
$cov=$num/$den;
print "Mean Coverage = $cov\n";
Just pipe samtools pileup to it like so:
/path/to/samtools pileup in.bam | awk '{print $4}' | perl mean_coverage.pl
It will also work easily for genomic regions. Just index the BAM file and use samtools view like so:
/path/to/samtools view -b in.bam <genomic region> | /path/to/samtools pileup - | awk '{print $4}' | perl mean_coverage.pl
Hope someone finds this useful.
One more comment on this pipeline... In my data, where I have relatively low coverage there are quite a few genomic positions not covered at all. These will not show up in the pileup command, so I guess that the mean coverage calculated this way will be the mean coverage of the covered positions. Am I right? Is there a way (in the "samtools pileup" command) to output all positions even if they are not covered by any read?
samtools depth bamfile | awk '{sum+=$3} END { print "Average = ",sum/NR}'
by @russdurrett also gives same output for my paired-end bam file from RNA-Seq:-)
Hello, I know this post is quite old, but I am trying to use this as it seems to work wonderfully for everyone else. I am new to bioinformatics and genomics, so please be patient with me. I have run the command 'samtools mpileup A186944_mtgenome.srt.dup.bam | awk '{print $4}' | perl coverage.pl' where my bam file is a sorted bam with duplicates removed and I get the following error:
Illegal division by zero at coverage.pl line 6.
Can someone please help me figure out the error? I think it's maybe because there are 0's where unmapped reads occur? How do I obtain coverage for just the mapped reads, and eliminate the 0 denominator?
A modern, very fast, solution for this would be mosdepth.
You can install the executable from their github release page: https://github.com/brentp/mosdepth/releases
Latest version at time of commenting (0.2.9):
(Download executable)
wget https://github.com/brentp/mosdepth/releases/download/v0.2.9/mosdepth
(Update permissions to run it)
chmod a+x mosdepth
(Good to go!)
<path_to_mosdepth> --help
Why not use bioconda?
More installation options and instructions can be found on the GitHub page.
@michael.james.clark, thank you. I wrote for myself a script that I can add to when the need arises:
#!/usr/bin/env perl
use strict;
use warnings;
# Usage: samtools pileup file.bam | bam_coverage.pl
my $num; # per residue coverage
my $len; # sequence length counter
my $min = 1000; # minimum coverage
my $max = 0; # maximum coverage
while (<>) {
my @a = split /\t/;
$num += $a[3];
$min = $a[3] if $min > $a[3];
$max = $a[3] if $max < $a[3];
$len++;
}
printf "Mean coverage : %.1f\n", $num/$len;
printf "Coverage range : %d - %d\n", $min, $max;
How about just running qualimap in GUI mode or to create a QC PDF? Your can run it on the whole genome or on 1 or multiple regions using a bed file.
use
samtools pileup accepted_hits.bam | awk '{ count++ ; SUM += $4 } END { print "Total: " SUM "\t" "Nucleotides: " count "\t" "Average_coverage: " SUM/count }'
(use $8 instead; if pileup with -cf option)
I used the above command to extract the average coverage on my aligned bam files and my statistics are: Total: 8011581243 Nucleotides: 510234050 Average_coverage: 15.7018
I want to know that the average_coverage is stating that on a distribution each of my base is read atleast 15 time right? What does the total signifies and the Nucleotides count: 510234050. Is this the total no of bases that are in the exonic region? then how will I calculate the number of reads from here that actually got mapped on in the exome region? Can we do this from the nucleotides count? I know for each reads consist of 100 bases so if I divide this 510234050 bases with 100 it gives me the reads but that will not give me the reads that actually got mapped in the exome region. How t extract that?
This is for people who have a sam/bam file consisting of reads mapped to multiple contigs (samtools needs to be in $PATH):
BBMap: http://sourceforge.net/projects/bbmap/
bbmap/pileup.sh -h
Written by Brian Bushnell
Last modified May 19, 2015
Description: Calculates per-scaffold coverage information from an unsorted sam file.
Usage: pileup.sh in=<input> out=<output>
Input may be stdin or a SAM file, compressed or uncompressed.
Output may be stdout or a file.
Input Parameters:
in=<file> The input sam file; this is the only required parameter.
ref=<file> Scans a reference fasta for per-scaffold GC counts, not otherwise needed.
fastaorf=<file> An optional fasta file with ORF header information in PRODIGAL's output format. Must also specify 'outorf'.
unpigz=t Decompress with pigz for faster decompression.
..
Use of this site constitutes acceptance of our User Agreement and Privacy Policy.
version 2.3.6
People seem to be equally devided between GATK, Samtools and BedTools.
The Depth of Coverage is depreciated in GATK4. It's only available in GATK3. See here
Most answers seems to be very old and hence would like to have updated suggestions.
I have 6 bam files and I have used
samtools depthto calculate chromosome wise depth for all chromosomes and then plotted in R. While looking at the plots at 2-3 places, depth shows upto 200-3500 and hence I would like to calculate average read depth of each chromosome from each bam file. For e.g.: 1) average depth of chr1 from bam1, average depth of chr1 from bam2, ...until bam6. 2) average depth of chr1 from all 6 bam together.Kindly guide.
Thanks.
Try the solutions if they don't work then you need to reframe the question with your commands and write a new query in a new post. This is advised. And yes paired end should work.
in addition, if you have already plotted depths it means that you have got a temporal R data.frame (or similar type) in there containing chromosome, position and depth information. splitting such data.frame by chromosome and calculating average values should be straight-forward. I would just suggest to take into account empty values not reported by samtools depth dividing each chromosome depths' sum by the size of the chromosome, and not by the number of depth values.
Why is it wrong to take the number of sequenced bases and divide by the genome size?
because you cannot align all the bases on the genome, some reads are optical duplicates, etc...
Okay, contamination is a fair point, but even if some reads fail to match, that could still be a sign of genetic divergence, and should still count towards genome coverage. Do any of the other suggestions deal with duplicates implicitly?
Qualimap is very fast and excellent tool. thanks
Hi,
I have multiple paired-end bam files from RNA-Seq data, already aligned and computed depth with
samtools depth > bam.depth. Would any one please guide for any straight forward/latest way to compute coverage plots, avg. coverage plots and other metrics from depth files?Picard CollectRNAseqmetrics,RSeQCandQoRTshave not been helpful. :-(Thanks.
Have you run any of the other tools mentioned in this thread (e.g. Qualimap, pileup from BBMap)?
@genomax My bam files are very big and
Qualimaphas not been helpful, so asbbmapwhich is Java-based (memory issues likePicardandQoRTs) and hence I usedsamtools depth. My depth files are also in GB.Both programs will happily accept larger allocations of RAM and can open very large files. If you don't have access to a machine with more RAM then that is a different problem.