Metadata-Version: 2.4
Name: ldetect-lite
Version: 0.2.2
Summary: Fast calculation of approximately independent LD blocks
Project-URL: Homepage, https://github.com/adamyhe/ldetect-lite
Project-URL: Repository, https://github.com/adamyhe/ldetect-lite
Project-URL: Issues, https://github.com/adamyhe/ldetect-lite/issues
Project-URL: Changelog, https://github.com/adamyhe/ldetect-lite/releases
Author: Adam He
License: MIT License
        
        Copyright (c) 2026 Adam Youlin He
        
        Permission is hereby granted, free of charge, to any person obtaining a copy
        of this software and associated documentation files (the "Software"), to deal
        in the Software without restriction, including without limitation the rights
        to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
        copies of the Software, and to permit persons to whom the Software is
        furnished to do so, subject to the following conditions:
        
        The above copyright notice and this permission notice shall be included in all
        copies or substantial portions of the Software.
        
        THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
        IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
        FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
        AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
        LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
        OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
        SOFTWARE.
License-File: LICENSE
Keywords: GWAS,LD,bioinformatics,linkage-disequilibrium
Classifier: Development Status :: 3 - Alpha
Classifier: Environment :: Console
Classifier: Intended Audience :: Science/Research
Classifier: License :: OSI Approved :: MIT License
Classifier: Operating System :: OS Independent
Classifier: Programming Language :: Python :: 3
Classifier: Programming Language :: Python :: 3 :: Only
Classifier: Programming Language :: Python :: 3.11
Classifier: Programming Language :: Python :: 3.12
Classifier: Topic :: Scientific/Engineering :: Bio-Informatics
Requires-Python: >=3.11
Requires-Dist: cyvcf2>=0.31
Requires-Dist: h5py>=3.10
Requires-Dist: hdf5plugin>=4.4
Requires-Dist: numba>=0.57
Requires-Dist: numpy>=1.24
Requires-Dist: scipy>=1.10
Provides-Extra: heatmap
Requires-Dist: matplotlib>=3.7; extra == 'heatmap'
Provides-Extra: profiling
Requires-Dist: psutil>=5.9; extra == 'profiling'
Description-Content-Type: text/markdown

# LDetect-lite

[![PyPI](https://img.shields.io/pypi/v/ldetect-lite)](https://pypi.org/project/ldetect-lite/)
[![Tests](https://github.com/adamyhe/ldetect-lite/actions/workflows/tests.yml/badge.svg)](https://github.com/adamyhe/ldetect-lite/actions/workflows/tests.yml)
[![PyPI Downloads](https://static.pepy.tech/personalized-badge/ldetect-lite?period=total&units=INTERNATIONAL_SYSTEM&left_color=BLACK&right_color=GREEN&left_text=downloads)](https://pepy.tech/projects/ldetect-lite)

A modern, fast re-implementation of [LDetect](https://bitbucket.org/nygcresearch/ldetect), a method for calculating approximately independent linkage disequilibrium (LD) blocks in the human genome. The algorithm is described in [Berisa & Pickrell, 2016](https://academic.oup.com/bioinformatics/article/32/2/283/1743626).

## Installation

LDetect-lite is available through PyPI:

```bash
pip install ldetect-lite
```

Or, with [uv](https://docs.astral.sh/uv/):

```bash
uv add 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](https://github.com/brentp/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](https://www.htslib.org/)/[bcftools](https://samtools.github.io/bcftools/) — since region-based partition reads require one.

**Optional** (`--generate-heatmap`): install matplotlib with `pip install "ldetect-lite[heatmap]"`, or use `uv sync --extra heatmap` from a source checkout. Generating covariance heatmaps requires a matplotlib install.

### Development

Install from source

```bash
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

```bash
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`.
- `--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)

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`.

### Interpolate genetic maps

Convert a recombination rate map (e.g. the [deCODE map](https://www.science.org/doi/10.1126/science.aau1043) or [HapMap-interpolated 1000G maps](https://github.com/joepickrell/1000-genomes-genetic-maps)) to per-SNP genetic positions required by steps 1 and 2:

```bash
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}` (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`), matching MacDonald et al.'s R interpolation scripts ([`interpolate.R`](https://github.com/jmacdon/LDblocks_GRCh38/blob/master/scripts/interpolate.R)/[`interpolate_pyhro.R`](https://github.com/jmacdon/LDblocks_GRCh38/blob/master/scripts/interpolate_pyhro.R)). Required for interval-rate maps such as the deCODE map — feeding those into `point` mode silently uses the *next* interval's rate for SNPs in the *current* interval, an off-by-one bug that produced a ~0.001–0.003 cM error per SNP in earlier testing (see `notes/findings/macdonald2022-reproduction.md`).

## 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) likelystem 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. See `notes/findings/ldetect-original-reproduction.md` for the full writeup, and `notes/findings/macdonald2022-reproduction.md` for the equivalent status reproducing MacDonald et al. (2022)'s GRCh38 blocks.

## 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](https://bitbucket.org/nygcresearch/ldetect-data) and in hg38 coordinates from a more recent effort by [MacDonald et al., 2022](https://www.biorxiv.org/content/10.1101/2022.03.04.483057v2).

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