Working with LD references
pecotmr authors
2026-10-22
Source:vignettes/ld-infrastructure.Rmd
ld-infrastructure.RmdOverview
Almost everything else in pecotmr needs LD. Fine-mapping from summary statistics needs it, RSS weight learning needs it, colocalization needs it, and the heritability estimators are built on nothing else. The pipelines handle it for you, so this vignette is about the layer underneath: the two kinds of LD reference pecotmr reads LD from, computing LD from genotypes, checking it is usable, and reducing a variant set that is too correlated to model.
Two kinds of LD reference
Every pipeline reads LD from one place: the ldReference
of the collection it is given. That reference is one of two things (or a
mixture of either, below).
- A genotype panel, read lazily from disk with
readGenotypes(). LD is computed on demand for whichever variants a step needs, so a panel far larger than any one region costs nothing until it is read. It is the same structure aQtlDatasetholds its genotypes in. - An
LdMatrix: a correlation matrix you built, held in memory. pecotmr does not read pre-computed LD from disk; read the matrix with whatever suits its format and wrap it.
A genotype panel
panel <- readGenotypes(
plink2Prefix = file.path(
system.file("extdata", "ld_reference", "chr22", package = "pecotmr"),
"protocol_example.LD.chr22"
)
)
panel## class: RangedSummarizedExperiment
## dim: 175 1000
## metadata(2): genotypeSource ldProjection
## assays(1): dosage
## rownames(175): chr22:10685239:C:T chr22:10761227:G:A ...
## chr22:18731027:T:C chr22:18850784:T:G
## rowData names(4): SNP A1 A2 AF
## colnames(1000): S1 S2 ... S999 S1000
## colData names(0):The bundled panel is a stochastic LD sketch: 1,000
pseudo-samples projected from 60 real people, which reproduces their LD
without carrying their genotypes. pecotmr recognises one by its
.afreq sidecar and records the real sample size and allele
frequencies on the panel when it is read, because the projected dosages
say nothing about either:
S4Vectors::metadata(panel)$ldProjection## $nSamples
## [1] 60
head(SummarizedExperiment::rowRanges(panel)$AF)## [1] 0.058333 0.150000 0.083333 0.058333 0.050000 0.050000A MAC cutoff or DENTIST’s truncation rank counts the 60 real samples, and ldsc’s unbiased r^2 and HDL use the effective size of the sketch’s LD, 1 / (1/60 + 1/1000) \approx 57. Neither is the 1,000 columns.
computeLd() correlates any subset of a panel’s
variants:
## [1] 175 175A panel can also be reached through an LD meta file: one row per
chromosome or block, giving its coordinates and the genotype files
behind it, with paths relative to the meta file.
readGenotypes(ldMeta = , region = ) opens the row covering
a region.
An LdMatrix
ld <- LdMatrix(R, nSamples = 57)
ld## LdMatrix: 175 variants, panel size 57
## no allele frequenciesThe matrix’s dimnames here are chrom:pos:A2:A1 ids, so
the variants are parsed from them; otherwise pass variants,
a GRanges or a data frame with chrom,
pos, A1 and A2 (and optionally
SNP and AF). Correlations must be signed by
the dosage of A1. nSamples is the number of
samples the LD was estimated from, and AF, the frequency of
A1, is what the frequency filters and cTWAS’s weight
scaling read.
LdMatrix() checks the matrix’s shape and labels, its
unit diagonal, that its entries lie in [-1,
1], and its symmetry. It does not check positive
semi-definiteness, which costs an eigendecomposition;
checkLd() below does.
A matrix has no genotypes, so the steps that need them work differently:
| Step | With an LdMatrix
|
|---|---|
| LD for any variant set | read from the matrix |
| MAF / MAC filters, cTWAS scaling | need AF (and nSamples for MAC) |
| DENTIST, ldsc’s unbiased r^2, HDL | need its nSamples
|
| Missingness filter | refused: there are no genotypes to be missing |
| RAiSS | inverts the matrix instead of the genotype SVD |
Mixtures of LD references
Summary statistics meta-analysed across ancestries are correlated
through the weighted average of the cohorts’ LD. LdMatrix()
builds that average from one LdMatrix per cohort, matching
the components by position and allele and keeping the variants they all
cover:
half <- seq_len(ncol(panel) / 2)
ldA <- LdMatrix(computeLd(panel[, half]), nSamples = 30)
ldB <- LdMatrix(computeLd(panel[, -half]), nSamples = 30)
LdMatrix(list(ldA, ldB), mixtureWeights = c(0.7, 0.3))## LdMatrix: 175 variants, panel size 52
## mixture of 2 panels, weights 0.7, 0.3
## no allele frequenciesThe weights are typically each cohort’s share of the GWAS sample size.
GenotypeMixture() builds the same average from genotype
panels, but lazily: the panels stay on disk and are read only for the
variants a step asks about, so each keeps answering for its own samples.
A frequency or missingness filter judges every panel on its own and
drops a variant that fails in any of them. The mixture is a
MultiAssayExperiment with one experiment per cohort, so the
panels are named after their cohorts, and experiments() and
sampleMap() show what it holds.
GenotypeMixture(
list(EUR = panel[, half], AFR = panel[, -half]),
mixtureWeights = c(EUR = 0.7, AFR = 0.3)
)## GenotypeMixture: 175 variants from 2 genotype panels
## panels EUR, AFR; weights 0.7, 0.3; panel sizes 500, 500Using a reference
Either kind goes in the same slot, and from there every step reads it the same way:
gwas <- GwasSumStats(..., ldReference = panel) # or ldReference = ldWhen two collections meet in one pipeline — a QTL and a GWAS side in
colocalization, say — their references must be the same: two panels must
share their samples and allele coding, two LdMatrix objects
must agree on the LD they share, and a panel never matches a matrix.
Computing LD from genotypes
computeLd() takes a dosage matrix and returns the
correlation matrix:
set.seed(1)
X <- matrix(rnorm(100 * 20), 100, 20, dimnames = list(NULL, paste0("v", 1:20)))
X[, 5] <- X[, 4] + rnorm(100, sd = 0.01) # a near-duplicate pair
R2 <- computeLd(X)
dim(R2)## [1] 20 20method selects the estimator — "sample"
(default), "population", or "gcta" — and
shrinkage applies shrinkage toward the identity.
backend can hand the work to SNPRelate or snpStats when
those are installed.
Checking that LD is usable
A correlation matrix estimated from finite samples is usually not
positive definite, and methods that invert it or take its
eigendecomposition will fail or return nonsense. checkLd()
reports the diagnosis. Run it on the real panel rather than the toy
matrix above, which is well conditioned by construction:
diagnosis <- checkLd(R)
diagnosis[c(
"isPd",
"isPsd",
"minEigenvalue",
"nNegative",
"conditionNumber"
)]## $isPd
## [1] FALSE
##
## $isPsd
## [1] TRUE
##
## $minEigenvalue
## [1] -1.043346e-15
##
## $nNegative
## [1] 0
##
## $conditionNumber
## [1] 1301696989That is what a real reference looks like: positive
semi-definite but not positive definite, with a smallest
eigenvalue at rounding scale (-7e-16) and a condition
number around 1e9. The negative value is floating-point
noise on a rank-deficient matrix, not structure — but it is exactly what
makes an unguarded sqrt() or solve() fail, and
why LDER clamps eigenvalues before forming its weights.
method = "shrink" or "eigenfix" returns a
repaired matrix instead of a report:
fixed <- checkLd(R, method = "eigenfix")
c(before = diagnosis$minEigenvalue, after = checkLd(fixed$R)$minEigenvalue)## before after
## -1.043346e-15 1.995988e-09The smallest eigenvalue is lifted off zero, though on a matrix this rank-deficient it stays tiny — repair makes the decomposition safe, it does not restore rank that the reference panel never had. Repairing also changes the data, so prefer a better-conditioned variant set where you can, which is what the next section is for.
Reducing a correlated variant set
Two variants in near-perfect LD carry the same information, and a design containing both is singular. Three tools address that, differing in what they optimise.
ldPruneByCorrelation() drops variants until no pair
exceeds a correlation threshold, keeping the choice arbitrary among
equivalents:
pruned <- ldPruneByCorrelation(X, corThres = 0.8)
ncol(pruned$X.new)## [1] 19
pruned$filter.id## [1] 1 2 3 4 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20It returns the reduced matrix and the retained column indices, so a parallel table of statistics can be subset the same way.
ldClumpByScore() keeps the best-scoring variant
in each LD clump rather than an arbitrary one — pass PIPs or −log10
p-values as score and it behaves like PLINK’s clumping:
ldClumpByScore(X, score = pips, chr = chrom, pos = position, r2 = 0.2)enforceDesignFullRank() is the targeted version: rather
than thresholding everything, it drops variants iteratively until the
design plus its covariates reaches full column rank, and stops
there.
Xdup <- cbind(X, dup = X[, 1]) # exactly collinear
dim(enforceDesignFullRank(Xdup, C = matrix(1, nrow(Xdup), 1)))## [1] 100 19The C argument is the covariate matrix the design will
sit alongside — rank has to be assessed on the combined matrix, since a
variant can be independent of the other variants yet collinear with a
covariate. strategy chooses what to sacrifice: the most
correlated variant, the lowest-variance one, or the one least correlated
with the response.
Aligning summary statistics to a panel
The counterpart to all of this is on the summary-statistics side:
variants must exist in the panel, and their alleles must be oriented the
same way. filterVariantsByLdReference() and
mergeVariantInfo() do that, and are covered in the QC vignette, where the motivation lives.
Next steps
- QC for summary statistics — matching statistics to a panel.
-
Heritability —
buildLdEigen()/buildLdScore()turn per-block LD into the statistics the h2 estimators consume. - Fine-mapping and TWAS weights — the RSS methods that consume LD.
Session info
## R version 4.5.3 (2026-03-11)
## Platform: x86_64-conda-linux-gnu
## Running under: Ubuntu 24.04.5 LTS
##
## Matrix products: default
## BLAS/LAPACK: /home/runner/work/pecotmr/pecotmr/.pixi/envs/default/lib/libopenblasp-r0.3.34.so; LAPACK version 3.12.0
##
## locale:
## [1] LC_CTYPE=C.UTF-8 LC_NUMERIC=C LC_TIME=C.UTF-8
## [4] LC_COLLATE=C.UTF-8 LC_MONETARY=C.UTF-8 LC_MESSAGES=C.UTF-8
## [7] LC_PAPER=C.UTF-8 LC_NAME=C LC_ADDRESS=C
## [10] LC_TELEPHONE=C LC_MEASUREMENT=C.UTF-8 LC_IDENTIFICATION=C
##
## time zone: Etc/UTC
## tzcode source: system (glibc)
##
## attached base packages:
## [1] stats graphics grDevices utils datasets methods base
##
## other attached packages:
## [1] pecotmr_0.9.3
##
## loaded via a namespace (and not attached):
## [1] tidyselect_1.2.1 pgenlibr_0.7.0
## [3] dplyr_1.2.1 farver_2.1.2
## [5] Biostrings_2.78.0 S7_0.2.2
## [7] bitops_1.1-0 gdsfmt_1.46.0
## [9] fastmap_1.2.0 reshape_0.8.10
## [11] mathjaxr_2.0-0 digest_0.6.39
## [13] lifecycle_1.0.5 magrittr_2.0.5
## [15] compiler_4.5.3 rlang_1.3.0
## [17] sass_0.4.10 tools_4.5.3
## [19] yaml_2.3.12 SNPRelate_1.44.0
## [21] knitr_1.52 S4Arrays_1.10.1
## [23] htmlwidgets_1.6.4 bit_4.6.0
## [25] DelayedArray_0.36.1 plyr_1.8.9
## [27] RColorBrewer_1.1-3 abind_1.4-8
## [29] BiocParallel_1.44.0 withr_3.0.3
## [31] purrr_1.2.2 numDeriv_2016.8-1.1
## [33] BiocGenerics_0.56.0 desc_1.4.3
## [35] grid_4.5.3 stats4_4.5.3
## [37] susieR_0.16.6 ggplot2_4.0.3
## [39] scales_1.4.0 MASS_7.3-66
## [41] MultiAssayExperiment_1.36.2 SummarizedExperiment_1.40.0
## [43] cli_3.6.6 rmarkdown_2.32
## [45] metafor_5.0-1 crayon_1.5.3
## [47] ragg_1.5.2 generics_0.1.4
## [49] otel_0.2.0 RcppParallel_5.1.11-2
## [51] tzdb_0.5.0 BiocBaseUtils_1.12.0
## [53] cachem_1.1.0 stringr_1.6.0
## [55] metadat_1.6-0 parallel_4.5.3
## [57] XVector_0.50.0 matrixStats_1.5.0
## [59] vctrs_0.7.3 Matrix_1.7-6
## [61] jsonlite_2.0.0 IRanges_2.44.0
## [63] hms_1.1.4 S4Vectors_0.48.1
## [65] bit64_4.8.6 mixsqp_0.3-54
## [67] irlba_2.4.1 archive_1.1.13
## [69] systemfonts_1.3.2 tidyr_1.3.2
## [71] jquerylib_0.1.4 glue_1.8.1
## [73] pkgdown_2.2.1 codetools_0.2-20
## [75] stringi_1.8.9 gtable_0.3.6
## [77] GenomicRanges_1.62.1 quadprog_1.5-8
## [79] tibble_3.3.1 pillar_1.11.1
## [81] htmltools_0.5.9 Seqinfo_1.0.0
## [83] R6_2.6.1 zigg_0.0.2
## [85] textshaping_1.0.5 vroom_1.7.1
## [87] evaluate_1.0.5 lattice_0.23-1
## [89] Biobase_2.70.0 readr_2.2.0
## [91] tictoc_1.2.1 backports_1.5.1
## [93] Rsamtools_2.26.0 RhpcBLASctl_0.23-42
## [95] Rfast_2.1.5.2 bslib_0.12.0
## [97] Rcpp_1.1.2 SparseArray_1.10.10
## [99] nlme_3.1-170 checkmate_2.3.3
## [101] xfun_0.61 fs_2.1.0
## [103] MatrixGenerics_1.22.0 pkgconfig_2.0.3