Bio Bam Statistics

input_read_count = flagstat_total - secondary - supplementary
                 = stats_raw_total_sequences

samtools flagstat input.bam

10000000 + 0 in total (QC-passed reads + QC-failed reads)
9950000 + 0 primary
0 + 0 secondary
50000 + 0 supplementary
0 + 0 duplicates
0 + 0 primary duplicates
9800000 + 0 mapped (98.00% : N/A)
9750000 + 0 primary mapped (97.99% : N/A)
9950000 + 0 paired in sequencing
4975000 + 0 read1
4975000 + 0 read2
9700000 + 0 properly paired (97.49% : N/A)
9720000 + 0 with itself and mate mapped
30000 + 0 singletons (0.30% : N/A)
15000 + 0 with mate mapped to a different chr
10000 + 0 with mate mapped to a different chr (mapQ>=5)

samtools flagstat -@ 4 input.bam

samtools flagstat input.bam > flagstat.txt

samtools idxstats input.bam

chr1    248956422    5000000    1000
chr2    242193529    4800000    800
chrM    16569        50000      100
*       0            0          150000

# Total mapped reads
samtools idxstats input.bam | awk '{sum += $3} END {print sum}'

# Mitochondrial percentage
samtools idxstats input.bam | awk '
    /^chrM/ {mt = $3}
    {total += $3}
    END {print mt/total*100 "% mitochondrial"}'

samtools stats input.bam > stats.txt

samtools stats input.bam | grep "^SN"

samtools stats input.bam > stats.txt
plot-bamstats -p plots/ stats.txt

samtools stats input.bam chr1:1000000-2000000 > region_stats.txt

samtools depth input.bam > depth.txt

samtools depth -r chr1:1000-2000 input.bam

samtools depth -a input.bam > depth_with_zeros.txt

# samtools mpileup historically capped depth at 8000 per position -- the cap was in mpileup, not depth.
# samtools depth -d/--max-depth is deprecated in 1.13+ (silently ignored).
# For mpileup, raise the cap explicitly when working with deep targeted/amplicon data:
samtools mpileup -d 1000000 -f ref.fa input.bam

# When fragment length < 2 * read_length, R1 and R2 overlap.
# Default samtools depth double-counts overlap; -s deducts:
samtools depth -s input.bam

mosdepth -t 4 sample input.bam                                                       # genome-wide per-base
mosdepth -t 4 --by exome.bed --thresholds 1,10,20,30,100 --no-per-base sample input.bam   # exome QC
mosdepth -t 4 --quantize 0:1:10:100: sample input.bam                                # CNV-style bands
mosdepth -t 4 -f ref.fa sample input.cram                                            # CRAM with reference

samtools depth -b regions.bed input.bam

samtools depth input.bam | awk '{sum += $3; n++} END {print sum/n}'

samtools coverage input.bam

samtools coverage -r chr1:1000000-2000000 input.bam

samtools coverage -b regions.bed input.bam

samtools coverage -m input.bam

import pysam

with pysam.AlignmentFile('input.bam', 'rb') as bam:
    total = mapped = paired = proper = 0
    for read in bam:
        total += 1
        if not read.is_unmapped:
            mapped += 1
        if read.is_paired:
            paired += 1
        if read.is_proper_pair:
            proper += 1

    print(f'Total: {total}')
    print(f'Mapped: {mapped} ({mapped/total*100:.1f}%)')
    print(f'Properly paired: {proper} ({proper/paired*100:.1f}%)')

import pysam

with pysam.AlignmentFile('input.bam', 'rb') as bam:
    for stat in bam.get_index_statistics():
        print(f'{stat.contig}: {stat.mapped} mapped, {stat.unmapped} unmapped')

import pysam

with pysam.AlignmentFile('input.bam', 'rb') as bam:
    for pileup in bam.pileup('chr1', 1000000, 1000001):
        print(f'Position {pileup.pos}: depth {pileup.n}')

import pysam

def mean_depth(bam_path, chrom, start, end):
    depths = []
    with pysam.AlignmentFile(bam_path, 'rb') as bam:
        for pileup in bam.pileup(chrom, start, end, truncate=True):
            depths.append(pileup.n)

    if depths:
        return sum(depths) / len(depths)
    return 0

depth = mean_depth('input.bam', 'chr1', 1000000, 2000000)
print(f'Mean depth: {depth:.1f}x')

import pysam

def coverage_stats(bam_path, chrom, start, end):
    covered = 0
    total_depth = 0

    with pysam.AlignmentFile(bam_path, 'rb') as bam:
        for pileup in bam.pileup(chrom, start, end, truncate=True):
            covered += 1
            total_depth += pileup.n

    length = end - start
    pct_covered = covered / length * 100
    mean_depth = total_depth / length if length > 0 else 0

    return {
        'length': length,
        'covered_bases': covered,
        'pct_covered': pct_covered,
        'mean_depth': mean_depth
    }

stats = coverage_stats('input.bam', 'chr1', 1000000, 2000000)
print(f'Coverage: {stats["pct_covered"]:.1f}%')
print(f'Mean depth: {stats["mean_depth"]:.1f}x')