Bio Bam Statistics

Version Compatibility

Reference examples tested with: pysam 0.22+, samtools 1.19+

Before using code patterns, verify installed versions match. If versions differ:

• Python: pip show then help(module.function) to check signatures
• CLI: --version then --help to confirm flags

If code throws ImportError, AttributeError, or TypeError, introspect the installed package and adapt the example to match the actual API rather than retrying.

BAM Statistics

"Get alignment statistics and coverage from my BAM file" -> Generate read counts, mapping rates, per-chromosome statistics, depth profiles, and coverage summaries.

• CLI: samtools flagstat, samtools stats, samtools depth, samtools coverage (samtools)
• Python: pysam.AlignmentFile with pileup() and get_index_statistics() (pysam)

Generate alignment statistics using samtools and pysam.

Quick Summary Commands

| Question | Best tool | Why | |----------|-----------|-----| | Quick read counts by FLAG category | samtools flagstat | Fast; counts secondary+supp in totals | | Per-chromosome counts | samtools idxstats | Fast (needs index); counts secondary+supp | | Insert size, MAPQ, error, GC | samtools stats -r ref.fa | Comprehensive; feeds MultiQC | | Per-position depth (small region) | samtools depth or pysam pileup | Slow on full genome | | Per-position depth (genome-wide) | mosdepth | 3-10x faster than samtools depth | | Per-region coverage (BED) | mosdepth --by regions.bed | Production default | | Coverage histogram / cumulative | mosdepth -t 4 --no-per-base | Single-pass histogram | | Breadth at depth thresholds | mosdepth --thresholds 1,10,30,100 | Standard exome QC | | Targeted enrichment QC | picard CollectHsMetrics | PCTOFFBAIT, FOLD80BASE_PENALTY, AT/GC dropout | | Cross-sample contamination | verifybamid2, somalier | FREEMIX =5)

### Multi-threaded
```bash
samtools flagstat -@ 4 input.bam

Output to File

samtools flagstat input.bam > flagstat.txt

samtools idxstats

Per-chromosome read counts (requires index).

samtools idxstats input.bam

Output format: chrom length mapped unmapped

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

Parse idxstats

# 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

Comprehensive statistics including insert size, base quality, and more.

samtools stats input.bam > stats.txt

View Summary Numbers

samtools stats input.bam | grep "^SN"

Key summary fields:

• raw total sequences - Total reads
• reads mapped - Mapped reads
• reads mapped and paired - Properly paired
• insert size average - Mean insert size
• insert size standard deviation - Insert size spread
• average length - Mean read length
• error rate - Mismatch rate

Generate Plots (with plot-bamstats)

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

Stats for Specific Region

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

samtools depth

Per-position read depth.

Basic Depth

samtools depth input.bam > depth.txt

Output: chrom position depth

Depth at Specific Positions

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

Include Zero-Depth Positions

samtools depth -a input.bam > depth_with_zeros.txt

Maximum Depth Cap (Critical Trap)

# 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

Pipelines that historically break the 8000 mpileup cap: targeted oncology hotspots (5000-50000x), mitochondrial DNA (small genome, large read share), amplicon viral (ARTIC: 1000-100000x per amplicon), UMI-deduped capture (14000-17000x post-collapse), highly expressed transcripts (rRNA, mt-RNA).

Overlapping Pair Correction

# When fragment 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')

Insert Size Distribution

Goal: Compute the insert size distribution to assess library preparation quality.

Approach: Iterate properly paired read1 records, accumulate template lengths into a Counter, then compute summary statistics.

Reference (pysam 0.22+):

import pysam
from collections import Counter

insert_sizes = Counter()

with pysam.AlignmentFile('input.bam', 'rb') as bam:
    for read in bam:
        if read.is_proper_pair and read.is_read1 and read.template_length > 0:
            insert_sizes[read.template_length] += 1

sizes = list(insert_sizes.keys())
mean_insert = sum(s * c for s, c in insert_sizes.items()) / sum(insert_sizes.values())
print(f'Mean insert size: {mean_insert:.0f}')
print(f'Min: {min(sizes)}, Max: {max(sizes)}')

Quick Reference

| Task | Command | |------|---------| | Quick counts | samtools flagstat input.bam | | Per-chrom counts | samtools idxstats input.bam | | Full stats | samtools stats input.bam | | Coverage summary | samtools coverage input.bam | | Per-position depth | samtools depth input.bam | | Mean depth | samtools depth input.bam \| awk '{sum+=$3;n++}END{print sum/n}' |

QC Thresholds Are Assay-Specific

A single "mapping rate > 95%" rule rejects valid ATAC, ChIP, RNA-seq, metagenomics, and aDNA samples. The threshold question is "is this rate normal for this assay?" not "is this rate above 95%?"

| Metric | WGS PCR-free | WGS PCR | WES | Targeted panel | Deep panel (UMI) | RNA-seq | scRNA (10x) | ATAC | ChIP | Long-read | aDNA | |--------|--------------|---------|-----|----------------|------------------|---------|-------------|------|------|-----------|------| | Mapping rate | >99% | >98% | >95% | >95% | >95% | >90% | >70% | >50% | >60% | >95% | 1-50% | | Duplicate rate | 95% | >95% | >85% | >80% | >80% | >70% | n/a | >50% | >70% | n/a | >60% | | Mean MAPQ | bimodal at 0/60 | bimodal | bimodal | bimodal | bimodal | bimodal incl 255 (STAR) | 0/3/255 | 30-55 | 30-55 | 30-50 | 20-40 | | Mt fraction | 0.1-2% | 0.1-2% | = 30 is more informative:

samtools view -c -F 2308 -q 30 in.bam   # primary, mapped, MAPQ>=30
samtools view -c -F 2308 in.bam          # primary, mapped (denominator)
# For STAR/STARsolo, use -q 255 instead of -q 30 (255 is the unique-mapping sentinel)

What Flagstat Does Not Reveal

A 99% flagstat mapping rate does NOT mean the data is usable. Common false-positive scenarios:

1. Adapter readthrough: short fragments (insert 5% suggests adapter contamination

```

1. Off-target enrichment (capture/WES): detect via picard CollectHsMetrics PCTOFFBAIT or PCTSELECTEDBASES.
2. Low-complexity pile-up: telomere/centromere reads mass at MAPQ-0; counted as mapped but useless. Detect via MAPQ distribution.
3. Cross-sample contamination: detect via verifybamid2 or somalier (FREEMIX > 1% degrades somatic calling; > 5% breaks germline calling).
4. Wrong reference build: a BAM aligned to GRCh37 viewed against GRCh38 looks fine to flagstat but produces nonsense pileups. Compare @SQ M5: from BAM header with samtools dict ref.fa -- see alignment-validation.

Insert Size Caveats

samtools stats reports the IS section only for FR-oriented properly paired reads. So:

• Mate-pair libraries (RF orientation): IS section empty -- proper-pair flag not set for RF
• ATAC-seq: bimodal/multimodal expected (nucleosome ladder ~50/~180/~340 bp). Unimodal suggests poor transposition.
• RNA-seq: TLEN includes intron span -- mean meaningless
• Bisulfite (PBAT): orientation reversed; samtools may not flag proper pair

Related Skills

• sam-bam-basics - View alignment files; aligner-aware MAPQ semantics
• alignment-indexing - idxstats requires index; secondary+supp counted
• alignment-validation - Insert size by library, contamination, sample-swap detection
• duplicate-handling - Library-aware duplicate rate expectations
• alignment-filtering - Filter before stats
• sequence-io/sequence-statistics - FASTA/FASTQ statistics

Source & license

This open-source skill is cataloged on AgentStack and links to its original source — we do not rehost the code.

• Author: GPTomics
• Source: GPTomics/bioSkills
• License: MIT

Install and usage instructions live in the source repository linked above.