ExP Heatmap

An ordered-pair workflow for regional cross-population genomic visualization. ExP Heatmap takes a VCF through filter → prepare → compute → plot and renders ordered population-pair matrices as empirical rank-score heatmaps.

ExP Heatmap is designed for visualizing cross-population selection signals, differentiation, and other statistical summaries across many populations at once. It is most useful when a regional view (e.g. a gene locus or an ordered-pair matrix over a candidate window) would be hard to interpret as dozens of separate per-pair tracks. The package supports the canonical 1000 Genomes Project panel out of the box and now also handles arbitrary custom panels (e.g. the Gambian Genome Variation Project).

ExP heatmap of the human lactose (LCT) gene showing population differences between 26 populations from the 1000 Genomes Project, displaying empirical rank scores for cross-population extended haplotype homozygosity (XPEHH) selection test. Create your own LCT heatmap with the Quick Start Guide

Developed by the Laboratory of Genomics and Bioinformatics, Institute of Molecular Genetics of the Academy of Sciences of the Czech Republic

Features

End-to-end workflow: Built-in filter-vcf → prepare → compute → plot pipeline with no external preprocessing dependencies required
Multiple statistical tests: XP-EHH, XP-nSL, Delta Tajima's D, and Hudson's Fst
Flexible input: Work from raw VCF, pre-computed statistics, or ready-to-plot TSVs
Custom population panels: Automatic detection of non-1000G panels (e.g. GGVP) with fallback to the canonical 1000G layout when the full 26-population set is present
Correct pair-specific missingness: A NaN in one population pair no longer silently removes that locus from all other pairs
Wide-region static rendering: Explicit, documented column downsampling (max/mean/median) instead of implicit raster compression
Interactive HTML: Plotly-based zoom, pan, hover, region comparison, and population-focused views
Reproducibility scaffolding: Conda environment, Docker image, and CI-tested pytest suite
Benchmark harness: Scripts for local display, full-pipeline, and population-scaling benchmarks

Installation

Requirements

Python ≥ 3.10 (tested on 3.10, 3.11, 3.12)
vcftools (optional; only needed if you prefer external SNP-only preprocessing)

ExP Heatmap now includes a built-in filter-vcf command for preparing biallelic SNP-only VCFs locally, so vcftools is no longer required for the standard preprocessing path.

Python Dependencies

ExP Heatmap requires the following Python packages (automatically installed):

Package	Version	Purpose
scikit-allel	latest	Population genetics computations
zarr	< 3.0.0	Efficient array storage (v3+ not supported)
numpy	latest	Numerical operations
pandas	latest	Data manipulation
matplotlib	latest	Static visualizations
seaborn	latest	Heatmap rendering
click	latest	Command-line interface
plotly	latest	Interactive visualizations
tqdm	latest	Progress bars

Important: ExP Heatmap requires zarr < 3.0.0. The package will automatically install a compatible version, but if you encounter issues, ensure you have the correct version:
pip install 'zarr<3.0.0'

Install from PyPI

pip install exp_heatmap

Install from GitHub (latest version)

pip install git+https://github.com/bioinfocz/exp_heatmap.git

Quick Start

Get started with ExP Heatmap in three simple steps:

Step 1: Download the prepared results of the extended haplotype homozygosity (XPEHH) selection test for the part of human chromosome 2, 1000 Genomes Project data either directly via Zenodo or via command:

wget "https://zenodo.org/records/16364351/files/chr2_output.tar.gz"

Step 2: Decompress the downloaded folder in your working directory:

tar -xzf chr2_output.tar.gz

Step 3: Run the exp_heatmap plot command:

exp_heatmap plot chr2_output/ --start 136108646 --end 137108646 --title "LCT gene" --out LCT_xpehh

The exp_heatmap package will read the files from chr2_output/ folder and create the ExP heatmap and save it as LCT_xpehh.png file.

Usage

Pipeline Overview

ExP Heatmap follows a simple three-step workflow: prepare → compute → plot. If your input VCF still contains indels or multiallelic sites, use filter-vcf first. Each step can be used independently depending on your data format.

┌─────────────┐      ┌─────────────┐      ┌─────────────┐      ┌─────────────┐
│   VCF File  │ ──── │ filter-vcf  │ ──── │ SNP-only VCF│ ──── │   prepare   │
└─────────────┘      └─────────────┘      └─────────────┘      └──────┬──────┘
                                                                      │
                                                                      ▼
                                                               ┌─────────────┐
                                                               │   ZARR Dir  │
                                                               └──────┬──────┘
                                                                      │
                                                                      ▼
┌─────────────┐      ┌─────────────┐      ┌─────────────┐      ┌─────────────┐
│ Panel File  │ ──── │   compute   │ ──── │  TSV Files  │ ──── │   Heatmap   │
└─────────────┘      └─────────────┘      └──────┬──────┘      │ PNG / HTML  │
                                                 │             └─────────────┘
                                                 ▼
                                          ┌─────────────┐
                                          │    plot     │
                                          └─────────────┘

Starting Points:

From raw VCF: Use filter-vcf first if needed, then prepare → compute → plot
From VCF: Use all three steps (prepare → compute → plot)
From ZARR: Skip prepare, use compute → plot
From TSV results: Skip to plot directly

Command-Line Interface

1. VCF Filtering - `filter-vcf`

Keep only biallelic SNP records before prepare.

exp_heatmap filter-vcf [OPTIONS] <input_vcf>

Argument/Option	Type	Default	Description
`<input_vcf>`	PATH	required	Input VCF or VCF.GZ file
`-o, --out`	PATH	required	Output VCF path (`.vcf` or `.vcf.gz`)
`--no-log`	flag	-	Disable logging to file
`--verbose`	flag	-	Show detailed debug output in console

Example:

exp_heatmap filter-vcf chr21.vcf.gz -o chr21_snps.vcf

2. Full Pipeline - `full`

Run the complete pipeline (prepare → compute → plot) in a single command.

exp_heatmap full [OPTIONS] <vcf_file> <panel_file>

Argument/Option	Type	Default	Description
`<vcf_file>`	PATH	required	Recoded VCF file (SNPs only recommended)
`<panel_file>`	PATH	required	Population panel file
`-o, --out`	PATH	`exp_heatmap`	Prefix for all output files
`-s, --start`	INT	required	Start position for displayed region
`-e, --end`	INT	required	End position for displayed region
`-t, --test`	choice	`xpehh`	Statistical test: `xpehh`, `xpnsl`, `delta_tajima_d`, `hudson_fst`
`-c, --chunked`	flag	-	Use chunked array to avoid memory exhaustion
`--title`	STR	-	Title of the heatmap
`--cmap`	STR	`Blues`	Colormap for visualization
`--interactive`	flag	-	Generate interactive HTML visualization
`--rank-scores`	choice	`directional`	Rank-score mode: `directional`, `2-tailed` (legacy alias), `ascending`, `descending` (see Rank Score Options)
`--populations`	STR	inferred	Comma-separated population codes declaring the expected panel (see Declaring the population panel)
`--max-columns`	INT	auto/static, `30000` interactive	Explicit column budget for wide regions
`--column-aggregation`	choice	`max`	Reducer for static downsampling: `max`, `mean`, `median`
`--dpi`	INT	`400`	DPI for saved static figures
`--no-log`	flag	-	Disable logging to file
`--verbose`	flag	-	Show detailed debug output

Example:

exp_heatmap full chr15_snps.recode.vcf genotypes.panel -s 47924019 -e 48924019 -o slc24a5_analysis --title "SLC24A5"

This creates: slc24a5_analysis_zarr/, slc24a5_analysis_compute/, and slc24a5_analysis_plot.png

3. Data Preparation - `prepare`

Convert VCF files to efficient Zarr format for faster computation.

exp_heatmap prepare [OPTIONS] <vcf_file>

Argument/Option	Type	Default	Description
`<vcf_file>`	PATH	required	VCF file (biallelic SNP-only recommended)
`-o, --out`	PATH	`zarr_output`	Directory for ZARR output files
`--no-log`	flag	-	Disable logging to file
`--verbose`	flag	-	Show detailed debug output in console

Example:

exp_heatmap prepare chr15_snps.recode.vcf -o chr15.zarr

4. Statistical Analysis - `compute`

Calculate population genetic statistics across all genomic positions.

exp_heatmap compute [OPTIONS] <zarr_dir> <panel_file>

Argument/Option	Type	Default	Description
`<zarr_dir>`	PATH	required	Directory with ZARR files from `prepare` step
`<panel_file>`	PATH	required	Population panel file (see Input File Formats)
`-o, --out`	PATH	`output`	Directory for output TSV files
`-t, --test`	choice	`xpehh`	Statistical test: `xpehh`, `xpnsl`, `delta_tajima_d`, `hudson_fst`
`-c, --chunked`	flag	-	Use chunked array to avoid memory exhaustion
`--no-log`	flag	-	Disable logging to file
`--verbose`	flag	-	Show detailed debug output in console

Statistical Tests:

xpehh: Cross-population Extended Haplotype Homozygosity - detects recent positive selection
xpnsl: Cross-population Number of Segregating sites by Length - robust to variation in recombination rate
delta_tajima_d: Delta Tajima's D - measures difference in allele frequency spectrum
hudson_fst: Hudson's Fst - genetic differentiation between populations

Note: The -t flag has different meanings for different commands: for compute it specifies the statistical test (--test), while for plot it specifies the heatmap title (--title).

Example:

exp_heatmap compute chr15.zarr genotypes.panel -o chr15_results -t xpehh

Allele-frequency filtering

compute discards low-frequency variants before any statistic is calculated. This changes which loci reach the output, so the exact rule is documented here.


Rule	A variant is kept when its total alternate-allele frequency is strictly greater than 0.05
Applied	Once, across the whole callset, before population pairs are formed — not per population and not per pair
Frequency source	`INFO/AF` from the VCF when that field is present, summed across all ALT alleles. If the field is absent, frequencies are derived from the genotypes as 1 − (reference allele frequency), and a warning is logged
Configurable	No. The threshold is fixed in the current release

This is alternate-allele frequency, not minor-allele frequency. The two differ whenever the alternate allele is the major allele: a site with ALT AF = 0.98 has a minor-allele frequency of 0.02 and is retained, where a MAF filter at the same threshold would drop it. Filter upstream with vcftools --maf or bcftools view -q if you need minor-allele semantics.

The two frequency sources are not interchangeable. INFO/AF in a public release such as 1000 Genomes is the frequency across that release's full cohort, whereas the genotype-derived fallback uses only the samples present in your file. If you subset samples, the same 0.05 threshold selects a different set of variants depending on whether INFO/AF survived your preprocessing. When reproducibility matters, check which source was used.

compute reports what the filter did on every run, so the effect on your data is visible without extra flags:

Using precomputed alternate allele frequencies from variants/AF
Allele-frequency filter: kept <retained> of <total> variants (<removed> removed) at total alternate-allele frequency > 0.05

The first line names the frequency source and is replaced by a warning when INFO/AF is absent and genotypes are used instead.

Both lines are also written to the run log file (unless --no-log), so the threshold, the frequency source and the retained/removed counts form part of the provenance record for an analysis.

5. Visualization - `plot`

Generate heatmap visualizations from computed statistics.

exp_heatmap plot [OPTIONS] <input_dir>

Argument/Option	Type	Default	Description
`<input_dir>`	PATH	required	Directory with TSV files from `compute` step
`-s, --start`	INT	required	Start genomic position
`-e, --end`	INT	required	End genomic position
`-t, --title`	STR	-	Title of the heatmap
`-o, --out`	PATH	`ExP_heatmap`	Output filename (without extension)
`-c, --cmap`	STR	`Blues`	Colormap name (see Colormap Options)
`--interactive`	flag	-	Generate interactive HTML visualization
`--rank-scores`	choice	`directional`	Rank-score mode: `directional`, `2-tailed` (legacy alias), `ascending`, `descending` (see Rank Score Options)
`--populations`	STR	inferred	Comma-separated population codes declaring the expected panel (see Declaring the population panel)
`--max-columns`	INT	auto/static, `30000` interactive	Explicit column budget for wide regions
`--column-aggregation`	choice	`max`	Reducer for static downsampling: `max`, `mean`, `median`
`--dpi`	INT	`400`	DPI for saved static figures
`--no-log`	flag	-	Disable logging to file
`--verbose`	flag	-	Show detailed debug output in console

Example:

# Basic usage
exp_heatmap plot chr15_results/ --start 47924019 --end 48924019 --title "SLC24A5" --out slc24a5

# Interactive HTML output
exp_heatmap plot chr15_results/ --start 47924019 --end 48924019 --interactive --out slc24a5_interactive

# Static output with explicit wide-region downsampling
exp_heatmap plot chr15_results/ --start 47924019 --end 48924019 --max-columns 6000 --column-aggregation max --out slc24a5_binned

# Declare the expected panel so an incomplete compute directory fails loudly
exp_heatmap plot ggvp_results/ --start 15000000 --end 16000000 --populations "FULA,GWD,JOLA,MANDINKA,WOLOFF" --out ggvp_5pop

Declaring the population panel

The row order of an ExP heatmap is n(n−1) ordered pairs, blocked by first population and ordered by second population within each block, with self-pairs omitted. The population order used for those blocks is resolved in one of two ways.

Inferred (default). With no --populations, the population set is recovered from the compute-output filenames. If that set is exactly the canonical 26-population 1000 Genomes panel, the standard 1000G row order and super-population annotations are used. Otherwise the populations are ordered by first appearance in the filenames read lexicographically, which is alphabetical for any complete set of pairwise outputs.

Declared. Passing --populations "GWD,MSL,ESN" states the panel explicitly. That order is used directly, and the plotter verifies that every pairwise output for the declared panel is present. Here a compute run was interrupted before it reached the MANDINKA pairs:

$ exp_heatmap plot ggvp_results/ -s 15000000 -e 16000000 \
    --populations "FULA,GWD,JOLA,MANDINKA,WOLOFF"
Error: The declared population panel is incomplete in this compute directory:
4 of 10 population pairs are missing, with no output at all for MANDINKA.
Missing pairs: FULA_MANDINKA, GWD_MANDINKA, JOLA_MANDINKA, MANDINKA_WOLOFF.
Re-run compute for the full panel, or omit the declared panel to plot only the
populations that are present.

Why this matters. Outputs for a subset of populations form a self-consistent smaller panel, so with an inferred panel an incomplete compute directory renders without warning — and if it drops below the canonical 26, the row order, super-population annotations and color-scale bounds all switch to the custom-panel defaults. The same rank score then maps to a different color, making the figure incomparable to others from the same project. Declare the panel whenever figures will be compared across runs, or whenever a compute run may have been interrupted.

A declared panel matching the canonical 1000 Genomes order receives the full 1000G treatment, so --populations is safe to use on standard 1000 Genomes analyses. A declared panel listing the same populations in a different order is respected as given. Passing fewer populations than are present on disk is a supported way to plot a subset.

6. Superpopulation Summary - `summary`

Collapse population-pair rows to superpopulation pairs, write the collapsed TSV, and optionally create an interactive summary heatmap.

exp_heatmap summary [OPTIONS] <input_dir>

Example:

exp_heatmap summary chr15_results/ --start 47924019 --end 48924019 --agg-func mean --out slc24a5_summary

This creates:

slc24a5_summary.tsv
slc24a5_summary.html unless --no-plot is used

7. Population-Focused View - `focus`

Generate an interactive view containing only rows involving one chosen population.

exp_heatmap focus [OPTIONS] <input_dir>

Example:

exp_heatmap focus chr15_results/ --start 47924019 --end 48924019 --population CEU --out slc24a5_ceu_focus

8. Region Comparison - `compare`

Generate an interactive side-by-side comparison of two genomic regions from the same computed result set.

exp_heatmap compare [OPTIONS] <input_dir>

Example:

exp_heatmap compare chr15_results/ --start-1 47924019 --end-1 48924019 --start-2 55000000 --end-2 57000000 --out region_comparison

9. Top-Region Extraction - `regions`

Extract top-scoring windows into a TSV.

exp_heatmap regions [OPTIONS] <input_dir>

Example:

exp_heatmap regions chr15_results/ --start 47924019 --end 48924019 --n-top 25 --out slc24a5_top_regions.tsv

Python Package

The Python API offers more flexibility and customization options. Choose the appropriate scenario based on your data format:

Scenario A: Ready-to-Plot Data

Use when: You have pre-computed rank scores in a TSV file.

Data format: TSV file with columns: CHROM, POS, followed by pairwise columns for population comparisons.

from exp_heatmap.plot import plot_exp_heatmap
import pandas as pd

# Load your data
data = pd.read_csv("rank_scores.tsv", sep="\t")

# Create heatmap
plot_exp_heatmap(
    data,
    start=135287850,
    end=136287850,
    title="Population Differences in LCT Gene",
    cmap="Blues",
    output="lct_analysis",
    populations="1000Genomes"  # Predefined population set
)

Scenario B: Statistical Results to Rank Scores

Use when: You have computed statistical test results that need conversion to rank scores.

from exp_heatmap.plot import plot_exp_heatmap, create_plot_input

# Convert statistical results to empirical rank scores
data_to_plot = create_plot_input(
    "results_directory/",      # Directory with test results
    start=135287850, 
    end=136287850, 
    populations="1000Genomes",
    rank_scores="directional"  # Options: "directional" (legacy alias: "2-tailed"), "ascending", "descending"
)

# Create heatmap
plot_exp_heatmap(
    data_to_plot,
    start=135287850,
    end=136287850,
    title="XP-NSL Test Results",
    cmap="expheatmap",         # Custom ExP colormap
    output="xpnsl_results"
)

Scenario C: Complete VCF Workflow

Use when: Starting from raw VCF files. Combine CLI commands with Python plotting:

import subprocess
from exp_heatmap.plot import plot_exp_heatmap, create_plot_input

# 1. Prepare data (using CLI)
subprocess.run(["exp_heatmap", "prepare", "data_snps.recode.vcf", "-o", "data.zarr"])

# 2. Compute statistics (using CLI) 
subprocess.run(["exp_heatmap", "compute", "data.zarr", "populations.panel", "-o", "results/"])

# 3. Create custom plots (using Python)
data_to_plot = create_plot_input("results/", start=47000000, end=49000000)
plot_exp_heatmap(data_to_plot, start=47000000, end=49000000, 
                 title="Custom Analysis", output="custom_plot")

Additional Python API Functions

Declare the expected population panel:

plot(), plot_interactive() and create_plot_input() all accept a populations= argument, the API equivalent of --populations. Pass an iterable of codes to declare the panel; omit it to infer the panel from the compute-output filenames.

from exp_heatmap.plot import plot

# Declared: raises ValueError naming the missing populations if any
# pairwise output for this panel is absent from ggvp_results/
plot("ggvp_results/", start=15000000, end=16000000, title="GGVP",
     populations=("FULA", "GWD", "JOLA", "MANDINKA", "WOLOFF"))

See Declaring the population panel for what this guards against.

Summarize by Superpopulation:

from exp_heatmap.plot import create_plot_input, summarize_by_superpopulation

# Load data
data = create_plot_input("results/", start=47000000, end=49000000)

# Aggregate to superpopulation level (AFR, EUR, EAS, SAS, AMR)
superpop_data = summarize_by_superpopulation(data, agg_func='mean')
# Result has 20 rows (5×4 superpopulation pairs) instead of 650

Extract Top Regions:

from exp_heatmap.plot import create_plot_input, extract_top_regions

# Load data
data = create_plot_input("results/", start=47000000, end=49000000)

# Find genomic windows with highest selection signals
top_regions = extract_top_regions(data, n_top=50, window_size=10000)
print(top_regions[['center', 'mean_score', 'top_population_pair']])

Custom Colorbar Parameters:

from exp_heatmap.plot import prepare_cbar_params

# Calculate optimal colorbar settings based on data range
cmin, cmax, cbar_ticks = prepare_cbar_params(data_to_plot, n_cbar_ticks=6)

Reproducibility

Two local reproducibility entry points are included:

environment.yml: conda environment with Python 3.11, optional vcftools, and the package installed in editable mode with the dev extra
Dockerfile: minimal container image that installs the CLI and optional external preprocessing tools

Quick start:

conda env create -f environment.yml
conda activate exp-heatmap-dev

Or:

docker build -t exp-heatmap-local .
docker run --rm exp-heatmap-local --help

Input File Formats

VCF File

Standard VCF format. For best results:

Filter to biallelic SNPs only before prepare
Use the built-in CLI filter or your preferred external normalization tool

exp_heatmap filter-vcf input.vcf.gz -o input_snps.vcf

Panel File

Tab-separated file defining population membership for each sample. Required columns:

Column	Description
`sample`	Sample identifier (must match VCF sample names exactly)
`pop`	Population code (e.g., "CEU", "YRI")
`super_pop`	Superpopulation code (e.g., "EUR", "AFR")

Example panel file:

sample	pop	super_pop	gender
HG00096	GBR	EUR	male
HG00097	GBR	EUR	female
HG00099	GBR	EUR	female
NA18486	YRI	AFR	male
NA18487	YRI	AFR	female
NA18488	YRI	AFR	male

Important: The sample order in the panel file must match the sample order in the VCF/ZARR file exactly.

1000 Genomes Panel File:

Download the official panel file:

wget "ftp://ftp.1000genomes.ebi.ac.uk/vol1/ftp/release/20130502/integrated_call_samples_v3.20130502.ALL.panel" -O genotypes.panel

Output File Formats

Compute Output (TSV Files)

The compute step generates one TSV file per population pair, named POP1_POP2.tsv.

Rows remain position-aligned across all pair files. If one population pair yields no valid value at a locus, that locus is preserved as NaN for that pair rather than being globally removed from every output file.

Columns:

Column	Description
`name`	Dataset identifier (derived from input filename)
`variant_pos`	Genomic position of the variant
`{test}`	Raw test statistic value (e.g., `xpehh`, `xpnsl`)
`-log10_p_value_ascending`	Empirical rank score (ascending sort)
`-log10_p_value_descending`	Empirical rank score (descending sort)

Example output (CEU_YRI.tsv):

name	variant_pos	xpehh	-log10_p_value_ascending	-log10_p_value_descending
chr15_snps	48120990	0.523	1.234	2.456
chr15_snps	48120995	-0.891	2.891	1.123
chr15_snps	48121003	1.245	0.891	3.012

Note: Column names contain "p_value" for backward compatibility, but these are empirical rank scores, not classical p-values. See Statistical Methodology for details.

Advanced Features

Interactive Visualizations

Generate HTML-based interactive heatmaps with zoom, pan, and hover tooltips:

from exp_heatmap.interactive import plot_interactive

# Create interactive HTML visualization
plot_interactive(
    "results_directory/",
    start=135287850,
    end=136287850,
    title="Interactive LCT Analysis",
    output="lct_interactive"  # Saves as lct_interactive.html
)

Or via CLI:

exp_heatmap plot results/ --start 135287850 --end 136287850 --interactive --out lct_interactive

For static PNG output, ExP Heatmap now supports explicit wide-region downsampling. When the genomic window is wider than the effective export width, the static renderer can aggregate neighboring SNP columns using max, mean, or median before drawing the heatmap. The default static behavior uses an automatic column budget derived from figure width and DPI.

Additional Interactive Functions:

from exp_heatmap.interactive import create_comparison_view, create_population_focus_view
from exp_heatmap.plot import create_plot_input

data = create_plot_input("results/", start=40000000, end=60000000)

# Compare two genomic regions side-by-side
create_comparison_view(
    data,
    region1=(47000000, 49000000),
    region2=(55000000, 57000000),
    title="Region Comparison",
    output="comparison"
)

# Focus on comparisons involving a specific population
create_population_focus_view(
    data,
    focus_population="CEU",
    start=47000000,
    end=49000000,
    title="CEU Selection Signals",
    output="ceu_focus"
)

Advanced Customization

Fine-tune your visualizations with advanced options:

from exp_heatmap.plot import (
    create_plot_input,
    plot_exp_heatmap,
    prepare_cbar_params,
    superpopulations,
)

# Load the region, restricted to one superpopulation's panel
data_to_plot = create_plot_input(
    "results/",
    start=135000000,
    end=137000000,
    populations=superpopulations["AFR"],  # 7 populations -> 42 ordered rows
)

# Colorbar bounds and ticks fitted to this data
cmin, cmax, cbar_ticks = prepare_cbar_params(data_to_plot, n_cbar_ticks=6)

# Advanced plot with multiple customizations
plot_exp_heatmap(
    data_to_plot,
    # Pass the region as rendered rather than as requested: create_plot_input
    # snaps the window to the variant positions actually present in the data.
    start=int(data_to_plot.columns[0]),
    end=int(data_to_plot.columns[-1]),
    title="Selection Signals in African Populations",

    # Must match the panel used for create_plot_input above
    populations=superpopulations["AFR"],
    # Available: superpopulations["AFR" | "AMR" | "EAS" | "EUR" | "SAS"], or a custom list

    # Visual customizations
    cmap="expheatmap",                    # Custom ExP colormap
    display_limit=1.60,                   # Filter noise (values below limit set to 0)
    display_values="higher",              # Keep values above display_limit

    # Annotations
    vertical_line=[                       # Mark important SNPs
        (135851073, "rs41525747"),        # (position, label)
        (135851081, "rs41380347"),
    ],

    # Colorbar customization
    cbar_vmin=cmin,
    cbar_vmax=cmax,
    cbar_ticks=cbar_ticks,

    # Output
    output="african_populations_analysis",
    xlabel="Custom region description"
)

Color Scale Range

The colormap controls which colors are used; the color scale controls what range they are mapped over. ExP Heatmap picks that range in one of two ways.

Recognized 1000 Genomes panel — fixed range, 1.301 to 4.833. Applied so that heatmaps of the same panel stay directly comparable between runs: a given rank score maps to the same color in every figure.

Bound	Value	Meaning
Lower	`1.301`	−log₁₀(0.05), the conventional 5% tail cutoff. Cells below it are rendered as background, which is the noise suppression that makes the display readable.
Upper	`4.833`	A fixed ceiling. Cells above it are drawn at the maximum color.

Any other panel — fitted to the data. The range is taken from the minimum and maximum of the plotted values, because the displayed statistic and its range are not known in advance (Hudson's F_ST is bounded 0–1, XP-EHH is unbounded, Delta Tajima's D is signed).

Values above the upper bound are clipped

Empirical rank scores can reach log₁₀(n) for n ranked loci, and ranks are computed over the whole compute run, not the plotted window. A whole-chromosome run therefore has a ceiling well above 4.833 — for 588,621 loci it is 5.771 — and any cell above the fixed bound is drawn at the maximum color, indistinguishable from a cell exactly at it.

ExP Heatmap reports this when it happens:

322 of 1266200 rendered cells exceed the color-scale maximum (4.833); the data
reaches 5.771. Clipped cells all render at the top color and cannot be told
apart. Widen the scale to show them.

The affected fraction is small — about 10^−4.833 ≈ 1.5 × 10⁻⁵ of each pair's loci — but they are that pair's most extreme values. Whether that matters depends on your figure: it does not change where a signal appears, but it does flatten the top of a peak.

To keep the full range, override the bounds. Static output:

from exp_heatmap.plot import create_plot_input, plot_exp_heatmap, prepare_cbar_params

data = create_plot_input("results/", start=136108646, end=137108646)
cmin, cmax, cbar_ticks = prepare_cbar_params(data, n_cbar_ticks=5)   # fit to the data

# Pass the region as rendered, not as requested: create_plot_input snaps the
# window to the nearest variant positions present in the data.
plot_exp_heatmap(data,
                 start=int(data.columns[0]),
                 end=int(data.columns[-1]),
                 cbar_vmin=cmin, cbar_vmax=cmax, cbar_ticks=cbar_ticks,
                 output="lct_full_range")

Interactive output uses zmin/zmax on plot_interactive_heatmap() instead. Overrides are honored on both the 1000 Genomes and custom paths, and suppress the warning once the range covers the data.

Note on comparability. Fitting the scale to the data is the right choice for a single figure, but it means two figures no longer share a color mapping. Keep the fixed range when figures will be compared side by side, and widen it only when the top of the distribution is what you are reading.

Colormap Options

ExP Heatmap supports all matplotlib colormaps plus a custom colormap:

Colormap	Description	Best For
`Blues`	Sequential blue gradient (default)	General use
`expheatmap`	Custom colormap based on `gist_ncar_r` with white background	Highlighting strong signals
`gist_heat`	Heat-style gradient	High contrast visualization
`Reds`	Sequential red gradient	Alternative sequential
`viridis`	Perceptually uniform	Color-blind friendly
`RdBu`	Diverging red-blue	Bidirectional data

See the full list of matplotlib colormaps.

The colormap only chooses the colors. The range they are stretched over is set separately — see Color Scale Range.

Custom expheatmap colormap:

Based on gist_ncar_r
White background for low values (better noise filtering)
Dark blue for highest values
Optimized for selection signal visualization

Workflow Examples

Complete Analysis: SLC24A5 Gene

This example demonstrates a full workflow analyzing the SLC24A5 gene, known for its role in human skin pigmentation using 1000 Genomes Project data. SLC24A5 is also known to show strong selection signals, which makes it a suitable example.

#!/bin/bash

# Download 1000 Genomes data
wget "ftp.1000genomes.ebi.ac.uk/vol1/ftp/release/20130502/ALL.chr15.phase3_shapeit2_mvncall_integrated_v5b.20130502.genotypes.vcf.gz" -O chr15.vcf.gz
wget "ftp.1000genomes.ebi.ac.uk/vol1/ftp/release/20130502/integrated_call_samples_v3.20130502.ALL.panel" -O genotypes.panel

# Filter to biallelic SNPs only
exp_heatmap filter-vcf chr15.vcf.gz -o chr15_snps.vcf

# Prepare data
exp_heatmap prepare chr15_snps.vcf -o chr15_snps.zarr

# Compute statistics
exp_heatmap compute chr15_snps.zarr genotypes.panel -o chr15_snps_output

# Generate heatmap for SLC24A5 region
exp_heatmap plot chr15_snps_output \
    --start 47924019 \
    --end 48924019 \
    --title "SLC24A5" \
    --cmap gist_heat \
    --out SLC24A5_heatmap

Gallery

Different Rank Score Computations

The same XP-EHH test data for the ADM2 gene region, showing different rank score calculation methods:

Directional rank scores (reciprocal ordered-pair ranking; legacy alias: "2-tailed"):

Ascending rank scores (one-sided, highlights lowest test values):

Descending rank scores (one-sided, highlights highest test values):

Noise Filtering

Using display_limit and display_values parameters to filter noisy data and highlight the most extreme regions:

Same data as above, but with display_limit=1.60 to filter noise and highlight the strongest signals.

Statistical Methodology

Empirical Rank Scores

ExP Heatmap uses empirical rank scores to visualize selection signals. Despite column names containing "p_value" for backward compatibility, these values are not inferential p-values.

How rank scores are computed:

Rank all variants: For each population pair, sort all genome-wide test statistics.
Convert the rank to a fraction: rank_fraction = rank / total_variants
Transform to the display scale: empirical_rank_score = -log10(rank_fraction)

Interpretation:

Higher values indicate more extreme variants within that empirical ranking.
A value of 3.0 means the variant falls in the top 0.1% of the ranked distribution.
A value of 2.0 means the variant falls in the top 1% of the ranked distribution.

Rank Score Options

Which rank scores are displayed is selected by the --rank-scores option on the full, plot, summary, focus, compare and regions commands, and by the rank_scores parameter in create_plot_input(), plot() and plot_interactive(). Both accept the same four values:

Option	Description	Use Case
`directional`	For POP1_POP2: use descending; for POP2_POP1: use ascending	Default - Captures reciprocal directionality
`2-tailed`	Legacy alias of `directional`	Backward-compatible legacy name
`ascending`	Lowest test values ranked first	Detect negative selection signals
`descending`	Highest test values ranked first	Detect positive selection signals

Both rank-score columns are written by compute for every population pair, so switching mode only changes which column the plotting stage reads — no recomputation is needed.

Example:

# Direction-aware reciprocal ranking (recommended for most analyses)
exp_heatmap plot results/ --start 47000000 --end 49000000 --rank-scores directional --out lct_directional

# One-sided ranking for a specific hypothesis, from the same compute output
exp_heatmap plot results/ --start 47000000 --end 49000000 --rank-scores descending --out lct_descending

# Direction-aware reciprocal ranking (recommended for most analyses)
data = create_plot_input("results/", start=47000000, end=49000000, rank_scores="directional")

# One-sided ranking for a specific hypothesis
data = create_plot_input("results/", start=47000000, end=49000000, rank_scores="descending")

Statistical Tests Explained

Test	Detects	Interpretation
XP-EHH	Recent positive selection	Positive values: selection in pop1; Negative: selection in pop2
XP-nSL	Selection (robust to recombination rate variation)	Similar to XP-EHH but more robust
Delta Tajima's D	Difference in allele frequency spectrum	Positive: excess rare variants in pop1
Hudson's Fst	Population differentiation	Higher values: greater genetic distance

Benchmarks

All benchmarks were run locally on a macOS 14 / arm64 (Apple M4 Pro) workstation with 14 logical CPUs and 48 GB RAM. Python 3.12, zarr 2.x. Full benchmark scripts are provided under scripts/benchmarks/.

Full pipeline on GGVP chr21

Input: Gambian Genome Variation Project integrated chromosome 21 VCF (553,906 input records), filtered to 509,924 biallelic SNPs across 505 samples and 5 population labels (FULA, GWD, JOLA, MANDINKA, WOLOFF). Statistic: XP-EHH. Default settings otherwise.

Stage	Small case (1 Mb, 11,620 SNPs)	Full chromosome (509,924 SNPs)
`prepare`	0.99 s / 0.28 GB peak RSS	12.64 s / 0.38 GB peak RSS
`compute`	2.76 s / 0.21 GB peak RSS	74.73 s / 0.99 GB peak RSS
`plot` (static overview)	2.34 s / 0.26 GB peak RSS	2.27 s / 0.74 GB peak RSS

compute dominates wall-clock cost at full-chromosome scale, while static overview plotting remains cheap once per-pair TSVs exist. Peak RSS stayed below 1 GB for every stage in both cases.

Display scaling with population count

Synthetic scaling test that holds the genomic window fixed and varies the number of ordered population-pair rows by symlinking existing per-pair outputs. All runs use the GGVP chr21 compute output as the seed.

Populations	Ordered rows	Static plot	Interactive HTML	Peak RSS (interactive)
5	20	1.80 s, 0.10 MB PNG	1.83 s, 5.4 MB HTML	0.43 GB
10	90	3.25 s	3.57 s, 10.3 MB HTML	0.62 GB
20	380	9.61 s	9.24 s, 32.4 MB HTML	1.08 GB
30	870	16.87 s	16.34 s, 54.1 MB HTML	1.48 GB
50	2,450	38.36 s, 2.47 MB PNG	36.09 s, 99.5 MB HTML	2.17 GB

The ordered-pair display stays computationally tractable up to 50 populations, but visual density and HTML size rise quickly. At 50 populations we recommend using the focus, summary, or regions workflows in addition to the full matrix view.

Reproducing the benchmarks

The VCF, panel and compute directory are positional arguments.

pip install -e '.[benchmarks]'

# Full pipeline benchmark (prepare -> compute -> plot) on your own VCF.
# prepare is measured once; compute and plot are replicated.
python scripts/benchmarks/run_pipeline_benchmark.py \
    path/to/input.vcf.gz \
    path/to/panel.tsv \
    --out-dir local_data/benchmarks/pipeline \
    --repeats 5 --warmup 1

# Display scaling with population count, seeded from an existing compute output directory
python scripts/benchmarks/run_population_scaling.py \
    path/to/compute_output \
    --out-dir local_data/benchmarks/population_scaling \
    --repeats 5 --warmup 1 --dpi 600

# Display and reporting cost across window sizes and population counts
python scripts/benchmarks/run_local_benchmark.py \
    path/to/compute_output \
    --out-dir local_data/benchmarks/display \
    --repeats 5 --warmup 1 --dpi 600

--repeats reports mean, sample standard deviation, coefficient of variation and a 95% confidence interval (Student-t) over the measured runs; --warmup discards leading runs that pay filesystem-cache and import costs. Every run writes its own log, the individual per-run timings are kept in the seconds_runs column, and machine_specs.json records the CPU model, core counts, memory, platform and Python version alongside the results.

Note on --dpi: the scripts default to --dpi 200, which times a screen-resolution render. Pass the DPI you publish figures at if you want the reported display cost to match your figures.

Validation on Public Data

ExP Heatmap ships with two reproducible 1000 Genomes Project validation pipelines under scripts/validation/. Both scripts start from public Phase 3 release URLs and run end-to-end through filter-vcf → prepare → compute → plot.

chr15 / SLC24A5 (full-pipeline validation)

The run_1kg_chr15_slc24a5.sh script reproduces the pigmentation-locus showcase using raw public 1000 Genomes chromosome 15 data:

Stage	Wall-clock
`filter-vcf`	99.59 s (biallelic-SNP filter; retains 6,456,568 / 6,477,157 records)
`prepare`	252.98 s (VCF → Zarr)
`compute`	3034.42 s (XP-EHH, all 650 ordered population pairs from 26 populations)
`plot` (static)	28.43 s
`plot` (interactive)	14.04 s

Output artifacts include a 21.95 GB filtered VCF, 499.69 MB Zarr store, 2.94 GB of pairwise TSV results, and a 442 KB PNG / 21.54 MB HTML heatmap of the SLC24A5 window (chr15:47,924,019-48,924,019).

chr2 / LCT locus (region-scoped reconstruction)

The run_1kg_chr2_lct_reconstruction.sh script reconstructs the canonical lactase-persistence locus as a region-scoped public-data run.The script does not try to re-run whole-chromosome compute; instead it filters chromosome 2 to the plotted LCT window plus 1 Mb of flanking sequence on each side before preparing and computing.

Stage	Wall-clock
`filter-vcf` (chr2:135,108,646-138,108,646)	268.18 s (retains 81,595 / 7,081,600 records)
`prepare`	11.18 s
`compute`	117.41 s
`plot` (static)	4.61 s
`plot` (interactive)	2.11 s

This scales the whole-locus reconstruction to ~7 minutes of wall time end-to-end while keeping the plotted window and interpretation identical. See scripts/validation/summarize_validation_run.py to regenerate the JSON/Markdown summaries from log files.

1000 Genomes Population Reference

ExP Heatmap is optimized for the 1000 Genomes Project Phase 3 data (26 populations).

Population Codes

Code	Population	Superpopulation
ACB	African Caribbean in Barbados	AFR
ASW	African Ancestry in SW USA	AFR
ESN	Esan in Nigeria	AFR
GWD	Gambian in Western Division	AFR
LWK	Luhya in Webuye, Kenya	AFR
MSL	Mende in Sierra Leone	AFR
YRI	Yoruba in Ibadan, Nigeria	AFR
BEB	Bengali in Bangladesh	SAS
GIH	Gujarati Indians in Houston	SAS
ITU	Indian Telugu in the UK	SAS
PJL	Punjabi in Lahore, Pakistan	SAS
STU	Sri Lankan Tamil in the UK	SAS
CDX	Chinese Dai in Xishuangbanna	EAS
CHB	Han Chinese in Beijing	EAS
CHS	Han Chinese South	EAS
JPT	Japanese in Tokyo	EAS
KHV	Kinh in Ho Chi Minh City, Vietnam	EAS
CEU	Utah residents (CEPH) with European ancestry	EUR
FIN	Finnish in Finland	EUR
GBR	British in England and Scotland	EUR
IBS	Iberian populations in Spain	EUR
TSI	Toscani in Italy	EUR
CLM	Colombian in Medellín	AMR
MXL	Mexican Ancestry in Los Angeles	AMR
PEL	Peruvian in Lima	AMR
PUR	Puerto Rican in Puerto Rico	AMR

Superpopulations

Code	Name	Populations
AFR	African	ACB, ASW, ESN, GWD, LWK, MSL, YRI
SAS	South Asian	BEB, GIH, ITU, PJL, STU
EAS	East Asian	CDX, CHB, CHS, JPT, KHV
EUR	European	CEU, FIN, GBR, IBS, TSI
AMR	American (Admixed)	CLM, MXL, PEL, PUR

Reference: 1000 Genomes Project

Troubleshooting

Common Errors and Solutions

Zarr Version Error

Error:

Unsupported zarr version: 3.x.x
Please downgrade to zarr version < 3.0.0

Solution:

pip install 'zarr<3.0.0'

No Data Found in Genomic Region

Error:

ValueError: No data found in the requested genomic region (X - Y). 
The data contains positions from A to B.

Cause: The specified --start/--end coordinates don't overlap with the data.

Solution:

Check that your coordinates match the chromosome in your data
Use coordinates within the available range shown in the error message
Verify you're using the correct output directory

Strongest Signals All Render as One Flat Block

Warning:

322 of 1266200 rendered cells exceed the color-scale maximum (4.833); the data
reaches 5.771. Clipped cells all render at the top color and cannot be told apart.
Widen the scale to show them.

Cause: The recognized 1000 Genomes panel uses a fixed color scale of 1.301 to 4.833. Rank scores can reach log₁₀(n) for n ranked loci, so a whole-chromosome run exceeds the upper bound and the most extreme cells are all drawn at the maximum color. Nothing is wrong with the data — only the top of the color range is saturated.

Solution:

Leave it as is if you are reading where signals fall. The fixed range keeps the figure comparable with your other heatmaps, and the region still renders dark.
Widen the range if you are reading how strong the strongest signals are: pass cbar_vmin/cbar_vmax for static output, or zmin/zmax for interactive output. See Color Scale Range.
prepare_cbar_params(data) returns bounds and tick marks fitted to your data.

Note that the colorbar is labeled with the scale maximum, not the data maximum, so read the warning rather than the colorbar when checking whether a figure is saturated.

Heatmap Has Fewer Populations Than Expected

Symptoms: The heatmap renders successfully, but the y-axis lists fewer populations than the panel file contains. On 1000 Genomes data the super-population labels and boundary lines are also missing, and the color scale differs from earlier figures.

Cause: The compute output directory is missing every pairwise file for one or more populations — usually an interrupted compute run, or files moved or deleted afterwards. Those remaining outputs form a self-consistent smaller panel, so with an inferred panel the plotter has no way to know populations are absent and renders them out.

Solution:

Re-run plot with --populations listing the panel you expect. The plotter then reports exactly which populations and pairs are missing instead of rendering a reduced heatmap.
Confirm the directory holds n(n−1)/2 TSV files for n populations — 325 for the 26-population 1000 Genomes panel.
Re-run compute for the full panel if any are missing.

Related error:

The declared population panel is incomplete in this compute directory: ...

This is the same condition, caught because a panel was declared. The message names the missing populations and pairs.

Sample Order Mismatch

Error:

Sample order differs! Found X mismatches

Cause: The sample order in the panel file doesn't match the VCF/ZARR file.

Solution:

Ensure you're using the correct panel file for your VCF
Check that both files are from the same data release/phase
Verify sample IDs match exactly (case-sensitive)

Memory Exhaustion During Compute

Symptoms: Process killed, system becomes unresponsive, or "MemoryError"

Solution: Use the --chunked flag to process data in smaller chunks:

exp_heatmap compute data.zarr panel.tsv -o output/ --chunked

Empty or All-NaN Results

Symptoms: One pair-specific TSV file contains only NaN values in the statistic and rank-score columns.

Cause: The statistical test produced no valid results for that specific population pair, often due to:

Too few variants
Insufficient allele frequency variation
For delta_tajima_d: window size too large

Solution:

Check your data has sufficient variants
For delta_tajima_d, the default window size is 13 SNPs; ensure your data has enough variants
Inspect logs to identify which population pair failed; other pairs will still be preserved in the output directory

Getting Help

Check the GitHub Issues for known problems
Enable verbose logging with --verbose flag for detailed debug output
Log files are saved automatically (disable with --no-log)

Contributing

We welcome contributions! Feel free to contact us or submit issues or pull requests.

Development Setup

git clone https://github.com/bioinfocz/exp_heatmap.git
cd exp_heatmap
pip install -e ".[dev]"

Running the test suite

python -m pytest

The suite currently contains 28 tests covering rank-score generation with ties and missing data, pair-specific output preservation, AF fallback behavior, static downsampling, VCF filtering, custom-panel inference, GGVP metadata handling, and CLI wiring. GitHub Actions runs the suite on Python 3.10, 3.11, and 3.12 for every push and pull request.

Building documentation-facing assets

Reproducibility scripts, benchmark drivers, and validation pipelines live under scripts/.

License

This project is licensed under the MIT License.

Contributors

Edvard Ehler (@EdaEhler) - Lead Developer & Corresponding Author
Adam Nógell (@AdamNogell) - Maintainer
Jan Pačes (@hpaces) - Developer
Mariana Komárková (@satrovam) - Developer
Ondřej Moravčík (@ondra-m) - Original Developer

Acknowledgments

If you use ExP Heatmap in your research, please cite our paper [citation details will be added upon publication].

Name		Name	Last commit message	Last commit date
Latest commit History 247 Commits
.github/workflows		.github/workflows
assets		assets
scripts		scripts
src/exp_heatmap		src/exp_heatmap
tests		tests
.dockerignore		.dockerignore
.gitignore		.gitignore
Dockerfile		Dockerfile
LICENSE		LICENSE
README.md		README.md
environment.yml		environment.yml
pyproject.toml		pyproject.toml

Uh oh!

Folders and files

Latest commit

History

Repository files navigation

ExP Heatmap

Features

Table of Contents

Installation

Requirements

Python Dependencies

Install from PyPI

Install from GitHub (latest version)

Quick Start

Usage

Pipeline Overview

Command-Line Interface

1. VCF Filtering - filter-vcf

2. Full Pipeline - full

3. Data Preparation - prepare

4. Statistical Analysis - compute

Allele-frequency filtering

5. Visualization - plot

Declaring the population panel

6. Superpopulation Summary - summary

7. Population-Focused View - focus

8. Region Comparison - compare

9. Top-Region Extraction - regions

Python Package

Scenario A: Ready-to-Plot Data

Scenario B: Statistical Results to Rank Scores

Scenario C: Complete VCF Workflow

Additional Python API Functions

Reproducibility

Input File Formats

VCF File

Panel File

Output File Formats

Compute Output (TSV Files)

Advanced Features

Interactive Visualizations

Advanced Customization

Color Scale Range

Values above the upper bound are clipped

Colormap Options

Workflow Examples

Complete Analysis: SLC24A5 Gene

Gallery

Different Rank Score Computations

Noise Filtering

Statistical Methodology

Empirical Rank Scores

Rank Score Options

Statistical Tests Explained

Benchmarks

Full pipeline on GGVP chr21

Display scaling with population count

Reproducing the benchmarks

Validation on Public Data

chr15 / SLC24A5 (full-pipeline validation)

chr2 / LCT locus (region-scoped reconstruction)

1000 Genomes Population Reference

Population Codes

Superpopulations

Troubleshooting

Common Errors and Solutions

Zarr Version Error

No Data Found in Genomic Region

Strongest Signals All Render as One Flat Block

Heatmap Has Fewer Populations Than Expected

Sample Order Mismatch

Memory Exhaustion During Compute

Empty or All-NaN Results

Getting Help

Contributing

Development Setup

Running the test suite

Building documentation-facing assets

License

Contributors

1. VCF Filtering - `filter-vcf`

2. Full Pipeline - `full`

3. Data Preparation - `prepare`

4. Statistical Analysis - `compute`

5. Visualization - `plot`

6. Superpopulation Summary - `summary`

7. Population-Focused View - `focus`

8. Region Comparison - `compare`

9. Top-Region Extraction - `regions`