def get_total_coverage_estimate(bamfile, nhead=1000000):
    """
    Estimate the total coverage across all bams and
    positions using the number of mapped reads and
    mean read length.
    """
    num_mapped_reads = bamfile.mapped
    mean_read_length = sum(map(lambda r: r.qlen, bamfile.head(nhead)))/nhead
    return num_mapped_reads * mean_read_length