
Tooluniverse Population Genetics
- 99 installs
- 1.6k repo stars
- Updated August 4, 2026
- mims-harvard/tooluniverse
Helps with ai & agent building tasks.
About
tooluniverse-population-genetics is a Claude Code skill in the AI & Agent Building category.
- tooluniverse-population-genetics
- AI & Agent Building
- AI-coding skill
Tooluniverse Population Genetics by the numbers
- 99 all-time installs (skills.sh)
- +5 installs in the week ending Aug 4, 2026 (Skillselion tracking)
- Ranked #4,392 of 16,546 AI & Agent Building skills by installs in the Skillselion catalog
- Data as of Aug 5, 2026 (Skillselion catalog sync)
npx skills add https://github.com/mims-harvard/tooluniverse --skill tooluniverse-population-geneticsAdd your badge
Show developers this skill is listed on Skillselion. Paste this into your README.
| Installs | 99 |
|---|---|
| repo stars | ★ 1.6k |
| Last updated | August 4, 2026 |
| Repository | mims-harvard/tooluniverse ↗ |
What it does
Helps with ai & agent building tasks.
Files
Population Genetics Analysis
MC Strategy: Population genetics MC questions often test whether you know a specific theorem or result. COMPUTE the answer first (use popgen_calculator.py or write Python), then match to options. Don't try to reason about which option "sounds right."
Analyze population-level genetic variation, allele frequencies, GWAS associations, clinical significance, and evolutionary constraints using ToolUniverse tools.
When to Use
Activate this skill when the user asks about:
- Allele frequencies across populations (gnomAD, 1000 Genomes)
- GWAS associations for diseases/traits
- Clinical variant interpretation (ClinVar, VEP)
- Gene-level constraint metrics (pLI, LOEUF, o/e ratios)
- Selection, drift, linkage disequilibrium, or population structure
- Variant annotation and functional consequences
LOOK UP, DON'T GUESS
Query gnomAD/1000Genomes/GWAS Catalog FIRST for allele frequencies and associations. Preferred: use the PopGen_hwe_test, PopGen_fst, PopGen_inbreeding, and PopGen_haplotype_count tools for HWE, Fst, inbreeding, and haplotype calculations. Fallback: run popgen_calculator.py directly. For theoretical problems (delta-q, drift, LD decay), apply the formulas in the Theoretical Reasoning section below.
---
Tool Quick Reference
| Tool | Key Parameters | Notes |
|---|---|---|
gnomad_search_variants | query (REQUIRED) | Resolve rsID to variant_id format "CHR-POS-REF-ALT" |
gnomad_get_variant | variant_id (REQUIRED), dataset | Population frequencies. Default dataset: gnomad_r3; use gnomad_r4 for latest |
gnomad_get_gene_constraints | gene_symbol (REQUIRED) | pLI, o/e ratios. May timeout -- retry once |
MyVariant_query_variants | query (REQUIRED) | Aggregated: ClinVar + dbSNP + gnomAD + CADD. Uses hg19 coordinates |
EnsemblVEP_annotate_rsid | variant_id (REQUIRED) | Functional consequence, SIFT, PolyPhen. Param is "variant_id" NOT "rsid" |
EnsemblVEP_variant_recoder | variant_id (REQUIRED) | Convert between rsID/HGVS/VCF/SPDI |
gwas_get_snps_for_gene | gene_symbol (REQUIRED) | All GWAS SNPs for a gene |
gwas_search_associations | query (REQUIRED) | GWAS for a disease/trait (NOT gene name -- use gwas_get_snps_for_gene for genes) |
gwas_get_variants_for_trait | trait (REQUIRED) | Variants associated with a trait |
ClinVar_search_variants | gene, condition, significance | At least one filter required |
RegulomeDB_query_variant | rsid (REQUIRED) | Regulatory scoring (1a=strongest to 7=minimal) |
Critical Gotchas
1. gnomAD variant_id: Format is "CHR-POS-REF-ALT" (no "chr" prefix). Always resolve rsIDs via gnomad_search_variants first. 2. gwas_search_associations: Takes disease/trait names ONLY. Gene names will fail. Use gwas_get_snps_for_gene for gene-based lookups. 3. gwas_search_snps: BROKEN (HTTP 500). Use gwas_get_snps_for_gene instead. 4. VEP/ClinVar responses: Format is variable (list, {data, metadata}, or {error}). Handle all three.
---
Workflow Patterns
Variant frequency: gnomad_search_variants -> gnomad_get_variant(dataset="gnomad_r4") -> MyVariant_query_variants (1000G pop breakdowns) -> EnsemblVEP_annotate_rsid
GWAS for disease: gwas_search_associations -> gwas_get_variants_for_trait -> gnomad_get_variant for top hits -> EuropePMC_search_articles
Gene characterization: gnomad_get_gene_constraints -> gwas_get_snps_for_gene -> ClinVar_search_variants -> PubMed_search_articles
Pathogenicity assessment: EnsemblVEP_annotate_rsid -> MyVariant_query_variants (CADD, ClinVar) -> gnomad_get_variant (frequency) -> RegulomeDB_query_variant (if non-coding)
---
Theoretical Reasoning (CRITICAL for computation problems)
These formulas are needed for quantitative population genetics problems. Work through step by step, showing intermediate values.
Allele Frequency Change Under Selection (delta-q)
For a recessive deleterious allele (fitness: AA=1, Aa=1, aa=1-s):
delta_q = -s * q^2 * p / (1 - s * q^2)where p = freq(A), q = freq(a), s = selection coefficient.
For dominant deleterious (AA=1, Aa=1-s, aa=1-s):
delta_q = -s * q * p / (1 - s * q * (2 - q))For heterozygote advantage (AA=1-s1, Aa=1, aa=1-s2):
equilibrium: q_hat = s1 / (s1 + s2)Example: plug in s1 and s2 from the question; q_hat = s1/(s1+s2).
Selection against recessives is slow at low q because most a alleles hide in heterozygotes. Time to reduce q from q0 to qt: t ~ (1/qt - 1/q0) / s generations.
Genetic Drift in Small Populations
Variance in allele frequency per generation: Var(delta_p) = pq / (2Ne)
Probability of fixation of a new neutral mutation: 1/(2*Ne)
Time to fixation (given it fixes): ~4*Ne generations for neutral alleles
Heterozygosity decay: H_t = H_0 (1 - 1/(2Ne))^t
After t generations, fraction of heterozygosity lost ~ 1 - e^(-t/(2*Ne))
Effective population size (Ne) adjustments:
- Unequal sex ratio: Ne = 4NfNm / (Nf + Nm)
- Fluctuating size: Ne = harmonic mean of N across generations
- Bottleneck: dominated by the smallest generation size
Drift vs selection: Drift dominates when |s| < 1/(2*Ne). A variant with s=0.01 behaves neutrally in a population of Ne < 50.
Linkage Disequilibrium (LD) Decay
D = freq(AB) - freq(A)*freq(B), where A and B are alleles at two loci.
Decay with recombination: D_t = D_0 * (1 - r)^t, where r = recombination fraction, t = generations.
Half-life of LD: t_half = -ln(2) / ln(1-r) ~ 0.693/r generations (for small r).
r-squared (normalized LD): r^2 = D^2 / (pA pa pB * pb). Range 0-1.
Expected r^2 in finite population at equilibrium: E[r^2] = 1 / (1 + 4Ner) (for drift-recombination balance).
Practical implications:
- Tightly linked loci (r < 0.01): LD persists for hundreds of generations
- Loosely linked (r = 0.5, independent assortment): LD halves every generation
- GWAS tag SNPs work because LD extends over blocks; block size depends on Ne and recombination rate
- African populations have shorter LD blocks (larger historical Ne) -> need denser SNP arrays
Hardy-Weinberg Equilibrium
For alleles A (freq p) and a (freq q=1-p): expected genotypes AA=p^2, Aa=2pq, aa=q^2.
Chi-square test: df=1 (2 alleles). Preferred: use PopGen_hwe_test tool. Fallback: popgen_calculator.py --type hwe --AA N1 --Aa N2 --aa N3.
Causes of HWE departure: non-random mating, selection, migration, drift, genotyping error. Excess homozygotes -> inbreeding or population structure (Wahlund effect). Excess heterozygotes -> overdominant selection or negative assortative mating.
Heritability
- H^2 (broad-sense) = V_G / V_P; h^2 (narrow-sense) = V_A / V_P
- V_G includes ALL genetic variance: additive + dominance + epistasis. Trap: "broad-sense" is not just additive.
- Under HWE with two alleles (p, q): genotype frequencies are p^2, 2pq, q^2
- Phenotype frequency from genotype: sum(genotype_freq * penetrance) for each genotype class
- For quantitative traits: V_P = V_G + V_E (no covariance assumed)
- With dominance: assign genotypic values (e.g., AA=a, Aa=d, aa=-a), compute mean, then V_G from freq-weighted squared deviations
- PGS vs SNP-h² trap: PGS R² is NOT necessarily ≤ h²_SNP. With large GWAS, PGS can exceed SNP-h² by tagging rare causal variants through LD with common SNPs. The word "necessarily" makes this claim False. h²_SNP is estimated from common variants; PGS can capture additional variance.
Path Analysis (Causal Diagrams)
- Trace ALL paths from cause to effect through the diagram (direct + indirect)
- Each path's contribution = product of path coefficients along that path
- Total effect (correlation) = sum of contributions from all paths
- Indirect effects can mask (suppression) or inflate (confounding) the direct effect
- Unanalyzed correlations (double-headed arrows) count as valid path segments
- Never ignore indirect paths — the total is rarely just the direct arrow
Genetic Combinatorics (F2 crosses, haplotype counting)
For n SNPs between two inbred (homozygous) strains:
- F1 is heterozygous at all n loci
- F2 distinct haplotypes = 2^n (each SNP contributes parental A or B allele)
- F2 distinct diploid genotypes = 3^n (AA, AB, BB at each locus)
- F2 unique chromosomes (distinct haplotypes) = 2^n (e.g., 5 SNPs → 2^5 = 32; but subtract the 2 parental haplotypes if "novel" is asked → 30)
- ALWAYS write and run Python code (
python3 -c "...") for these counts. Never enumerate by hand. - For specimens/counting from field data: parse the data into a structure and compute programmatically.
Mutation-Selection Balance
Equilibrium frequency of a deleterious allele:
- Recessive lethal: q_hat = sqrt(mu/s) ~ sqrt(mu) when s=1
- Dominant lethal: q_hat = mu/s
- Example: mu=1e-5, s=1 (recessive lethal) -> q_hat = 0.003 (carrier freq ~ 0.006)
F-statistics and Population Structure
- Fis: Inbreeding within subpopulations (heterozygote deficit within demes)
- Fst: Differentiation between subpopulations. Fst = Var(p) / (p_bar * q_bar)
- Fit: Total inbreeding. (1-Fit) = (1-Fis)(1-Fst)
- Fst interpretation: <0.05 little, 0.05-0.15 moderate, 0.15-0.25 great, >0.25 very great differentiation
- Preferred: use
PopGen_fsttool. Fallback:popgen_calculator.py --type fst --p1 X --p2 Y --n1 N1 --n2 N2
---
Mendelian Genetics Reasoning Framework
For any genetics cross problem, follow these steps IN ORDER. Do not skip steps.
Step 1: Identify genes, locations, and allele relationships
- List every gene involved in the cross
- Determine chromosomal location: autosomal vs X-linked (X-linked genes show different inheritance in males vs females)
- Determine allele relationships: dominant/recessive, codominant, incomplete dominance
- Note any epistasis, suppressor, or modifier interactions between genes
Step 2: Write parental genotypes explicitly
- Use standard notation (e.g., Aa Bb for autosomal; X^w X^+ for X-linked)
- For X-linked genes, males are hemizygous (X^w Y), not homozygous
- If parental genotypes are not given, deduce them from phenotypes and pedigree context
Step 3: Draw Punnett square(s) for each gene
- For multi-gene crosses, handle each gene independently (if unlinked) then combine
- For linked genes, use recombination frequency to adjust gamete ratios
- For X-linked genes, remember that fathers pass X to all daughters and Y to all sons
Step 4: Calculate expected phenotypic ratios
- Multiply independent gene ratios (e.g., 3:1 x 3:1 = 9:3:3:1)
- For X-linked: calculate male and female ratios separately, then combine or report separately as required
Step 5: Verify ratios sum to 1.0
- Convert all ratios to fractions and confirm they sum to 1
- If they don't sum to 1, there is an error in the Punnett square or gamete calculation
Step 6: Apply phenotype modification rules AFTER computing genotypic ratios
- For epistasis: first compute the full genotypic ratios (e.g., 9:3:3:1), then collapse genotype classes that produce the same phenotype
- For suppressor genes: a suppressor homozygote (su/su) restores wild-type in an otherwise mutant background. Apply suppression AFTER determining which individuals carry the mutant allele
- Example: 9 A_B_ : 3 A_bb : 3 aaB_ : 1 aabb with recessive epistasis (aa masks B) becomes 9:3:4
---
E. coli Hfr Mapping Framework
For bacterial conjugation and Hfr mapping problems:
Core Principles
- In Hfr x F- crosses, the Hfr chromosome is transferred linearly starting from the origin of transfer (oriT)
- Gene transfer order = chromosomal order from the origin
- Early markers (entering first) are closest to the origin of transfer
- Late markers (entering last) are farthest from the origin
Interrupted Mating Experiments
- Genes that appear in recombinants at earlier time points are closer to oriT
- The time of entry gives the order and approximate distance between genes
- Recombinants require integration by homologous recombination (double crossover)
Recombination Frequency Between Markers
- KEY TRAP: Highest recombination frequency occurs between markers that are FARTHEST APART on the transferred segment
- This is because more time elapses between entry of distant markers, providing more opportunity for recombination events between them
- Conversely, markers that enter close together in time show LOW recombination between them
- Do NOT confuse "highest recombination frequency" with "first markers to enter" -- these are opposite concepts
Ordering Markers from Hfr Data
1. Use time-of-entry data to establish gene order relative to oriT 2. Use recombination frequency data between pairs of selected markers to confirm/refine order 3. Multiple Hfr strains with different origins can be used to build a circular map
---
MCQ Elimination Strategy for Genetics
General MCQ Protocol
1. ALWAYS evaluate ALL options before choosing an answer 2. Never select the first option that seems correct -- there may be a better or more precise answer 3. Read the question stem carefully for qualifiers: "MOST likely", "LEAST likely", "NOT true", "ALWAYS", "NEVER"
"Which is NOT true" Questions
- Evaluate EACH statement independently as True or False
- Mark each option with T or F before selecting
- The answer is the statement marked F
- Double-check: verify the "false" statement is genuinely false, not just misleadingly worded
"Which mechanism" Questions
- Test each proposed mechanism against ALL observations given in the question
- A correct mechanism must explain every observation, not just some
- Eliminate mechanisms that contradict even one observation
Specific Traps to Watch For
- Subfunctionalization vs neofunctionalization: Subfunctionalization = partitioning of EXISTING ancestral functions between duplicates (both copies needed to perform original function). Neofunctionalization = one copy acquires a genuinely NEW function not present in the ancestor
- Copy-neutral LOH: Caused by mitotic recombination (segmental, affects part of a chromosome), NOT uniparental disomy (UPD, which is whole-chromosome). The question may try to conflate these
- Penetrance vs expressivity: Penetrance = fraction of individuals with genotype who show ANY phenotype. Expressivity = degree/severity of phenotype among those who show it. These are distinct concepts
- Complementation vs recombination: Complementation = two mutations in DIFFERENT genes restore wild-type in trans. Recombination = exchange between two mutations in the SAME or different genes. Complementation is tested in F1 (heterozygote); recombination is tested in progeny
---
Common Genetics Reasoning Traps
These are specific patterns that have caused reasoning failures in hard genetics questions. Review before answering genetics MCQs.
Suppressor Genetics
- A suppressor mutation, when homozygous, restores wild-type phenotype in an otherwise mutant background
- In F2 crosses involving both the original mutation and an autosomal recessive suppressor:
- Treat as a dihybrid cross — the primary mutation and the suppressor segregate independently
- Only 1/4 of F2 are homozygous for the suppressor
- The suppressor only acts in individuals that are also homozygous for the primary mutation
- Use a Punnett square to enumerate all genotypic classes, then apply the suppression rule to determine phenotypes
Non-disjunction (Bridges' Experiments)
- Bridges used non-disjunction to prove the chromosome theory of inheritance
- X0 males arise from female meiosis non-disjunction events
- Meiosis I non-disjunction: both X chromosomes go to one pole -> XX egg + O egg (nullo-X)
- Meiosis II non-disjunction: sister chromatids fail to separate -> XX egg from one secondary oocyte
- The classic Bridges result: exceptional white-eyed females (X^w X^w) and red-eyed males (from nullo-X eggs + Y sperm = X0, but these are typically sterile)
- Key distinction: know which type of non-disjunction (MI vs MII) produces which specific gamete types
GWAS LD Blocks
- SNPs WITHIN the same LD block are correlated and can inflate false positive associations (one causal SNP drags along non-causal tag SNPs)
- SNPs ACROSS different LD blocks are largely independent and do NOT create misleading cross-locus associations
- LD block structure varies by population (shorter in African populations due to larger historical Ne)
- Fine-mapping within an LD block is needed to distinguish the causal variant from hitchhiking tag SNPs
Gene Retention After Whole-Genome Duplication
- Neofunctionalization: One copy acquires a NEW function -> most commonly cited reason for gene RETENTION after duplication (preserves both copies because each is now essential)
- Subfunctionalization: Ancestral functions are PARTITIONED between copies -> explains DIVERGENCE of duplicate copies, but both copies must be retained to maintain the full ancestral function
- Dosage balance: Some genes are retained in duplicate to maintain stoichiometric balance in protein complexes
- Trap: Questions may ask "what explains retention" vs "what explains divergence" -- these have different best answers
- For retention: neofunctionalization (new function makes both copies essential)
- For divergence of expression/function: subfunctionalization (partitioning of ancestral roles)
---
Advanced Genetics Traps v2
PGS vs Heritability: "Necessarily True" Logic
For "necessarily true" questions about PGS and heritability: a statement is necessarily true only if it holds when V_D=0 AND when V_D=V_G. Test the extremes.
Path Diagram Sign Assignment Protocol
Do NOT guess path signs from general knowledge. Signs may differ from well-known systems. Follow this protocol: 1. Establish reference direction: What varies? What is increasing? 2. For each path X→Y: Ask ONLY "when X increases, does Y increase (+) or decrease (-)?" 3. Use the question's experimental context (knockout/control comparisons, provided data) to determine signs — not intuition 4. Expect negative paths: Path diagrams test your ability to identify negative relationships. All-positive is almost always wrong. Direct residual paths (e) often have opposite sign from expectation.
Chi-Square: "Most Likely to Reject" Protocol
Compute chi-square from the expected ratio given in the question. Compare to chi-square-critical at df = (number of phenotype classes - 1). Pick the answer choice with the highest chi-square, but also check which pattern is biologically diagnostic of the alternative hypothesis.
LD and Misleading GWAS Associations
LD block boundaries at recombination hotspots are a source of GWAS false localization — strong signal in the block does not guarantee the causal variant is in the block.
Low-Frequency Allele Detection
Duplex sequencing (unique molecular identifiers + double-strand consensus) detects alleles at 0.01% frequency — far below standard NGS even at 80X depth. Simply increasing read depth does NOT help for ultra-rare variants because the Illumina error rate (~0.1%) masks variants rarer than ~1% regardless of depth. Error correction methods (UMIs, duplex consensus) are needed to distinguish true rare variants from sequencing errors.
---
Bundled Computation Script
Script: skills/tooluniverse-population-genetics/scripts/popgen_calculator.py
Preferred: Use ToolUniverse tools (via MCP/SDK) instead of the script when possible:
PopGen_hwe_testtool -- HWE chi-square test. Fallback:popgen_calculator.py --type hwePopGen_fsttool -- Weir-Cockerham Fst. Fallback:popgen_calculator.py --type fstPopGen_inbreedingtool -- Inbreeding coefficient from pedigree. Fallback:popgen_calculator.py --type inbreedingPopGen_haplotype_counttool -- Expected haplotype diversity. Fallback:popgen_calculator.py --type haplotypes
Fallback script modes (all require --type):
hwe:--AA N --Aa N --aa N-- chi-square HWE test with p-valuefst:--p1 F --p2 F --n1 N --n2 N-- Weir-Cockerham Fstinbreeding:--pedigree TYPE --generations G-- F from pedigree (self, full-sib, half-sib, first-cousin, etc.)haplotypes:--snps N --generations G --recomb_rate R-- expected haplotype diversity
---
Key Concepts
- MAF: Minor allele frequency. Common: >5%. Rare: <1%. Ultra-rare: <0.01%.
- pLI: P(LoF intolerant). >0.9 = haploinsufficient gene.
- LOEUF: LoF o/e upper fraction. <0.35 = highly constrained.
- CADD PHRED: >=10 top 10%, >=20 top 1%, >=30 top 0.1% most deleterious.
- Genome-wide significance: GWAS p < 5e-8 (Bonferroni for ~1M independent tests).
- Effect size: OR > 1 = risk allele, < 1 = protective. Beta > 0 = increases trait.
Evidence Grading
- T1: ClinVar pathogenic/likely pathogenic, FDA pharmacogenomics
- T2: gnomAD frequencies, GTEx eQTLs, GWAS genome-wide significant
- T3: CADD/SIFT/PolyPhen predictions, RegulomeDB, constraint metrics
- T4: VEP consequence terms, dbSNP annotations, literature mentions
#!/usr/bin/env python3
"""Population genetics calculator for HWE, Fst, inbreeding coefficients, and haplotypes.
Usage:
python popgen_calculator.py --type hwe --AA 100 --Aa 50 --aa 25
python popgen_calculator.py --type fst --p1 0.8 --p2 0.3 --n1 100 --n2 100
python popgen_calculator.py --type inbreeding --pedigree "half-sib" --generations 3
python popgen_calculator.py --type haplotypes --snps 5 --generations 3 --recomb_rate 0.5
"""
import argparse
import math
import sys
# ---------------------------------------------------------------------------
# HWE chi-square test
# ---------------------------------------------------------------------------
def hwe_test(obs_AA: int, obs_Aa: int, obs_aa: int) -> dict:
"""Test Hardy-Weinberg equilibrium for a biallelic locus.
Args:
obs_AA: Observed count of homozygous dominant genotypes.
obs_Aa: Observed count of heterozygous genotypes.
obs_aa: Observed count of homozygous recessive genotypes.
Returns:
dict with chi2, p_value, in_HWE flag, allele frequencies, expected counts.
"""
n = obs_AA + obs_Aa + obs_aa
if n == 0:
return {"error": "Total count is zero; no data to test."}
p = (2 * obs_AA + obs_Aa) / (2 * n)
q = 1.0 - p
exp_AA = n * p ** 2
exp_Aa = n * 2 * p * q
exp_aa = n * q ** 2
# Chi-square with 1 degree of freedom (HWE test)
chi2 = 0.0
for obs, exp in [(obs_AA, exp_AA), (obs_Aa, exp_Aa), (obs_aa, exp_aa)]:
if exp > 0:
chi2 += (obs - exp) ** 2 / exp
# Approximate p-value using chi-square CDF (1 df)
p_value = _chi2_pvalue(chi2, df=1)
return {
"obs_AA": obs_AA,
"obs_Aa": obs_Aa,
"obs_aa": obs_aa,
"exp_AA": round(exp_AA, 2),
"exp_Aa": round(exp_Aa, 2),
"exp_aa": round(exp_aa, 2),
"allele_freq_A": round(p, 4),
"allele_freq_a": round(q, 4),
"chi2": round(chi2, 4),
"p_value": round(p_value, 6),
"df": 1,
"in_HWE": p_value > 0.05,
"interpretation": _hwe_interpretation(chi2, p_value, obs_Aa, exp_Aa),
}
def _hwe_interpretation(chi2: float, p_value: float, obs_het: int, exp_het: float) -> str:
if p_value > 0.05:
return "Locus is in Hardy-Weinberg equilibrium (p > 0.05)."
direction = "excess" if obs_het > exp_het else "deficit"
return (
f"Significant deviation from HWE (chi2={chi2:.2f}, p={p_value:.4f}). "
f"Heterozygote {direction} observed. "
"Possible causes: selection, inbreeding, assortative mating, or genotyping error."
)
# ---------------------------------------------------------------------------
# Fst (Weir-Cockerham simplified)
# ---------------------------------------------------------------------------
def weir_cockerham_fst(p1: float, p2: float, n1: int, n2: int) -> dict:
"""Simplified Weir-Cockerham Fst between two populations.
Args:
p1: Allele frequency in population 1.
p2: Allele frequency in population 2.
n1: Sample size (individuals) for population 1.
n2: Sample size (individuals) for population 2.
Returns:
dict with Fst value and interpretation.
"""
for val, name in [(p1, "p1"), (p2, "p2")]:
if not 0.0 <= val <= 1.0:
return {"error": f"{name}={val} is out of range [0, 1]."}
if n1 <= 0 or n2 <= 0:
return {"error": "Sample sizes must be positive integers."}
p_bar = (n1 * p1 + n2 * p2) / (n1 + n2)
if p_bar == 0.0 or p_bar == 1.0:
return {
"p1": p1, "p2": p2, "n1": n1, "n2": n2,
"p_bar": round(p_bar, 4),
"Fst": 0.0,
"interpretation": "Allele is fixed or absent in both populations; Fst is undefined, set to 0.",
}
msp = (n1 * (p1 - p_bar) ** 2 + n2 * (p2 - p_bar) ** 2) / 1.0
msg = (n1 * p1 * (1 - p1) + n2 * p2 * (1 - p2)) / (n1 + n2 - 2)
denom = msp + msg
fst = max(0.0, (msp - msg) / denom) if denom > 0 else 0.0
return {
"p1": p1,
"p2": p2,
"n1": n1,
"n2": n2,
"p_bar": round(p_bar, 4),
"Fst": round(fst, 4),
"interpretation": _fst_interpretation(fst),
}
def _fst_interpretation(fst: float) -> str:
if fst < 0.05:
return f"Fst={fst:.4f}: Little genetic differentiation between populations."
if fst < 0.15:
return f"Fst={fst:.4f}: Moderate genetic differentiation."
if fst < 0.25:
return f"Fst={fst:.4f}: Great genetic differentiation."
return f"Fst={fst:.4f}: Very great genetic differentiation (unusual for human populations)."
# ---------------------------------------------------------------------------
# Inbreeding coefficient from pedigree type
# ---------------------------------------------------------------------------
# Base inbreeding coefficients per generation for common pedigree types.
# F for one mating of that type; raised to the power of generations if requested.
PEDIGREE_BASE_F = {
"self": 0.5,
"full-sib": 0.25,
"half-sib": 0.125,
"first-cousin": 0.0625,
"double-first-cousin": 0.125,
"uncle-niece": 0.125,
"aunt-nephew": 0.125,
"half-first-cousin": 0.03125,
"second-cousin": 0.015625,
}
def inbreeding_coefficient(pedigree: str, generations: int = 1) -> dict:
"""Compute inbreeding coefficient F for a known pedigree type.
For iterated mating (generations > 1) uses the recurrence:
F_t = base_F + (1 - base_F) * F_{t-1}
which converges toward 1 with repeated inbreeding.
Args:
pedigree: One of the pedigree type strings (e.g., "half-sib").
generations: Number of successive inbreeding generations.
Returns:
dict with F value per generation and cumulative F.
"""
pedigree_key = pedigree.lower().strip()
if pedigree_key not in PEDIGREE_BASE_F:
available = ", ".join(sorted(PEDIGREE_BASE_F))
return {"error": f"Unknown pedigree type '{pedigree}'. Available: {available}"}
if generations < 1:
return {"error": "generations must be >= 1."}
base_f = PEDIGREE_BASE_F[pedigree_key]
# Compute F for each generation using the recurrence
f_per_gen = []
f_prev = 0.0
for g in range(1, generations + 1):
# Each generation: offspring of two relatives with F = f_prev
f_g = base_f + (1 - base_f) * f_prev
f_per_gen.append(round(f_g, 6))
f_prev = f_g
cumulative_f = f_per_gen[-1]
return {
"pedigree": pedigree_key,
"base_F_per_mating": base_f,
"generations": generations,
"F_per_generation": f_per_gen,
"cumulative_F": round(cumulative_f, 6),
"heterozygosity_retained": round(1.0 - cumulative_f, 6),
"interpretation": _inbreeding_interpretation(pedigree_key, generations, cumulative_f),
}
def _inbreeding_interpretation(pedigree: str, generations: int, f: float) -> str:
if f < 0.0625:
severity = "low"
elif f < 0.25:
severity = "moderate"
else:
severity = "high"
return (
f"After {generations} generation(s) of {pedigree} mating, "
f"F = {f:.4f} ({severity} inbreeding). "
f"{(1 - f) * 100:.1f}% of heterozygosity retained."
)
# ---------------------------------------------------------------------------
# Haplotype counting after recombination
# ---------------------------------------------------------------------------
def haplotype_count(n_snps: int, generations: int, recomb_rate: float) -> dict:
"""Estimate distinct haplotypes after recombination.
Starting from 2^n_snps possible haplotypes, each generation recombination
reshuffles combinations. We track:
- Maximum theoretical haplotypes = 2^n_snps
- Expected recombinants per generation per chromosome pair = (n_snps - 1) * recomb_rate
- Cumulative probability of having recombined at least once across generations
Args:
n_snps: Number of biallelic SNP loci in the region.
generations: Number of generations of random mating.
recomb_rate: Per-adjacent-interval recombination rate (0-1) per generation.
Returns:
dict with theoretical max haplotypes, expected recombinants, and estimates.
"""
if n_snps < 1:
return {"error": "n_snps must be >= 1."}
if generations < 0:
return {"error": "generations must be >= 0."}
if not 0.0 <= recomb_rate <= 1.0:
return {"error": "recomb_rate must be between 0 and 1."}
max_haplotypes = 2 ** n_snps
n_intervals = max(n_snps - 1, 0)
# Expected number of crossovers per gamete per generation
expected_crossovers_per_gen = n_intervals * recomb_rate
# Probability that no recombination occurs anywhere across all intervals in one generation
p_no_recomb_one_gen = (1.0 - recomb_rate) ** n_intervals
# Over G generations: P(at least one recombination event somewhere)
p_at_least_one_recomb = 1.0 - p_no_recomb_one_gen ** generations
# Rough haplotype diversity estimate using decay of linkage disequilibrium.
# Each generation, LD between flanking markers decays as (1 - r)^G.
# Haplotypes diversify as: H(G) ~ H_max * (1 - (1-r)^G) + 2 (at least 2 ancestral)
if generations == 0:
estimated_haplotypes = 2 # ancestral
else:
ld_decay = (1.0 - recomb_rate) ** generations
diversity_fraction = 1.0 - ld_decay
estimated_haplotypes = max(2, int(round(2 + (max_haplotypes - 2) * diversity_fraction)))
estimated_haplotypes = min(estimated_haplotypes, max_haplotypes)
# Table of cumulative expected crossovers per generation
gen_table = []
cumulative_crossovers = 0.0
for g in range(1, generations + 1):
cumulative_crossovers += expected_crossovers_per_gen
p_no_recomb_cum = p_no_recomb_one_gen ** g
gen_table.append({
"generation": g,
"cumulative_expected_crossovers": round(cumulative_crossovers, 4),
"p_no_recomb_so_far": round(p_no_recomb_cum, 6),
})
return {
"n_snps": n_snps,
"n_intervals": n_intervals,
"recomb_rate_per_interval": recomb_rate,
"generations": generations,
"max_theoretical_haplotypes": max_haplotypes,
"expected_crossovers_per_gamete_per_gen": round(expected_crossovers_per_gen, 4),
"p_at_least_one_recombination": round(p_at_least_one_recomb, 6),
"estimated_distinct_haplotypes_after_G_gen": estimated_haplotypes,
"per_generation_table": gen_table,
"interpretation": _haplotype_interpretation(
n_snps, generations, recomb_rate, estimated_haplotypes, max_haplotypes
),
}
def _haplotype_interpretation(
n_snps: int, gens: int, r: float, estimated: int, max_h: int
) -> str:
pct = 100.0 * estimated / max_h if max_h > 0 else 0.0
if gens == 0:
return f"Before recombination: 2 ancestral haplotypes out of {max_h} theoretically possible."
return (
f"After {gens} generation(s) with r={r} per interval: "
f"~{estimated} of {max_h} possible haplotypes expected ({pct:.1f}% of maximum diversity). "
f"{'High LD retained.' if r < 0.1 else 'Substantial LD decay expected.'}"
)
# ---------------------------------------------------------------------------
# Chi-square p-value approximation (regularized incomplete gamma function)
# ---------------------------------------------------------------------------
def _chi2_pvalue(x: float, df: int = 1) -> float:
"""Approximate p-value for chi-square distribution using series expansion.
Uses the regularized upper incomplete gamma function Q(a, x) = 1 - P(a, x).
For df=1 this simplifies to: p = erfc(sqrt(x/2)).
"""
if x <= 0:
return 1.0
if df == 1:
return math.erfc(math.sqrt(x / 2.0))
# General case: use series / continued fraction via log-gamma
# P(df/2, x/2) approximated via regularized incomplete gamma
return _upper_incomplete_gamma_regularized(df / 2.0, x / 2.0)
def _upper_incomplete_gamma_regularized(a: float, x: float) -> float:
"""Q(a, x) = 1 - P(a, x) via series for small x, continued fraction for large x."""
if x < 0:
return 1.0
if x == 0:
return 1.0
if x < a + 1.0:
return 1.0 - _gamma_series(a, x)
return _gamma_continued_fraction(a, x)
def _gamma_series(a: float, x: float, max_iter: int = 200, eps: float = 1e-9) -> float:
"""Series representation of regularized lower incomplete gamma P(a,x)."""
log_gamma_a = math.lgamma(a)
ap = a
delta = term = 1.0 / a
total = term
for _ in range(max_iter):
ap += 1.0
term *= x / ap
total += term
if abs(term) < abs(total) * eps:
break
return total * math.exp(-x + a * math.log(x) - log_gamma_a)
def _gamma_continued_fraction(a: float, x: float, max_iter: int = 200, eps: float = 1e-9) -> float:
"""Continued fraction representation of Q(a, x)."""
log_gamma_a = math.lgamma(a)
fpmin = 1e-30
b = x + 1.0 - a
c = 1.0 / fpmin
d = 1.0 / b
h = d
for i in range(1, max_iter + 1):
an = -i * (i - a)
b += 2.0
d = an * d + b
if abs(d) < fpmin:
d = fpmin
c = b + an / c
if abs(c) < fpmin:
c = fpmin
d = 1.0 / d
delta = d * c
h *= delta
if abs(delta - 1.0) < eps:
break
return math.exp(-x + a * math.log(x) - log_gamma_a) * h
# ---------------------------------------------------------------------------
# CLI
# ---------------------------------------------------------------------------
def build_parser() -> argparse.ArgumentParser:
p = argparse.ArgumentParser(
description="Population genetics calculator: HWE, Fst, inbreeding, haplotypes.",
formatter_class=argparse.RawDescriptionHelpFormatter,
epilog=__doc__,
)
p.add_argument(
"--type",
required=True,
choices=["hwe", "fst", "inbreeding", "haplotypes"],
help="Calculation type.",
)
# HWE arguments
p.add_argument("--AA", type=int, default=0, help="[hwe] Observed AA genotype count.")
p.add_argument("--Aa", type=int, default=0, help="[hwe] Observed Aa genotype count.")
p.add_argument("--aa", type=int, default=0, help="[hwe] Observed aa genotype count.")
# Fst arguments
p.add_argument("--p1", type=float, default=None, help="[fst] Allele freq in population 1.")
p.add_argument("--p2", type=float, default=None, help="[fst] Allele freq in population 2.")
p.add_argument("--n1", type=int, default=None, help="[fst] Sample size for population 1.")
p.add_argument("--n2", type=int, default=None, help="[fst] Sample size for population 2.")
# Inbreeding arguments
p.add_argument(
"--pedigree",
type=str,
default=None,
help=(
"[inbreeding] Pedigree type: self, full-sib, half-sib, first-cousin, "
"double-first-cousin, uncle-niece, aunt-nephew, half-first-cousin, second-cousin."
),
)
p.add_argument("--generations", type=int, default=1, help="Number of generations (default 1).")
# Haplotype arguments
p.add_argument("--snps", type=int, default=None, help="[haplotypes] Number of SNP loci.")
p.add_argument(
"--recomb_rate",
type=float,
default=0.5,
help="[haplotypes] Recombination rate per adjacent interval per generation (default 0.5).",
)
return p
def _print_result(result: dict) -> None:
if "error" in result:
print(f"ERROR: {result['error']}", file=sys.stderr)
sys.exit(1)
for key, value in result.items():
if isinstance(value, list) and value and isinstance(value[0], dict):
print(f"{key}:")
for row in value:
row_str = " " + ", ".join(f"{k}={v}" for k, v in row.items())
print(row_str)
else:
print(f" {key}: {value}")
def main() -> None:
parser = build_parser()
args = parser.parse_args()
calc_type = args.type
if calc_type == "hwe":
result = hwe_test(args.AA, args.Aa, args.aa)
elif calc_type == "fst":
missing = [n for n, v in [("p1", args.p1), ("p2", args.p2), ("n1", args.n1), ("n2", args.n2)] if v is None]
if missing:
parser.error(f"--type fst requires: {', '.join('--' + m for m in missing)}")
result = weir_cockerham_fst(args.p1, args.p2, args.n1, args.n2)
elif calc_type == "inbreeding":
if args.pedigree is None:
parser.error("--type inbreeding requires --pedigree")
result = inbreeding_coefficient(args.pedigree, args.generations)
elif calc_type == "haplotypes":
if args.snps is None:
parser.error("--type haplotypes requires --snps")
result = haplotype_count(args.snps, args.generations, args.recomb_rate)
else:
parser.error(f"Unknown type: {calc_type}")
print(f"\n=== popgen_calculator: {calc_type.upper()} ===")
_print_result(result)
if __name__ == "__main__":
main()