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)}'

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).

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.

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.

A modern, very fast, solution for this would be mosdepth.

@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;

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)

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.

http://qualimap.bioinfo.cipf.es/

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.

..

To make this thread a little more complete. adding pandepth as a current tool : https://github.com/HuiyangYu/PanDepth

BioPerl has a library named Bio-Samtools which provides a familiar Perl/BioPerl access to samtools. I heard the author, Lincoln Stein, give a talk on on how he wove the two together use XS. As I recall, one of his examples was exactly the question you posed

Now the new release of Samtools for coverage analysis.

samtools coverage input.bam

Genozip calculates coverage a bit differently, based on the number of bases rather than reads - ignoring hard-clipped and soft-clipped bases, as well as unmapped, secondary, supplementary and duplicate alignments. It can also calculate approximate coverage directly on a FASTQ file.

 $ genocat --coverage my-sample.bam.genozip

--coverage for: my-sample.bam.genozip
Contig         LN        Reads        Bases       Bases   Depth
1              249.3 Mb  61,215,549   9.0 Gb      7.1 %   36.10
2              243.2 Mb  63,342,488   9.3 Gb      7.3 %   38.26
3              198.0 Mb  50,819,357   7.5 Gb      5.9 %   37.76
4              191.2 Mb  48,131,726   7.1 Gb      5.6 %   37.04
5              180.9 Mb  46,495,225   6.8 Gb      5.4 %   37.81
6              171.1 Mb  43,659,445   6.4 Gb      5.0 %   37.53
7              159.1 Mb  41,171,496   6.0 Gb      4.7 %   38.01
8              146.4 Mb  38,081,676   5.6 Gb      4.4 %   38.27
9              141.2 Mb  31,755,198   4.7 Gb      3.7 %   33.05
10             135.5 Mb  34,954,121   5.1 Gb      4.0 %   37.90
11             135.0 Mb  35,234,129   5.2 Gb      4.1 %   38.38
12             133.9 Mb  34,523,219   5.1 Gb      4.0 %   37.91
13             115.2 Mb  24,526,664   3.6 Gb      2.8 %   31.32
14             107.3 Mb  23,525,557   3.5 Gb      2.7 %   32.21
15             102.5 Mb  22,613,925   3.3 Gb      2.6 %   32.41
16             90.4 Mb   23,675,707   3.5 Gb      2.7 %   38.41
17             81.2 Mb   22,148,439   3.2 Gb      2.5 %   40.01
18             78.1 Mb   19,542,517   2.9 Gb      2.3 %   36.81
19             59.1 Mb   16,374,766   2.4 Gb      1.9 %   40.50
20             63.0 Mb   16,530,922   2.4 Gb      1.9 %   38.52
21             48.1 Mb   9,811,790    1.4 Gb      1.1 %   29.93
22             51.3 Mb   10,238,901   1.5 Gb      1.2 %   29.26
X              155.3 Mb  20,421,365   3.0 Gb      2.3 %   19.29
Y              59.4 Mb   3,971,007    582.0 Mb    0.5 %    9.80
MT             16.6 Kb   105,906      15.5 Mb     0.0 %  937.41
Other contigs            28,989,258   4.1 Gb      3.2 %
-----
All contigs              771,860,353  113.3 Gb    88.8%   36.60
Soft clip                             1.2 Gb      0.9 %
Unmapped                 17,071,958   2.4 Gb      1.9 %
Supplementary            1,028,655    52.2 Mb     0.0 %
Duplicate                71,068,615   10.6 Gb     8.3 %
TOTAL                    861,029,581  127.5 Gb    100.0%

Documentation: https://genozip.com/coverage.html

Installing: https://genozip.com/installing.html

Publication: https://www.researchgate.net/publication/349347156_Genozip_-_A_Universal_Extensible_Genomic_Data_Compressor

Also tiwih meandepth

Nice tool to get coverage stats: https://github.com/brentp/goleft/tree/master/covstats