snputils

The toolkit for population genetics workflows in Python

pip install snputils

You can also install with uv, conda, or pixi. Available on Bioconda.

snputils is a Python library for reading, processing, analyzing, and visualizing genetic variation, ancestry, phenotype data, and relatedness.

Read

Read and write genetic variation, local and global ancestry, phenotype, and IBD data. Supported formats include VCF, BCF, BGEN, PLINK BED and PGEN, MSP, FLARE, LANC, and ADMIXTURE.

Analyze

Explore population structure with PCA, mdPCA, and maasMDS. Calculate allele frequencies and F-statistics, perform genotype QC and relatedness analysis, and run GWAS, admixture mapping, and more.

Visualize

Visualization options include embedding plots, local ancestry displays, chromosome paintings, admixture charts, Manhattan plots, and Q–Q plots, all created directly from snputils objects.

Analysis and visualization

Analyze population structure, visualize local ancestry, and run admixture mapping in Python.

Population structure from 1KGP

Example that uses 150 real samples from five 1000 Genomes populations on chromosome 22, then estimates pairwise Hudson FST. The chromosome-22 estimates are illustrative rather than genome-wide. The first run downloads and caches about 206 MB.

                        import snputils as su
from snputils.processing import embedding_dataframe_from_model
from snputils.stats import fst

populations = ["YRI", "CEU", "GIH", "CHB", "PEL"]
snp = su.load_dataset(
    "1kgp",
    resource="phase3",
    chromosomes=22,
    populations=populations,
    samples_per_population=30,
    max_variants=15_000,
    maf=0.05,
    require_complete=True,
    snv_only=True,
)

# Or load your own VCF, BCF, PGEN, BED, or BGEN with:
# snp = su.read_snp("cohort.vcf.gz")

pca = su.PCA(snp, n_components=2, fitting="exact")
coords = embedding_dataframe_from_model(pca)
coords["population"] = snp.sample_fid

su.viz.plot_embedding(coords, hue="population", show=True)
fst_results = fst(snp, sample_labels=snp.sample_fid, block_size=1_000)
fst_results.rename(columns={
    "pop1": "population 1", "pop2": "population 2",
    "est": "FST estimate", "se": "standard error",
}).round(3).head()
                    

PCA of 150 samples from five 1000 Genomes Project populations on chromosome 22

Pairwise Hudson FST on chromosome 22

Pairwise Hudson FST estimates for selected 1000 Genomes populations on chromosome 22
population 1	population 2	FST estimate	standard error
CEU	CHB	0.108	0.014
CEU	GIH	0.031	0.006
CEU	PEL	0.094	0.014
CEU	YRI	0.137	0.012
CHB	GIH	0.072	0.009

Local ancestry across the genome

Create cohort and individual chromosome views from local ancestry data. Example with synthetic data, no files are downloaded.

                        import snputils as su

demo = su.build_synthetic_chromosome_painting_dataset(
    n_samples=10,
    windows_per_chromosome=40,
    seed=42,
    ancestry_map={"0": "AFR", "1": "EUR", "2": "EAS"},
)
lai = demo["laiobj"]

# Or load your own MSP, ancestry VCF, or LANC file with:
# lai = su.read_lai("cohort.msp.tsv")

su.viz.plot_lai(
    lai,
    sort=False,
    figsize=(18, 5.2),
    legend=True,
    legend_kwargs={"bbox_to_anchor": (0.5, 1.01), "frameon": False},
    fontsize={"xticks": 7, "yticks": 8, "xlabel": 12, "ylabel": 12, "legend": 10},
)
su.viz.chromosome_painting(
    lai,
    "paintings",
    sample_id="sample0",
)
                    

Synthetic cohort-level local ancestry plot for ten samples across chromosomes 1 through 22 and X

Individual chromosome painting for sample0 showing two haplotypes across chromosomes 1 through 22 and X

Admixture mapping

Synthetic admixture mapping example with 2,000 samples and 22,000 ancestry windows, the plot shows ancestry-specific associations.

                        import snputils as su

demo = su.build_synthetic_admixture_dataset(
    n_samples=2_000, n_windows=22_000, n_covariates=3, seed=42
)
phenotype = su.PhenotypeObject(
    samples=demo["sample_ids"],
    values=demo["phenotype"],
    phenotype_name="trait",
    quantitative=True,
)
results = su.run_admixture_mapping(
    phenotype,
    demo["laiobj"],
    covar=demo["covariates"],
    covar_variance_standardize=True,
    ci=0.95,
    adjust=True,
    keep_hla=True,
)

afr_results = results.query("ANCESTRY == 'AFR'")
su.viz.manhattan_plot(afr_results, figsize=(10.5, 4.8))
                    

Manhattan plot of AFR-specific local ancestry associations from a synthetic admixture-mapping analysis

Reader benchmark

snputils is the fastest and most memory-efficient Python reader that suports BED, PGEN, VCF, BCF, and BGEN.
Benchmark on chromosome 22 of the 1000 Genomes Project.

Wall-clock read time and peak memory on Python 3.12 with eight AMD EPYC 9684X CPU cores and 200 GB RAM. Methodology and code.

Developed in collaboration between

and more collaborators worldwide.