LDetect-lite

A modern, fast re-implementation of LDetect, a method for calculating approximately independent linkage disequilibrium (LD) blocks in the human genome. The algorithm is described in Berisa & Pickrell, 2016.

Installation

LDetect-lite is available through PyPI:

pip install ldetect-lite

Or, with uv:

uv tool install ldetect-lite

This installs three equivalent CLI entry points — ldetect-lite, ldetect, and ldl — so pick whichever is most convenient; examples below use ldetect.

The main ldetect run pipeline reads the VCF/BCF reference panel via cyvcf2, a core dependency installed automatically — no separate tabix binary or htslib system package is required to run the pipeline. However, the VCF/BCFs must be indexed before running ldetect run — tabix -p vcf (for .vcf.gz) or bcftools index (for .bcf), from htslib/bcftools — since region-based partition reads require one.

Optional (--generate-heatmap): install matplotlib with pip install "ldetect-lite[heatmap]"/uv tool install ldetect-lite[heatmap]. Generating covariance heatmaps requires a matplotlib install.

Development

Install from source

git clone https://github.com/adamyhe/ldetect-lite.git
cd ldetect-lite
uv sync --group dev

From a development checkout, run CLI commands through uv run so they use the managed environment.

Usage

End-to-end pipeline

ldetect run \
  --genetic-map chr2.interpolated_genetic_map.gz \
  --reference-panel 1000G.chr2.vcf.gz \
  --individuals eurinds.txt \
  --chromosome chr2 \
  --output-dir results/chr2/

This writes results/chr2/chr2-ld-blocks.bed — a BED file of approximately independent LD blocks.

Global options (before the subcommand):

• -v / --verbosity {debug,info,warning,error} — logging verbosity (default: info; use warning to silence progress messages, debug for full detail)

Options:

• --ne FLOAT — effective population size Ne used by the Wen & Stephens shrinkage estimator (default: 11418.0, the CEU/HapMap II value; reproduction configs may override this for non-European populations)
• --cov-cutoff FLOAT — LD pairs with absolute shrinkage correlation below this threshold are not written to disk, reducing storage (default: 1e-7)
• --covariance-cache {compact,full} — partition cache schema for ldetect run (default: compact). Compact caches write only canonical position pairs, shrink_ld, diagonals, and lookup indexes, which is enough for restartable matrix-to-vector, metric, and local-search steps. Use full when debugging or when later running full-matrix/heatmap readers.
• --covariance-compression {lzf,zstd} — HDF5 compression codec for covariance partitions (default: zstd). zstd is smaller and faster to read/write than lzf at equal precision — see docs/optimizations.md.
• --ld-kernel {bitpacked,uint8} — compact covariance pair-count backend (default: bitpacked). bitpacked uses packed haplotypes and popcounts. uint8 keeps the older array-sum backend available for reference and diagnostics; use it when requesting --covariance-cache full.
• --n-snps-bw-bpoints N — target mean number of SNPs between consecutive breakpoints; controls block granularity (default: 10000, following Berisa & Pickrell 2016). The target breakpoint count is ceil(n_snps / N - 1). Mutually exclusive with --n-bpoints.
• --n-bpoints N — directly specify the number of breakpoints, bypassing the --n-snps-bw-bpoints formula; useful when replicating a published analysis with a known block count
• --subset {fourier,fourier_ls,uniform,uniform_ls} — which of the four breakpoint sets to write to the BED file (default: fourier_ls; see docs/pipeline-steps.md step 4)
• --all-breakpoint-subsets — compute all four breakpoint sets in the JSON output. By default, run computes only the requested --subset and its dependencies to avoid unused local-search work.
• --workers N — parallel workers for the pipeline (default: 1); set to the number of available cores to speed up covariance calculation (step 2) and, unless overridden below, matrix-to-vector, local search, and metric scoring as well
• --matrix-workers N — override parallel workers for matrix-to-vector partition processing (default: inherit --workers)
• --local-search-workers N — override parallel workers for local search (default: inherit --workers). Higher values can multiply RAM use because each worker loads its own covariance window.
• --metric-workers N — override parallel workers for streaming metric row passes during breakpoint scoring (default: inherit --workers)
• --high-precision — use 50-digit Decimal arithmetic for local search instead of the default float path (slower; mainly useful for exact reference comparisons)
• --filter-window {scipy-periodic,symmetric} — Hanning window mode for breakpoint filtering (default: scipy-periodic). scipy-periodic matches modern SciPy's periodic Hann window (scipy.signal.get_window(..., fftbins=True)) and is correct for new analyses. symmetric uses np.hanning — not a legacy fallback, but the setting needed to reproduce published reference blocks generated under scipy <1.1 (e.g. Berisa & Pickrell's original 2015 analysis), where a scipy defect made periodic odd-length windows bit-identical to symmetric ones until it was fixed in scipy 1.1.0 (2018). See notes/findings/ldetect-original-reproduction.md for the full root-cause writeup.
• --filter-workers N — minimum Numba threads for adaptive single-filter convolutions during the Step 4 filter-width search (default: 1)

If --workers is greater than 1 and none of OMP_NUM_THREADS/OPENBLAS_NUM_THREADS/MKL_NUM_THREADS/NUMEXPR_NUM_THREADS/NUMBA_NUM_THREADS are set, run prints a startup warning: numpy/BLAS/numba otherwise size their own internal thread pools to the whole machine's core count, not --workers, which oversubscribes real CPUs when multiple ldetect run processes share a node (e.g. several jobs on the same Slurm allocation). Export those five variables to match --workers to avoid this — see examples/ldetect_original/Snakefile's run_ldetect rule for a worked example.

Each of the five stages (partition, covariance, matrix-to-vector, find-minima, extract-bpoints) can also be run individually, along with a covariance-summary inspection utility — see docs/pipeline-steps.md.

The overview command writes pipeline-overview.svg; the per-step schematic command writes SVG/PDF figures under schematics/plots/.

Interpolate genetic maps

Convert a recombination rate map (e.g. the deCODE map or HapMap-interpolated 1000G maps) to per-SNP genetic positions required by steps 1 and 2:

ldetect interpolate-maps \
  --snp-file snps.bed.gz \
  --genetic-map recombination_map.gz \
  --output chr2.interpolated_genetic_map.gz

Arguments:

• --snp-file PATH — bgzipped BED file of SNP positions (columns: chrom start end rs_id); typically extracted from a filtered VCF with bcftools query -f '%CHROM\t%POS0\t%POS\t%ID\n'
• --genetic-map PATH — gzipped recombination map; interpolation is used to assign a cM value to each SNP position
• --output PATH — gzipped output map in the 3-column format expected by steps 1 and 2 (rs_id position cM)
• --mode {point,interval,hapmap,macdonald-decode,macdonald-pyrho} (default: point) — interpolation algorithm:
  • point — treats --genetic-map as discrete (position, cM) points and linearly interpolates between the two points bracketing each SNP. Correct for point-sampled maps (e.g. HapMap-interpolated 1000G maps).
  • interval — treats each map row as the start of a genomic interval with its own recombination rate (Begin, rate_cM_Mb, cumulative_cM_at_End). Correct for interval-rate maps such as deCODE; feeding those into point mode silently uses the next interval's rate for SNPs in the current interval.
  • hapmap — treats each map row's cM as the cumulative genetic position at that row's physical position, with the row's rate applying to the following interval. Correct for pyrho/HapMap-format maps.
  • macdonald-decode / macdonald-pyrho — compatibility modes for reproducing MacDonald et al.'s R interpolation scripts, including their dataframe/indexing conventions. Use these only for replication diagnostics; use interval/hapmap for corrected coordinates.

Algorithm

The pipeline detects LD block boundaries by finding local minima in a smoothed diagonal-sum signal derived from the shrinkage LD covariance matrix:

1. Partition — chromosome split into ~5000-SNP overlapping windows at low-recombination boundaries
2. Covariance — Wen & Stephens shrinkage estimator applied to phased haplotypes; shrinks sample correlations toward the expected LD decay to reduce finite-sample noise
3. Matrix → vector — each covariance matrix reduced to a [position, diagonal_sum] signal; troughs correspond to LD block boundaries
4. Find minima — binary search for optimal Hanning-window filter width; scipy.signal.argrelextrema finds local minima; local search refines each breakpoint using sum of squared inter-block correlations as the quality metric
5. Extract — chosen breakpoint set written as BED

The available breakpoint sets are fourier and uniform (raw minima from Fourier-filtered and uniformly-spaced candidates), plus fourier_ls and uniform_ls (after local search refinement). fourier_ls is the recommended output.

Known limitations

ldetect-lite reproduces the published Berisa & Pickrell (2016) 1000 Genomes LD blocks exactly for ASN (all 22 autosomes) and AFR (all chromosomes except chr22), and matches EUR block counts and coverage exactly but with shifted internal boundaries on chr8–chr12. These two residual divergences (EUR chr8-12, AFR chr22) likely stem from an unidentified upstream input/provenance difference from the original authors' pipeline, not a bug in this implementation — an extensive diagnostic effort ruled out VCF release-version provenance, SNP filtering, genetic map family, Ne assignment, duplicate/cross-partition handling, and reference-BED integrity as causes. Reproducing this result requires --filter-window symmetric (see above); the scipy-periodic default reflects the window shape the original ldetect code actually intended, and reproduces exactly what a from-scratch analysis on modern scipy would produce, but not the original 2015-era output, due to a since-fixed scipy defect. MacDonald et al. (2022) reproduction is exact for substantive deCODE/EUR boundaries and high-concordance but not exact for pyrho maps; the remaining pyrho gap is documented as local-search near-tie sensitivity plus map-desert effects, not a planned code change. See notes/findings/ldetect-original-reproduction.md and notes/findings/macdonald2022-reproduction.md for the full writeups.

Pre-computed LD blocks

Pre-computed BED files for 1000 Genomes reference populations are available from in hg19 coordinates from the original LDetect data repository and in hg38 coordinates from a more recent effort by MacDonald et al., 2022.

BED files produced by our work will be released once the code base leaves alpha/we finish major breaking updates.