Bio Population Genetics Scikit Allel Analysis

scikit-allel Analysis

Python library for population genetics analysis with efficient array data structures.

Installation

pip install scikit-allel
# Optional: zarr for chunked storage
pip install zarr

Reading VCF Files

Load VCF

import allel

callset = allel.read_vcf('data.vcf.gz')

print(callset.keys())
# dict_keys(['samples', 'calldata/GT', 'variants/CHROM', 'variants/POS', 'variants/REF', 'variants/ALT', ...])

samples = callset['samples']
genotypes = callset['calldata/GT']
positions = callset['variants/POS']
chroms = callset['variants/CHROM']

Specify Fields

callset = allel.read_vcf('data.vcf.gz',
    fields=['samples', 'calldata/GT', 'variants/POS', 'variants/CHROM', 'variants/QUAL'])

callset = allel.read_vcf('data.vcf.gz', fields='*')  # All fields

callset = allel.read_vcf('data.vcf.gz',
    region='chr1:1000000-2000000',
    samples=['sample1', 'sample2'])

Large Files (Chunked)

import zarr

allel.vcf_to_zarr('large.vcf.gz', 'data.zarr', fields='*', overwrite=True)
callset = zarr.open('data.zarr', mode='r')
gt = allel.GenotypeArray(callset['calldata/GT'])

Genotype Arrays

GenotypeArray

gt = allel.GenotypeArray(callset['calldata/GT'])

print(gt.shape)  # (n_variants, n_samples, ploidy)
print(gt.n_variants)
print(gt.n_samples)

print(gt[0])  # Genotypes at first variant
print(gt[:, 0])  # All variants for first sample

Basic Operations

ac = gt.count_alleles()
print(ac.shape)  # (n_variants, n_alleles)

af = ac.to_frequencies()
is_segregating = ac.is_segregating()
gt_filtered = gt.compress(is_segregating, axis=0)

Missing Data

is_called = gt.is_called()
is_missing = gt.is_missing()

miss_per_variant = (~is_called).sum(axis=1)
miss_per_sample = (~is_called).sum(axis=0)

call_rate_variant = is_called.mean(axis=1)
call_rate_sample = is_called.mean(axis=0)

Allele Counts and Frequencies

ac = gt.count_alleles()
ac_ref = ac[:, 0]
ac_alt = ac[:, 1]

af = ac.to_frequencies()
maf = af.min(axis=1)

n_singletons = (ac[:, 1] == 1).sum()
n_doubletons = (ac[:, 1] == 2).sum()

By Population

subpops = {
    'pop1': [0, 1, 2, 3, 4],
    'pop2': [5, 6, 7, 8, 9]
}

ac_subpops = gt.count_alleles_subpops(subpops)

ac_pop1 = ac_subpops['pop1']
ac_pop2 = ac_subpops['pop2']

Haplotype Arrays

h = gt.to_haplotypes()
print(h.shape)  # (n_variants, n_haplotypes)
print(h.n_haplotypes)

ac_hap = h.count_alleles()

PCA

import allel
import numpy as np

gn = gt.to_n_alt(fill=-1)
gn_filtered = gn[is_segregating]
gn_imputed = np.where(gn_filtered < 0, 0, gn_filtered)

coords, model = allel.pca(gn_imputed, n_components=10, scaler='patterson')
print(coords.shape)  # (n_samples, n_components)

Plot PCA

import matplotlib.pyplot as plt

plt.figure(figsize=(8, 6))
plt.scatter(coords[:, 0], coords[:, 1], c=population_labels)
plt.xlabel('PC1')
plt.ylabel('PC2')
plt.savefig('pca.png')

Diversity Statistics

Heterozygosity

ho = allel.heterozygosity_observed(gt)
he = allel.heterozygosity_expected(ac, ploidy=2)

mean_ho = np.mean(ho)
mean_he = np.mean(he)

Nucleotide Diversity (Pi)

pi = allel.sequence_diversity(positions, ac)
print(f'Pi = {pi:.6f}')

windows = allel.moving_statistic(positions, statistic=lambda x: allel.sequence_diversity(x, ac), size=10000, step=5000)

Watterson's Theta

theta_w = allel.watterson_theta(positions, ac)
print(f'Theta_W = {theta_w:.6f}')

Site Frequency Spectrum

sfs = allel.sfs(ac[:, 1])

plt.figure(figsize=(10, 5))
allel.plot_sfs(sfs)
plt.savefig('sfs.png')

Folded SFS

sfs_folded = allel.sfs_folded(ac)

plt.figure(figsize=(10, 5))
allel.plot_sfs_folded(sfs_folded)
plt.savefig('sfs_folded.png')

Windowed Statistics

pos = np.array(positions)
windows = np.arange(0, pos.max(), 100000)

pi_windowed, windows_used, n_bases, counts = allel.windowed_diversity(pos, ac, size=100000, step=50000)

plt.figure(figsize=(14, 4))
plt.plot(windows_used[:, 0], pi_windowed)
plt.xlabel('Position')
plt.ylabel('Pi')
plt.savefig('pi_windows.png')

Sample Subsetting

pop1_idx = np.array([0, 1, 2, 3, 4])
pop2_idx = np.array([5, 6, 7, 8, 9])

gt_pop1 = gt.take(pop1_idx, axis=1)
gt_pop2 = gt.take(pop2_idx, axis=1)

ac_pop1 = gt_pop1.count_alleles()
ac_pop2 = gt_pop2.count_alleles()

Filter Variants

is_snp = callset['variants/is_snp']
is_biallelic = ac.max_allele() == 1
is_segregating = ac.is_segregating()
qual = callset['variants/QUAL']
is_high_qual = qual > 30

flt = is_snp & is_biallelic & is_segregating & is_high_qual

gt_filtered = gt.compress(flt, axis=0)
pos_filtered = positions[flt]

Complete Workflow Example

import allel
import numpy as np

callset = allel.read_vcf('data.vcf.gz', fields=['samples', 'calldata/GT', 'variants/POS'])
gt = allel.GenotypeArray(callset['calldata/GT'])
pos = callset['variants/POS']
samples = callset['samples']

ac = gt.count_alleles()
flt = ac.is_segregating() & (ac.max_allele() == 1)
gt = gt.compress(flt, axis=0)
pos = pos[flt]
ac = gt.count_alleles()

print(f'Variants after filtering: {gt.n_variants}')
print(f'Samples: {gt.n_samples}')
print(f'Nucleotide diversity: {allel.sequence_diversity(pos, ac):.6f}')
print(f'Mean Het observed: {allel.heterozygosity_observed(gt).mean():.4f}')

gn = gt.to_n_alt(fill=-1)
gn = np.where(gn < 0, 0, gn)
coords, model = allel.pca(gn, n_components=10, scaler='patterson')

Related Skills

• selection-statistics - Fst, Tajima's D, iHS with scikit-allel
• linkage-disequilibrium - LD calculations in Python
• variant-calling/vcf-basics - VCF format and bcftools

scikit-allel Analysis - Usage Guide

Overview

scikit-allel provides Python data structures and algorithms for population genetics analysis. It's ideal for custom analyses, interactive exploration in Jupyter, and integration with other Python tools.

Prerequisites

pip install scikit-allel

# Optional for large files
pip install zarr

Note: scikit-allel is in maintenance mode. For new projects, consider sgkit for long-term support.

Quick Start

Tell your AI agent what you want to do:

• "Calculate nucleotide diversity from my VCF"
• "Compute allele frequencies per population"
• "Run PCA on my genetic data in Python"
• "Calculate Fst between populations"
• "Analyze haplotype structure"

Example Prompts

Loading and Basic Statistics

"Load my VCF into scikit-allel and calculate allele frequencies"

"Compute nucleotide diversity (pi) across the genome"

"Calculate per-site heterozygosity for each sample"

Population Comparisons

"Calculate pairwise Fst between my three populations"

"Run PCA and plot the first two components"

"Compute Watterson's theta in sliding windows"

Selection Analysis

"Calculate Tajima's D in 10kb windows"

"Compute iHS scores across chromosome 2"

"Find regions with unusual allele frequency differentiation"

Large Data Handling

"Convert my large VCF to Zarr format for efficient access"

"Calculate statistics on a Zarr-backed dataset"

"Process my VCF in chunks to avoid memory issues"

What the Agent Will Do

1. Load VCF data into appropriate array structures 2. Subset data by samples/populations if specified 3. Calculate requested statistics 4. Handle windowing or genome-wide aggregation 5. Generate visualizations if requested 6. Return results as DataFrames or arrays

Tips

• Use zarr backend for VCFs larger than available RAM
• GenotypeArray is for diploid data; HaplotypeArray for phased haplotypes
• Convert to AlleleCountsArray early for faster frequency calculations
• Filter missing data before calculating statistics
• Many functions accept pos arrays for windowed calculations

Data Structures

Quick Reference

import allel
import numpy as np

callset = allel.read_vcf('data.vcf.gz')
gt = allel.GenotypeArray(callset['calldata/GT'])
ac = gt.count_alleles()

pi = allel.sequence_diversity(callset['variants/POS'], ac)
print(f'Nucleotide diversity: {pi:.6f}')

Memory Management

For large VCFs, use Zarr:

allel.vcf_to_zarr('large.vcf.gz', 'data.zarr', fields='*')

import zarr
callset = zarr.open('data.zarr', mode='r')

Resources

• scikit-allel Documentation
• scikit-allel GitHub
• sgkit (successor)