R/variant_processing.R
filter_variants_ld.RdIdentifies pairs of variants with absolute correlation exceeding a threshold and removes one from each pair using a greedy algorithm. Optionally also removes linearly dependent columns, in column order.
filter_variants_ld(
G,
ld_threshold = 0.9,
remove_lindep = TRUE,
prefer_keep = "lower_maf"
)Numeric matrix (n x p) of genotype dosage.
Numeric correlation threshold. Pairs with |r| > ld_threshold are considered redundant. Default 0.9.
Logical. If TRUE (default), also remove columns
that are exact linear combinations of other columns after LD pruning.
The columns are visited in order and a column is dropped when its
residual after projection onto the kept columns before it has norm
below 1e-07 times its own norm, which is the rule of qr()
with its default tolerance. Prevents singularity in downstream Z-score
computation and correlation matrices. Requires a matrix without
missing values.
Character controlling which variant to keep when breaking ties in LD pruning. One of:
"lower_maf" (default): keep the rarer variant
"higher_maf": keep the more common variant
"first": keep the earlier column (positional order)
Integer vector of column indices to KEEP.
Cost: for \(N\) samples and \(m\) variable columns, the correlation is
one BLAS cross-product, \(O(N m^2)\), on the centered and scaled columns
(the two-pass column variance, equal to var(), decides the
zero-variance exclusion and scales the columns). The greedy prune thresholds
the correlation once and keeps each variant's count of above-threshold
partners up to date as variants are removed, \(O(m^2 + k m)\) for
\(k\) removals. Memory holds the \(m \times m\) correlation and its
logical threshold mask (about 130 MB and 65 MB at \(m = 4000\)). A matrix
with missing values takes the pairwise-complete cor() path instead,
which is far slower; extract_variant_set imputes before this
step, so the scan never does.
The linear-dependence step is a blocked in-order Gram-Schmidt applied twice
(.independent_columns()), \(O(N r m)\) in BLAS calls for \(r\)
kept columns among the \(m\) visited, and it stops once \(r\) reaches
the number of rows with a nonzero entry. It keeps the columns that
LINPACK qr() keeps, which shifted every remaining column once per
dropped column and so dominated genes with more variants than samples.
set.seed(42)
n <- 100
# Three variants: v1 and v2 highly correlated, v3 independent
v1 <- rbinom(n, 2, 0.1)
v2 <- v1 # perfect correlation
v2[1:3] <- 2 - v2[1:3] # introduce small differences
v3 <- rbinom(n, 2, 0.05)
G <- cbind(v1, v2, v3)
keep <- filter_variants_ld(G, ld_threshold = 0.9)
keep # one of v1/v2 removed, v3 kept
#> v1 v2 v3
#> 1 2 3
G_pruned <- G[, keep, drop = FALSE]