Generate alignment statistics using samtools flagstat, stats, depth, coverage, and mosdepth. Use when assessing alignment quality, calculating coverage, or generating QC reports.
Installer avec Codex ou Claude Copiez ce prompt, collez-le dans Codex, Claude ou un autre assistant, puis laissez-le vérifier la page du skill et l'installer pour vous.
Une commande directe contourne le prompt de vérification. Examinez la source avant de l'exécuter.
Generate alignment statistics using samtools flagstat, stats, depth, coverage, and mosdepth. Use when assessing alignment quality, calculating coverage, or generating QC reports.
Before using code patterns, verify installed versions match. If versions differ:
Python: pip show <package> then help(module.function) to check signatures
CLI: <tool> --version then <tool> --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.
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 < 2 * read_length, R1 and R2 overlap.# Default samtools depth double-counts overlap; -s deducts:
samtools depth -s input.bam
Without -s, doubled support inflates somatic VAFs at sites covered by overlapping pairs (especially in fragmented samples: FFPE, cfDNA). mosdepth does not double-count overlap. samtools mpileup and bcftools mpileup both enable overlap detection by default; pass -x to disable (long form --disable-overlap-removal in samtools, --ignore-overlaps in bcftools).
mosdepth excludes unmapped, secondary, QC-fail, and duplicate reads by default (--flag 1796); supplementary reads are NOT excluded (use --flag 3844 to drop them too). Configurable via --flag. Memory ~ 4 bytes x longest chrom (1 GB for human chr1, 12+ GB for axolotl). Does not honor base quality; use samtools depth -q INT if needed.
import pysam
with pysam.AlignmentFile('input.bam', 'rb') as bam:
total = mapped = paired = proper = 0for read in bam:
total += 1ifnot read.is_unmapped:
mapped += 1if read.is_paired:
paired += 1if read.is_proper_pair:
proper += 1print(f'Total: {total}')
print(f'Mapped: {mapped} ({mapped/total*100:.1f}%)')
print(f'Properly paired: {proper} ({proper/paired*100:.1f}%)')
Per-Chromosome Counts
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')
Calculate Depth at Position
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}')
Mean Depth in Region
import pysam
defmean_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:
returnsum(depths) / len(depths)
return0
depth = mean_depth('input.bam', 'chr1', 1000000, 2000000)
print(f'Mean depth: {depth:.1f}x')
Coverage Statistics
Goal: Compute coverage breadth and depth for a genomic region from a BAM file.
Approach: Iterate pileup columns in the region, count covered positions and accumulate depth, then derive percentages and means.
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)}')
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
<5%
5-15%
20-50%
20-50%
50-90% pre-consensus
(skip)
(use UMI)
10-30%
5-30%
n/a
20-60%
Proper pair 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/1/3/255
30-55
30-55
30-50
20-40
Mt fraction
0.1-2%
0.1-2%
<1%
<0.1%
<0.1%
varies
varies
<10% (Omni-ATAC goal; original Buenrostro-2013 libraries were often majority-mito)
<2%
n/a
varies
Mean MAPQ is misleading; the distribution is bimodal (0 and aligner-max). The fraction at MAPQ >= 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:
Adapter readthrough: short fragments (insert < 2 * read_length) sequence into adapter; aligners soft-clip the adapter portion and flag the read as MAPPED. Detect:
Off-target enrichment (capture/WES): detect via picard CollectHsMetrics PCT_OFF_BAIT or PCT_SELECTED_BASES.
Low-complexity pile-up: telomere/centromere reads mass at MAPQ-0; counted as mapped but useless. Detect via MAPQ distribution.
Cross-sample contamination: detect via verifybamid2 or somalier (FREEMIX > 1% degrades somatic calling; > 5% breaks germline calling).
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
