vignettes/split-intuition.Rmd
split-intuition.RmdLANTERN splits variant dosages by local ancestry using one of two methods:
| Mode | Input | Splitting rule |
|---|---|---|
| Unphased (dosage) | Genotype matrix + per-variant ancestry code | Probabilistic (p1 … pK proportions estimated per variant) |
| Phased (haplotype) | Per-haplotype GT + per-haplotype ancestry | Deterministic — no probability estimation |
This vignette builds intuition for both methods using small, hand-built matrices — no VCF or MSP files involved. Part A covers the two-population case; Part B generalizes to K ≥ 3 populations (the same functions, just with K-population inputs).
vignette("toy-vcf-msp-example").vignette("real-data-chr19").split_diploid()
split_diploid() takes two matrices of identical
dimensions:
gt_genotype — integer matrix, rows =
variants, cols = samples. Values: 0 (hom-ref), 1 (het), 2
(hom-alt).ancestry — integer matrix, same shape.
Ancestry codes: 1 = EUR/EUR, 2 =
AFR/EUR (mixed het), 3 = AFR/AFR.The ancestry code summarises both haplotypes at the locus (e.g. from an RFMix MSP file collapsed into a per-variant diploid code).
set.seed(42)
n_var <- 10 # variants
n_samp <- 6 # samples
# Simulate integer genotypes (0/1/2)
gt <- matrix(
sample(0:2, n_var * n_samp, replace = TRUE, prob = c(0.6, 0.3, 0.1)),
nrow = n_var, ncol = n_samp
)
rownames(gt) <- paste0("var", seq_len(n_var))
colnames(gt) <- paste0("S", seq_len(n_samp))
# Simulate ancestry codes (1=EUR, 2=MIX, 3=AFR)
anc <- matrix(
sample(1:3, n_var * n_samp, replace = TRUE, prob = c(0.4, 0.2, 0.4)),
nrow = n_var, ncol = n_samp
)
result <- split_diploid(gt, anc)
cat("African dosage (first 4 variants × 6 samples):\n")
#> African dosage (first 4 variants × 6 samples):
round(result$african[1:4, ], 2)
#> [,1] [,2] [,3] [,4] [,5] [,6]
#> [1,] 2 0 2 1.00 0 0
#> [2,] 1 0 0 0.00 0 0
#> [3,] 0 0 0 0.00 0 0
#> [4,] 1 0 2 0.67 0 1
cat("\nEuropean dosage (first 4 variants × 6 samples):\n")
#>
#> European dosage (first 4 variants × 6 samples):
round(result$european[1:4, ], 2)
#> [,1] [,2] [,3] [,4] [,5] [,6]
#> [1,] 0 0 0 0.00 0 0
#> [2,] 1 1 0 1.00 0 0
#> [3,] 0 2 2 0.00 0 0
#> [4,] 0 0 0 0.33 2 0For each variant the function computes ancestry-specific proportions p1 (African) and p2 (European) from the diploid counts across all samples:
p1 = (2·N_AA_hom + N_AA_het + N_MIX_hom) / denominator
p2 = (N_MIX_hom + 2·N_EUR_hom + N_EUR_het) / denominatorwhere denominator = total_alt_alleles − N_MIX_het
(mixed-het samples are assigned stochastically via p1/p2 rather than
excluded from the denominator).
Ambiguous-dominated edge-case: when all alt
alleles are carried by mixed-het individuals (the extreme case),
split_diploid() alone (no gla argument, its
default) falls back to p1 = p2 = 0.5 to avoid dividing by zero. But the
more general problem isn’t just this all-ambiguous extreme — it’s any
variant where mixed-het carriers dominate the observations,
even with one or two unambiguous carriers still present to compute a
(misleadingly confident) ratio from. Part C below explains why
that’s worth improving on, and how;
ancestry_split() applies that improvement (GLA shrinkage)
automatically.
# A variant where only mixed-het (code 2) samples carry the alt allele
gt_edge <- matrix(c(1L, 0L, 1L), nrow = 1) # three samples: het, ref, het
anc_edge <- matrix(c(2L, 1L, 2L), nrow = 1) # all alt carriers are MIX
r <- split_diploid(gt_edge, anc_edge)
cat("Singleton MIX case — p1 = p2 = 0.5\n")
#> Singleton MIX case — p1 = p2 = 0.5
cat("African: ", r$african, "\nEuropean:", r$european, "\n")
#> African: 0.5 0 0.5
#> European: 0.5 0 0.5For matrices with mismatched samples or variant sets,
ancestry_split_dosage() wraps split_diploid()
with automatic overlap handling and monomorphic filtering — see
?ancestry_split_dosage.
split_haplotype()
When RFMix (or equivalent) provides per-haplotype local ancestry calls, LANTERN can deterministically assign each haplotype’s allele to the correct ancestry background. Inputs are four matrices of the same shape (variants × samples):
| Argument | Description |
|---|---|
gt_hap0 |
integer 0/1 allele for haplotype 0 |
gt_hap1 |
integer 0/1 allele for haplotype 1 |
anc_hap0 |
ancestry code for haplotype 0 |
anc_hap1 |
ancestry code for haplotype 1 |
pop_codes |
named vector mapping c(AFR=0L, EUR=1L) (RFMix
default) |
# 3 variants × 4 samples
gt0 <- matrix(c(1L, 0L, 1L, # variant 1: hap0 alleles for 3 samples
0L, 1L, 0L, # variant 2
1L, 1L, 0L), # variant 3
nrow = 3, ncol = 3)
gt1 <- matrix(c(0L, 1L, 0L,
1L, 0L, 1L,
0L, 1L, 1L),
nrow = 3, ncol = 3)
# Haplotype ancestry: 0 = AFR, 1 = EUR (RFMix convention)
a0 <- matrix(c(0L, 0L, 1L,
1L, 0L, 0L,
0L, 1L, 0L),
nrow = 3, ncol = 3)
a1 <- matrix(c(1L, 1L, 0L,
0L, 1L, 1L,
1L, 0L, 1L),
nrow = 3, ncol = 3)
phased_result <- split_haplotype(gt0, gt1, a0, a1)
cat("African dosage (phased):\n"); print(phased_result$african)
#> African dosage (phased):
#> [,1] [,2] [,3]
#> [1,] 1 1 1
#> [2,] 0 1 1
#> [3,] 0 0 0
cat("European dosage (phased):\n"); print(phased_result$european)
#> European dosage (phased):
#> [,1] [,2] [,3]
#> [1,] 0 0 0
#> [2,] 1 0 1
#> [3,] 1 1 1The rule is simple: for each cell, hap0 contributes its allele to the
AFR count if anc_hap0 == 0, to EUR if
anc_hap0 == 1, and nothing otherwise. Same for hap1 — no
probability estimation, just bookkeeping.
split_diploid_multi() and
split_haplotype_multi() generalize Part A’s functions to K
≥ 2 populations. K = 2 is numerically identical to
split_diploid()/split_haplotype() — see
?split_diploid_multi for a worked equivalence check.
split_diploid_multi()
The diploid ancestry code for K populations encodes which combination of two parental backgrounds is present. You supply two objects:
pure_codes — named integer vector,
one entry per population:
c(AFR = 3L, EUR = 1L, NAT = 4L)mixed_codes — data frame with
columns code, pop1, pop2. One row
per mixed-ancestry diploid type:
data.frame(
code = c(2L, 5L, 6L),
pop1 = c("AFR", "AFR", "EUR"),
pop2 = c("EUR", "NAT", "NAT")
)For K populations there are K(K−1)/2 mixed-pair types:
| K | Pure types | Mixed pairs | Total codes |
|---|---|---|---|
| 2 | 2 | 1 | 3 |
| 3 | 3 | 3 | 6 |
| 4 | 4 | 6 | 10 |
| 5 | 5 | 10 | 15 |
set.seed(2025)
N_VAR <- 12 # variants
N_SAMP <- 30 # samples
# Ancestry code system
# Pure: AFR=3, EUR=1, NAT=4
# Mixed: AFR/EUR=2, AFR/NAT=5, EUR/NAT=6
pure_codes <- c(AFR = 3L, EUR = 1L, NAT = 4L)
mixed_codes <- data.frame(
code = c(2L, 5L, 6L ),
pop1 = c("AFR", "AFR", "EUR"),
pop2 = c("EUR", "NAT", "NAT")
)
# Simulate rare-variant genotype matrix (mostly 0, rare 1/2)
gt3 <- matrix(
sample(0:2, N_VAR * N_SAMP, replace = TRUE, prob = c(0.88, 0.10, 0.02)),
nrow = N_VAR, ncol = N_SAMP
)
rownames(gt3) <- paste0("chr19:", seq(1e6, by = 500L, length.out = N_VAR))
colnames(gt3) <- paste0("S", sprintf("%02d", seq_len(N_SAMP)))
# Simulate ancestry codes reflecting a tri-admixed cohort
# Roughly: 30% AFR, 20% EUR, 10% NAT, 20% AFR/EUR, 10% AFR/NAT, 10% EUR/NAT
all_codes <- c(pure_codes, mixed_codes$code)
code_probs <- c(0.30, 0.20, 0.10, 0.20, 0.10, 0.10)
anc3 <- matrix(
sample(all_codes, N_VAR * N_SAMP, replace = TRUE, prob = code_probs),
nrow = N_VAR, ncol = N_SAMP
)
out <- split_diploid_multi(gt3, anc3, pure_codes, mixed_codes)
cat("Output list names:", names(out), "\n")
#> Output list names: AFR EUR NAT
cat("Dimensions (each matrix):", nrow(out$AFR), "variants x",
ncol(out$AFR), "samples\n\n")
#> Dimensions (each matrix): 12 variants x 30 samples
cat("AFR dosage (first 4 variants, 6 samples):\n")
#> AFR dosage (first 4 variants, 6 samples):
round(out$AFR[1:4, 1:6], 3)
#> [,1] [,2] [,3] [,4] [,5] [,6]
#> [1,] 0 0 0.000 0 0 0
#> [2,] 0 0 0.667 0 0 0
#> [3,] 0 0 0.000 1 0 0
#> [4,] 0 0 0.000 0 0 0
cat("EUR dosage (first 4 variants, 6 samples):\n")
#> EUR dosage (first 4 variants, 6 samples):
round(out$EUR[1:4, 1:6], 3)
#> [,1] [,2] [,3] [,4] [,5] [,6]
#> [1,] 0 0 0.000 0 0 0
#> [2,] 0 0 0.333 0 0 0
#> [3,] 0 0 0.000 1 0 1
#> [4,] 0 0 0.000 1 0 0
# Conservation check: per-element sum across pools never exceeds raw gt
total <- out$AFR + out$EUR + out$NAT
excess <- max(total - gt3)
cat("\nMax excess of pooled dosage over raw gt:", excess,
"(must be <= 0)\n")
#>
#> Max excess of pooled dosage over raw gt: 0 (must be <= 0)How p1, p2, p3 are estimated. For each variant, the function estimates population allele proportions from unambiguous observations only (pure-ancestry samples of any genotype, plus mixed hom-alts which each contribute one alt allele to each parent pool). Ambiguous mixed-ancestry heterozygotes are excluded from the estimation and then split using the result.
A mixed-ancestry heterozygote in pair (i, j) then receives: of its dosage in pool i and in pool j.
# Trace p_k for a single variant by hand
var_idx <- 3 # pick variant 3
gt_v <- gt3[var_idx, ]
anc_v <- anc3[var_idx, ]
# Pure AFR (code 3)
N_AA2 <- sum(anc_v == 3 & gt_v == 2)
N_AA1 <- sum(anc_v == 3 & gt_v == 1)
# Pure EUR (code 1)
N_EE2 <- sum(anc_v == 1 & gt_v == 2)
N_EE1 <- sum(anc_v == 1 & gt_v == 1)
# Pure NAT (code 4)
N_NN2 <- sum(anc_v == 4 & gt_v == 2)
N_NN1 <- sum(anc_v == 4 & gt_v == 1)
# Mixed hom-alts
N_AE2 <- sum(anc_v == 2 & gt_v == 2)
N_AN2 <- sum(anc_v == 5 & gt_v == 2)
N_EN2 <- sum(anc_v == 6 & gt_v == 2)
D <- 2*N_AA2 + N_AA1 + 2*N_EE2 + N_EE1 + 2*N_NN2 + N_NN1 +
2*N_AE2 + 2*N_AN2 + 2*N_EN2
p1 <- if (D > 0) (2*N_AA2 + N_AA1 + N_AE2 + N_AN2) / D else 0.5
p2 <- if (D > 0) (2*N_EE2 + N_EE1 + N_AE2 + N_EN2) / D else 0.5
p3 <- if (D > 0) (2*N_NN2 + N_NN1 + N_AN2 + N_EN2) / D else 0.5
cat(sprintf("Variant %d D=%d p_AFR=%.3f p_EUR=%.3f p_NAT=%.3f\n",
var_idx, as.integer(D), p1, p2, p3))
#> Variant 3 D=5 p_AFR=0.600 p_EUR=0.400 p_NAT=0.000
cat("p1+p2+p3 =", p1 + p2 + p3, "\n")
#> p1+p2+p3 = 1Ambiguous-dominated edge case. When the denominator
D = 0 for a given variant (all alt alleles are carried exclusively by
mixed-ancestry heterozygotes, with no pure-ancestry or mixed hom-alt alt
alleles), split_diploid_multi() alone (no gla
argument, its default) falls back to a 50/50 split within each affected
pair. That’s the extreme case, but the same unreliability shows up
whenever mixed-het carriers dominate the count even with D > 0 — Part
C revisits both.
# One sample, one variant: AFR/EUR mixed het, no other data
gt_s <- matrix(1L, nrow = 1, ncol = 1)
anc_s <- matrix(2L, nrow = 1, ncol = 1) # code 2 = AFR/EUR
r_s <- split_diploid_multi(gt_s, anc_s, pure_codes, mixed_codes)
cat("Singleton AFR/EUR het -> AFR:", r_s$AFR[1,1],
" EUR:", r_s$EUR[1,1],
" NAT:", r_s$NAT[1,1], "\n")
#> Singleton AFR/EUR het -> AFR: 0.5 EUR: 0.5 NAT: 0split_haplotype_multi()
When per-haplotype ancestry calls are available (e.g. from RFMix),
the split is deterministic: each haplotype’s allele
goes directly to its population pool. No probability estimation is
needed. split_haplotype_multi() takes a single named
integer vector pop_codes of length K, matching the codes
used in the MSP file (RFMix default: 0-based integers).
pop_codes_phased <- c(AFR = 0L, EUR = 1L, NAT = 2L)
set.seed(2025)
N_VAR2 <- 8
N_SAMP2 <- 6
# Per-haplotype allele matrices (0/1)
gt_hap0 <- matrix(
sample(0:1, N_VAR2 * N_SAMP2, replace = TRUE, prob = c(0.9, 0.1)),
nrow = N_VAR2, ncol = N_SAMP2
)
gt_hap1 <- matrix(
sample(0:1, N_VAR2 * N_SAMP2, replace = TRUE, prob = c(0.9, 0.1)),
nrow = N_VAR2, ncol = N_SAMP2
)
# Per-haplotype ancestry codes (0=AFR, 1=EUR, 2=NAT)
anc_hap0 <- matrix(
sample(0:2, N_VAR2 * N_SAMP2, replace = TRUE, prob = c(0.4, 0.4, 0.2)),
nrow = N_VAR2, ncol = N_SAMP2
)
anc_hap1 <- matrix(
sample(0:2, N_VAR2 * N_SAMP2, replace = TRUE, prob = c(0.4, 0.4, 0.2)),
nrow = N_VAR2, ncol = N_SAMP2
)
phased_out <- split_haplotype_multi(gt_hap0, gt_hap1,
anc_hap0, anc_hap1,
pop_codes_phased)
cat("Output names:", names(phased_out), "\n")
#> Output names: AFR EUR NAT
cat("AFR dosage (first 4 variants):\n")
#> AFR dosage (first 4 variants):
print(phased_out$AFR[1:4, ])
#> [,1] [,2] [,3] [,4] [,5] [,6]
#> [1,] 0 0 0 0 0 0
#> [2,] 0 0 0 1 0 0
#> [3,] 0 0 0 0 0 0
#> [4,] 0 1 0 0 0 0
# Exact reconstruction: AFR + EUR + NAT == hap0 + hap1 everywhere
raw_diploid <- gt_hap0 + gt_hap1
recon_error <- max(abs((phased_out$AFR + phased_out$EUR + phased_out$NAT) -
raw_diploid))
cat("Max reconstruction error:", recon_error, "(must be 0)\n")
#> Max reconstruction error: 0 (must be 0)The phased function is exact because every haplotype allele goes to exactly one pool; there is no probabilistic component.
split_haplotype_multi() with K = 2 is numerically
identical to split_haplotype():
pop2 <- c(AFR = 0L, EUR = 1L)
gt0b <- matrix(c(1L,0L,1L,0L), 2, 2)
gt1b <- matrix(c(0L,1L,0L,1L), 2, 2)
a0b <- matrix(c(0L,1L,0L,1L), 2, 2)
a1b <- matrix(c(1L,0L,1L,0L), 2, 2)
m2 <- split_haplotype_multi(gt0b, gt1b, a0b, a1b, pop2)
p2 <- split_haplotype(gt0b, gt1b, a0b, a1b, pop2)
cat("AFR identical:", identical(unname(m2$AFR), p2$african), "\n")
#> AFR identical: TRUE
cat("EUR identical:", identical(unname(m2$EUR), p2$european), "\n")
#> EUR identical: TRUEThe same interface extends directly to four or five populations by
adding entries to pure_codes and rows to
mixed_codes. There are K(K−1)/2 mixed-pair codes for K
populations — 6 for K = 4, 10 for K = 5 — and mixed_codes
must list all of them.
# Four populations: AFR, EUR, NAT, EAS
pure_4 <- c(AFR = 3L, EUR = 1L, NAT = 4L, EAS = 5L)
mixed_4 <- data.frame(
code = c(2L, 6L, 7L, 8L, 9L, 10L),
pop1 = c("AFR","AFR","AFR","EUR","EUR","NAT"),
pop2 = c("EUR","NAT","EAS","NAT","EAS","EAS")
)
set.seed(11)
n_var4 <- 6; n_samp4 <- 16
gt4 <- matrix(sample(0:2, n_var4 * n_samp4, replace = TRUE,
prob = c(0.85, 0.12, 0.03)),
nrow = n_var4, ncol = n_samp4)
anc4 <- matrix(sample(c(pure_4, mixed_4$code), n_var4 * n_samp4, replace = TRUE),
nrow = n_var4, ncol = n_samp4)
out4 <- split_diploid_multi(gt4, anc4, pure_4, mixed_4)
cat("Populations:", paste(names(out4), collapse = ", "), "\n")
#> Populations: AFR, EUR, NAT, EAS
# Conservation check: pooled dosage never exceeds the raw genotype
total4 <- Reduce("+", out4)
cat("Max excess of pooled dosage over raw gt:", max(total4 - gt4), "(must be <= 0)\n")
#> Max excess of pooled dosage over raw gt: 0 (must be <= 0)
# Five populations: AFR, EUR, NAT, EAS, SAS — 10 mixed-pair codes
mixed_5 <- data.frame(
code = 7:16,
pop1 = c("AFR","AFR","AFR","AFR","EUR","EUR","EUR","NAT","NAT","EAS"),
pop2 = c("EUR","NAT","EAS","SAS","NAT","EAS","SAS","EAS","SAS","SAS")
)
cat("K=5 mixed_codes has", nrow(mixed_5), "rows (expected 5*4/2 =", 5 * 4 / 2, ")\n")
#> K=5 mixed_codes has 10 rows (expected 5*4/2 = 10 )Parts A and B each flagged an edge case and deferred it here: the pure singleton (D = 0, every alt allele carried by mixed-ancestry heterozygotes) is only the extreme end of a broader problem — the raw formula gets unreliable well before D actually hits zero, any time ambiguous mixed-het carriers dominate the count. Take 10 AFR/EUR mixed-het carriers (code 2) and just one pure-EUR hom-alt carrier (code 1) at the same variant — nothing else:
gt_dom <- matrix(c(rep(1L, 10), 2L), nrow = 1) # 10 mixed hets + 1 pure EUR hom-alt
anc_dom <- matrix(c(rep(2L, 10), 1L), nrow = 1) # code 2 = AFR/EUR, code 1 = EUR
r_raw <- split_diploid_multi(gt_dom, anc_dom, pure_codes, mixed_codes) # gla=NULL (default)
cat("Raw split -> AFR total:", sum(r_raw$AFR),
" EUR total:", sum(r_raw$EUR), " NAT total:", sum(r_raw$NAT), "\n")
#> Raw split -> AFR total: 0 EUR total: 12 NAT total: 0D > 0 here, so this isn’t the singleton branch — the raw formula happily computes p_AFR = 0, p_EUR = 1 from that one unambiguous carrier, and confidently assigns all ten ambiguous carriers’ alt alleles to EUR. That’s not a “no answer, guess 50/50” problem; it’s a confidently wrong-looking answer extrapolated from a single data point. Nothing in the raw formula distinguishes “1 unambiguous carrier informing 2 total alt alleles” from “1 unambiguous carrier informing 100” — it takes whatever ratio it computes at face value, regardless of how little evidence backs it.
GLA (global local ancestry) shrinkage fixes both
this case and the pure singleton with one mechanism: discount the raw
pairwise ratio in proportion to how much of a variant’s
evidence is actually ambiguous for that specific pair,
and fill in the rest from the chromosome arm’s independently-estimated
ancestry mixture, renormalised to just the two populations in the pair.
For the mixed pair (i,j) behind a given ambiguous het
(here, AFR/EUR):
Blending against the pair-renormalised GLA ratio — rather than each
population’s raw GLA_k[arm] directly — matters once K >
2: it guarantees an ambiguous AFR/EUR het never receives any dosage
toward a third population like NAT, however strongly shrinkage applies.
The example below confirms this: NAT total stays exactly 0
even though GLA_NAT[arm] = 0.05 > 0 — a formula that
blended each p_k toward its raw GLA_k[arm]
directly (instead of the pair-renormalised ratio above) would
incorrectly leak some dosage to NAT here. For K = 2, where
GLA_AFR[arm] + GLA_EUR[arm] = 1 by construction, the
pair-renormalised target collapses to plain GLA_AFR[arm]
and this formula reduces to the simpler two-population case from Part
A.
w_ij = 0 (pair (i,j)’s own unambiguous
carriers dominate) reduces to exactly the raw pairwise ratio from Parts
A/B. w_ij = 1 (the pure singleton for that pair) falls back
entirely to the pair’s GLA-conditioned ratio — not a coin flip; that’s
just one end of this continuum, not a separate mechanism.
GLA_k[arm] itself is computed once per chromosome arm
(split at the centromere) directly from the MSP ancestry
tracts, not from genotypes or phase, so it stays
well-estimated regardless of how rare any individual variant is.
split_diploid()/split_diploid_multi()
accept this as an optional gla (proportions) +
arm_id (which row of gla each variant uses)
argument pair; leaving them at their default NULL
reproduces Parts A/B exactly, as seen above. Revisiting the example
above with a gla argument supplied — say, this arm’s cohort
is mostly AFR/EUR admixed with little NAT ancestry:
# This arm's actual ancestry mixture (from MSP tracts, not this variant)
gla_arm <- matrix(c(0.70, 0.25, 0.05), nrow = 1,
dimnames = list("p", c("AFR", "EUR", "NAT")))
r_gla <- split_diploid_multi(gt_dom, anc_dom, pure_codes, mixed_codes,
gla = gla_arm, arm_id = 0L)
cat("Without GLA -> AFR total:", sum(r_raw$AFR),
" EUR total:", sum(r_raw$EUR), " NAT total:", sum(r_raw$NAT), "\n")
#> Without GLA -> AFR total: 0 EUR total: 12 NAT total: 0
cat("With GLA -> AFR total:", round(sum(r_gla$AFR), 2),
" EUR total:", round(sum(r_gla$EUR), 2),
" NAT total:", round(sum(r_gla$NAT), 2), "\n")
#> With GLA -> AFR total: 6.14 EUR total: 5.86 NAT total: 0Here w = 10/11 ≈ 0.91 — the near-total ambiguity is
almost fully discounted, and the ten mixed-het carriers now split
roughly 2:1 toward AFR (matching the arm’s 0.70/0.25 mixture) instead of
being handed to EUR outright. Same ten carriers, same one pure-EUR data
point — but the GLA-informed split reflects what’s actually known about
this arm, this cohort, instead of extrapolating from one
sample. The pure singleton from B.1 is the same mechanism at its extreme
(w = 1, no unambiguous carriers at all):
gt_s <- matrix(1L, nrow = 1, ncol = 1)
anc_s <- matrix(2L, nrow = 1, ncol = 1) # code 2 = AFR/EUR
r_flat <- split_diploid_multi(gt_s, anc_s, pure_codes, mixed_codes) # gla=NULL (default)
r_gla_s <- split_diploid_multi(gt_s, anc_s, pure_codes, mixed_codes,
gla = gla_arm, arm_id = 0L)
cat("Without GLA (flat fallback) -> AFR:", r_flat$AFR[1,1],
" EUR:", r_flat$EUR[1,1], " NAT:", r_flat$NAT[1,1], "\n")
#> Without GLA (flat fallback) -> AFR: 0.5 EUR: 0.5 NAT: 0
cat("With GLA (this arm's ancestry) -> AFR:", round(r_gla_s$AFR[1,1], 3),
" EUR:", round(r_gla_s$EUR[1,1], 3), " NAT:", round(r_gla_s$NAT[1,1], 3), "\n")
#> With GLA (this arm's ancestry) -> AFR: 0.737 EUR: 0.263 NAT: 0In practice you won’t build gla/arm_id by
hand — ancestry_split(mode = "dosage") computes both
automatically from the MSP file’s tracts and applies them for you. Pass
use_gla = FALSE to turn shrinkage off and get the flat
Parts A/B behavior instead, e.g. for comparing against the original
estimator:
split_gla <- ancestry_split(vcf_path, msp_path, mode = "dosage") # use_gla=TRUE (default)
split_flat <- ancestry_split(vcf_path, msp_path, mode = "dosage", use_gla = FALSE)See ?split_diploid_multi and
?split_haplotype_multi for full argument documentation.
Next: vignette("toy-vcf-msp-example") builds these same
inputs from real VCF + MSP file formats.
sessionInfo()
#> R version 4.6.1 (2026-06-24)
#> Platform: x86_64-pc-linux-gnu
#> Running under: Ubuntu 24.04.5 LTS
#>
#> Matrix products: default
#> BLAS: /usr/lib/x86_64-linux-gnu/openblas-pthread/libblas.so.3
#> LAPACK: /usr/lib/x86_64-linux-gnu/openblas-pthread/libopenblasp-r0.3.26.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: UTC
#> tzcode source: system (glibc)
#>
#> attached base packages:
#> [1] stats graphics grDevices utils datasets methods base
#>
#> other attached packages:
#> [1] lantern_0.99.0
#>
#> loaded via a namespace (and not attached):
#> [1] digest_0.6.39 desc_1.4.3 R6_2.6.1 fastmap_1.2.0
#> [5] xfun_0.61 cachem_1.1.0 knitr_1.52 htmltools_0.5.9
#> [9] rmarkdown_2.32 lifecycle_1.0.5 cli_3.6.6 sass_0.4.10
#> [13] pkgdown_2.2.1 textshaping_1.0.5 jquerylib_0.1.4 systemfonts_1.3.2
#> [17] compiler_4.6.1 tools_4.6.1 ragg_1.5.2 bslib_0.12.0
#> [21] evaluate_1.0.5 yaml_2.3.12 otel_0.2.0 jsonlite_2.0.0
#> [25] rlang_1.3.0 fs_2.1.0