Generate QC statistics from alignment files including mapping rates, read counts, coverage depth, and per-chromosome distributions.
# samtools
conda install -c bioconda samtools
# pysam and plotting
pip install pysam matplotlib
# For plot generation
conda install -c bioconda gnuplot
Tell your AI agent what you want to do: - "Get quick statistics for my BAM file" - "Calculate coverage depth across my genome" - "Generate a comprehensive QC report with plots" - "Check per-chromosome read distribution"
"Run flagstat on my BAM file and interpret the results"
"How many reads are mapped in sample.bam?"
"What percentage of reads are properly paired?"
"Show read counts per chromosome using idxstats"
"Calculate the mitochondrial contamination percentage"
"Check X/Y ratio for sex determination"
"Calculate mean coverage depth for my WGS sample"
"What percentage of the genome is covered at 20x?"
"Generate depth statistics for target regions"
"Generate samtools stats and create QC plots"
"Create a summary table of statistics for all my samples"
| Need | Command | Speed |
|---|---|---|
| Quick read counts | flagstat | <1 sec |
| Per-chromosome counts | idxstats | <1 sec |
| Full QC report | stats | 1-10 min |
| Coverage summary | coverage | 1-5 min |
| Per-position depth | depth | 5-30 min |
samtools flagstat sample.bam
Output interpretation: - mapped: Percentage of reads that aligned - properly paired: Both mates aligned in expected orientation/distance - duplicates: PCR/optical duplicates (after markdup) - singletons: Paired reads where only one mate mapped
samtools idxstats sample.bam
# Total mapped reads
samtools idxstats sample.bam | awk '{sum += $3} END {print sum}'
# Mitochondrial fraction
samtools idxstats sample.bam | awk '/^chrM/ {mt=$3} {total+=$3} END {printf "MT: %.2f%%\n", mt/total*100}'
# Sex check (X/Y ratio)
samtools idxstats sample.bam | awk '/^chrX/ {x=$3} /^chrY/ {y=$3} END {printf "X:Y = %.2f\n", x/(y+1)}'
samtools stats sample.bam > sample_stats.txt
grep "^SN" sample_stats.txt | cut -f 2-
# Generate plots
plot-bamstats -p stats_plots/ sample_stats.txt
samtools coverage sample.bam
samtools coverage -r chr1:1000000-2000000 sample.bam
samtools coverage -m sample.bam # ASCII histogram
samtools depth sample.bam > depth.txt
samtools depth -a sample.bam > depth_with_zeros.txt # Include zero positions
samtools depth -r chr1:1000000-1001000 sample.bam
# Mean depth
samtools depth sample.bam | awk '{sum+=$3; n++} END {print "Mean:", sum/n}'
# Coverage at thresholds
samtools depth sample.bam | awk '{total++; if($3>=10) c10++; if($3>=20) c20++} END {printf ">=10x: %.1f%% >=20x: %.1f%%\n", c10/total*100, c20/total*100}'
import pysam
def flagstat(bam_path):
stats = {'total': 0, 'mapped': 0, 'paired': 0, 'proper': 0, 'duplicate': 0}
with pysam.AlignmentFile(bam_path, 'rb') as bam:
for read in bam:
stats['total'] += 1
if not read.is_unmapped:
stats['mapped'] += 1
if read.is_paired:
stats['paired'] += 1
if read.is_proper_pair:
stats['proper'] += 1
if read.is_duplicate:
stats['duplicate'] += 1
stats['map_rate'] = stats['mapped'] / stats['total'] * 100
return stats
import pysam
def region_coverage(bam_path, chrom, start, end):
with pysam.AlignmentFile(bam_path, 'rb') as bam:
depths = [0] * (end - start)
for pileup in bam.pileup(chrom, start, end, truncate=True):
if start <= pileup.pos < end:
depths[pileup.pos - start] = pileup.n
covered = sum(1 for d in depths if d > 0)
mean_depth = sum(depths) / len(depths)
return {'length': len(depths), 'covered': covered, 'pct_covered': covered / len(depths) * 100, 'mean_depth': mean_depth}
Per-assay thresholds (mapping rate, duplicate rate, proper-pair rate, MAPQ, mitochondrial fraction across WGS / WES / RNA-seq / scRNA / ATAC / ChIP / long-read / aDNA) are in SKILL.md "QC Thresholds Are Assay-Specific" -- avoid duplicating tables across SKILL.md and usage-guide.md.
echo -e "Sample\tTotal\tMapped\tPaired\tDuplicates" > summary.tsv
for bam in *.bam; do
sample=$(basename "$bam" .bam)
samtools flagstat "$bam" | awk -v s="$sample" '
/in total/ {total=$1}
/^[0-9]+ \+ [0-9]+ mapped \(/ {mapped=$1}
/properly paired/ {paired=$1}
/^[0-9]+ \+ [0-9]+ duplicates$/ {dup=$1}
END {print s"\t"total"\t"mapped"\t"paired"\t"dup}
' >> summary.tsv
# Anchored patterns target the totals; samtools 1.13+ also emits primary mapped / primary duplicates lines.
done
idxstats requires an index:
samtools index sample.bam
samtools idxstats sample.bam
Use region filtering or sampling:
samtools depth -r chr1:1-10000000 sample.bam # Single chromosome
samtools view -s 42.1 -b -o sub.bam sample.bam && samtools depth sub.bam # Reproducible 10% sample