coverage_plot — genome coverage / depth QC from BAMs

Compute genome-wide mean depth and breadth from BAM files with samtools coverage (length-weighted across contigs), then plot per-sample depth, per-sample breadth, and a per-chromosome depth heatmap.

See theory.md for depth vs breadth, why the genome-wide mean must be length-weighted, the quality filters, and how to read the numbers (including the low-depth → use angsd_het handoff).

Two scripts

script	role
coverage_qc.sh	BAM(s) → per-contig tables + aggregated coverage_summary.tsv
plot_coverage.py	summary (+ per-contig tables) → depth / breadth / per-chrom figures

Dependencies

samtools (≥1.10) for the compute step; Python with pandas, numpy, matplotlib for plotting. Covered by the repo environment.yml.

Usage

# compute — accepts explicit BAMs (-b, repeatable), a directory (-d), or a list (-l)
bash coverage_qc.sh -d /path/to/bams -o coverage_qc -q 20 -Q 20
bash coverage_qc.sh -l bam_list.txt -o coverage_qc -r autosomes.txt   # restrict summary

# plot
python plot_coverage.py --summary coverage_qc/coverage_summary.tsv --outdir figs \
    --coverage-dir coverage_qc/persample \
    --metadata samples.tsv --id-col sample --group-col population \
    --chrom-map chrom_map.tsv

-r autosomes.txt (one contig per line) restricts the genome-wide summary to autosomes; per-contig tables keep every contig so a Z/A or X/A depth ratio can still be computed. Run either script with -h / --help for all options.

Output

• coverage_qc/persample/<sample>.coverage.tsv — raw samtools coverage per contig.
• coverage_qc/coverage_summary.tsv — sample mean_depth breadth_pct covbases total_bp n_contigs.
• figs/coverage.mean_depth.{png,pdf}, .breadth.{png,pdf}, .per_chrom_depth.{png,pdf}.

Smoke test

bash example/make_example.sh          # builds a tiny 2-contig BAM (needs samtools)
bash coverage_qc.sh -b example/toy.bam -o example/out
python plot_coverage.py --summary example/out/coverage_summary.tsv \
    --outdir example/figs --coverage-dir example/out/persample