lantern splits each sample’s alternate-allele dosage at a variant into one dosage per ancestral population. There are two ways to do this:
| Mode | Uses | How an allele is assigned |
|---|---|---|
Phased ("haplotype") |
the allele and the ancestry of each haplotype | to the population of the haplotype that carries it |
Unphased ("dosage") |
the genotype (0/1/2) and the pair of ancestries at the site | shared by estimated proportions when the genotype is ambiguous |
Most analyses only need ancestry_split(), which reads a
VCF and an MSP file and returns the split dosages. Underneath it are
four matrix-level functions, exported so you can split matrices you
already have in memory, or see exactly what the split does, as this
vignette does:
| Unphased | Phased | |
|---|---|---|
| Any number of populations | split_diploid_multi() |
split_haplotype_multi() |
| Two populations only | split_diploid() |
split_haplotype() |
ancestry_split(mode = "dosage") calls
split_diploid_multi(), and
ancestry_split(mode = "haplotype") calls
split_haplotype_multi(). It does so for any number of
populations, including two.split_diploid() and split_haplotype() are
two-population shortcuts with fixed AFR/EUR codes that return
african and european matrices. With two
populations they give the same results as the _multi
versions (checked in Part B). ancestry_split() does not use
them.The unphased and phased functions take different inputs, all of them variants x samples matrices:
| Unphased functions | Phased functions | |
|---|---|---|
| Genotype | one matrix of 0/1/2 alternate-allele counts | two matrices of 0/1, one per haplotype |
| Ancestry | one matrix of codes for each sample’s pair of ancestries, e.g. AFR/EUR | two matrices of population codes, one per haplotype, as in the MSP file |
ancestry_split() builds both from the VCF and MSP files.
For the unphased split it adds the two haplotypes to get the genotype
and combines the two haplotype ancestries into a pair code, so you never
choose codes yourself.
One more difference: ancestry_split(mode = "dosage")
always shrinks the unphased split toward 1/2 (Part C). The unphased
functions shrink only with shrink = TRUE. Parts A and B use
the default shrink = FALSE to show the plain estimate, so
their numbers differ slightly from what ancestry_split()
would give on the same data.
Part A covers two populations, Part B any number of populations, and
Part C the shrinkage. vignette("toy-vcf-msp-example") shows
the same splits starting from VCF and MSP files.
split_diploid()Without phase we only know, for each sample at a variant, its
genotype (0, 1 or 2 alternate alleles) and its pair of ancestries.
split_diploid() expects these codes:
(split_diploid_multi() lets you choose the codes; see
B.1.)
For most samples the split is obvious. An AFR/AFR sample’s alleles are all African, and an AFR/EUR sample with genotype 2 has one allele on each haplotype, so one African and one European. The only ambiguous case is an AFR/EUR heterozygote: it has one alternate allele, and we cannot tell which haplotype carries it. The full rule:
| Ancestry | Genotype | African dosage | European dosage |
|---|---|---|---|
| AFR/AFR | 2 | 2 | 0 |
| AFR/AFR | 1 | 1 | 0 |
| AFR/EUR | 2 | 1 | 1 |
| AFR/EUR | 1 | \(P_\text{AFR}\) | \(P_\text{EUR}\) |
| EUR/EUR | 2 | 0 | 2 |
| EUR/EUR | 1 | 0 | 1 |
Genotype 0 contributes nothing to either dosage. Only the highlighted row needs the proportions \(P_\text{AFR}\) and \(P_\text{EUR} = 1 - P_\text{AFR}\), which are estimated from the other rows as shown below.
Consider a single variant in six samples:
gt <- matrix(c(1, 2, 1, 1, 0, 2), nrow = 1,
dimnames = list("var1", paste0("S", 1:6)))
anc <- matrix(c(3, 3, 1, 2, 2, 2), nrow = 1,
dimnames = list("var1", paste0("S", 1:6)))
rbind(genotype = gt[1, ], ancestry = anc[1, ])
#> S1 S2 S3 S4 S5 S6
#> genotype 1 2 1 1 0 2
#> ancestry 3 3 1 2 2 2S1 and S2 are AFR/AFR carriers, S3 is a EUR/EUR carrier, S6 is an AFR/EUR homozygote and S4 is the one ambiguous AFR/EUR heterozygote. To split S4’s allele we estimate how this variant’s alternate alleles divide between African and European haplotypes, using only the unambiguous carriers. Write \(N_{kk}^{(g)}\) for the number of samples with ancestry \(k/k\) and genotype \(g\), and \(N_{\text{AFR,EUR}}^{(g)}\) for those with ancestry AFR/EUR. The unambiguous alternate alleles on each ancestry number
\[ n_\text{AFR} = 2N_\text{AFR,AFR}^{(2)} + N_\text{AFR,AFR}^{(1)} + N_\text{AFR,EUR}^{(2)}, \qquad n_\text{EUR} = 2N_\text{EUR,EUR}^{(2)} + N_\text{EUR,EUR}^{(1)} + N_\text{AFR,EUR}^{(2)}, \]
since an AFR/EUR homozygote contributes one allele to each ancestry. Their total \(D = n_\text{AFR} + n_\text{EUR}\) is all alternate alleles at the variant minus those of the ambiguous heterozygotes, and the proportions are
\[ P_\text{AFR} = \frac{n_\text{AFR}}{D}, \qquad P_\text{EUR} = \frac{n_\text{EUR}}{D}. \]
The ambiguous heterozygote then gets \(P_\text{AFR}\) of its allele as African dosage and \(P_\text{EUR}\) as European dosage. For our variant:
N <- function(a, g) sum(anc == a & gt == g) # codes: 3 = AFR/AFR, 2 = AFR/EUR, 1 = EUR/EUR
n_afr <- 2 * N(3, 2) + N(3, 1) + N(2, 2)
n_eur <- 2 * N(1, 2) + N(1, 1) + N(2, 2)
D <- n_afr + n_eur
c(n_AFR = n_afr, n_EUR = n_eur, D = D, P_AFR = n_afr / D, P_EUR = n_eur / D)
#> n_AFR n_EUR D P_AFR P_EUR
#> 4.0000000 2.0000000 6.0000000 0.6666667 0.3333333split_diploid() returns the two dosage matrices. S4’s
allele is split 2/3 African and 1/3 European, and every other sample’s
split is whole:
res <- split_diploid(gt, anc)
split_table <- rbind(AFR = res$african[1, ], EUR = res$european[1, ])
colnames(split_table) <- colnames(gt)
split_table
#> S1 S2 S3 S4 S5 S6
#> AFR 1 2 0 0.6666667 0 1
#> EUR 0 0 1 0.3333333 0 1If every carrier of a variant is an ambiguous heterozygote, \(D = 0\) and the proportions cannot be
estimated. split_diploid() then splits each allele 1/2 and
1/2:
gt_edge <- matrix(c(1, 0, 1), nrow = 1)
anc_edge <- matrix(c(2, 1, 2), nrow = 1)
split_diploid(gt_edge, anc_edge)
#> $african
#> [,1] [,2] [,3]
#> [1,] 0.5 0 0.5
#>
#> $european
#> [,1] [,2] [,3]
#> [1,] 0.5 0 0.5Part C shows that estimates from very few unambiguous carriers are
also unreliable, and how ancestry_split() handles that.
split_haplotype()With phased genotypes and local ancestry for each haplotype, nothing
needs to be estimated: each haplotype’s allele goes to that haplotype’s
population. split_haplotype() takes four matrices of the
same shape (variants x samples): the alleles (0/1) on haplotype 0 and on
haplotype 1, and the ancestry of each haplotype, by default with RFMix’s
codes 0 = AFR and 1 = EUR (set pop_codes if your MSP file
numbers them differently).
ids <- list(paste0("var", 1:3), paste0("S", 1:3))
gt0 <- matrix(c(1, 0, 1, 0, 1, 0, 1, 1, 0), nrow = 3, dimnames = ids)
gt1 <- matrix(c(0, 1, 0, 1, 0, 1, 0, 1, 1), nrow = 3, dimnames = ids)
a0 <- matrix(c(0, 0, 1, 1, 0, 0, 0, 1, 0), nrow = 3, dimnames = ids)
a1 <- matrix(c(1, 1, 0, 0, 1, 1, 1, 0, 1), nrow = 3, dimnames = ids)Take S1 at var1: haplotype 0 carries the allele
(gt0 is 1) and is African (a0 is 0), and
haplotype 1 does not carry it. So S1 has an African dosage of 1 and a
European dosage of 0 at var1.
phased <- split_haplotype(gt0, gt1, a0, a1)
phased$african
#> [,1] [,2] [,3]
#> [1,] 1 1 1
#> [2,] 0 1 1
#> [3,] 0 0 0
phased$european
#> [,1] [,2] [,3]
#> [1,] 0 0 0
#> [2,] 1 0 1
#> [3,] 1 1 1The two dosages always add up to the genotype,
gt0 + gt1:
split_diploid_multi() and
split_haplotype_multi() do the same for any number of
populations K. They are the functions ancestry_split()
calls. With two populations they give the same results as the functions
in Part A; the only difference is that their results are named by
population (AFR, EUR, …) rather than
african and european.
split_diploid_multi()With K populations a sample’s pair of ancestries is either pure (both haplotypes from the same population) or mixed (from two different populations). There are K pure and K(K-1)/2 mixed combinations:
| K | Pure | Mixed | Total codes |
|---|---|---|---|
| 2 | 2 | 1 | 3 |
| 3 | 3 | 3 | 6 |
| 4 | 4 | 6 | 10 |
You choose a code for each combination and pass two descriptions to
the function: pure_codes, the code of each pure
combination, and mixed_codes, a table with the code and the
two populations of each mixed combination. For three populations
(African, European and Native American):
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"))
mixed_codes
#> code pop1 pop2
#> 1 2 AFR EUR
#> 2 5 AFR NAT
#> 3 6 EUR NATWe simulate 12 rare variants in 30 samples from a three-way admixed group:
set.seed(2025)
gt3 <- matrix(sample(0:2, 12 * 30, replace = TRUE, prob = c(0.88, 0.10, 0.02)),
nrow = 12, ncol = 30)
anc3 <- matrix(sample(c(pure_codes, mixed_codes$code), 12 * 30, replace = TRUE,
prob = c(0.30, 0.20, 0.10, 0.20, 0.10, 0.10)),
nrow = 12, ncol = 30)
out <- split_diploid_multi(gt3, anc3, pure_codes, mixed_codes)
names(out)
#> [1] "AFR" "EUR" "NAT"The estimate works as in A.1, with one proportion per population. Pure carriers count for their population and mixed homozygotes contribute one allele to each of their two populations. For a variant, let \(N_{kk}^{(g)}\) count carriers with homozygous local ancestry \(k/k\) and alternate-allele dosage \(g\), and \(N_{ij}^{(g)} = N_{ji}^{(g)}\) those with heterozygous local ancestry \(i/j\), \(i \ne j\). The number of unambiguous alternate alleles on ancestry \(k\) is
\[ n_k = 2N_{kk}^{(2)} + N_{kk}^{(1)} + \sum_{i < j:\; k \in \{i,j\}} N_{ij}^{(2)}, \]
and their total \(D = \sum_k n_k\) is the \(D\) of A.1. An ambiguous carrier with local ancestry \(i/j\) is split within its own pair. Without shrinkage (Part C) it gets
\[ \tilde{\pi}_{ij} = \frac{n_i}{n_i + n_j} \]
of its allele in ancestry \(i\) and \(1 - \tilde{\pi}_{ij}\) in ancestry \(j\).
To check this, we take the first variant with an ambiguous AFR/EUR heterozygote and count the unambiguous alleles by hand:
v <- which(rowSums(anc3 == 2 & gt3 == 1) > 0)[1]
g <- gt3[v, ]
a <- anc3[v, ]
N <- function(code, genotype) sum(a == code & g == genotype)
# codes: 3 = AFR/AFR, 1 = EUR/EUR, 4 = NAT/NAT, 2 = AFR/EUR, 5 = AFR/NAT, 6 = EUR/NAT
n_k <- c(AFR = 2 * N(3, 2) + N(3, 1) + N(2, 2) + N(5, 2),
EUR = 2 * N(1, 2) + N(1, 1) + N(2, 2) + N(6, 2),
NAT = 2 * N(4, 2) + N(4, 1) + N(5, 2) + N(6, 2))
n_k
#> AFR EUR NAT
#> 2 1 2The heterozygote’s African share should be \(\tilde{\pi}_\text{AFR,EUR} = n_\text{AFR} / (n_\text{AFR} + n_\text{EUR})\), which matches the function’s output:
s <- which(a == 2 & g == 1)[1]
c(by_hand = unname(n_k["AFR"] / (n_k["AFR"] + n_k["EUR"])),
lantern = out$AFR[v, s])
#> by_hand lantern
#> 0.6666667 0.6666667With its default codes, which are the two-population codes of A.1,
split_diploid_multi() gives the same split as
split_diploid() on the variant from A.1:
split_haplotype_multi()The phased split needs only the code of each population, as used in the MSP file (RFMix numbers populations from 0):
set.seed(2025)
hap0 <- matrix(rbinom(8 * 6, 1, 0.1), nrow = 8)
hap1 <- matrix(rbinom(8 * 6, 1, 0.1), nrow = 8)
anc0 <- matrix(sample(0:2, 8 * 6, replace = TRUE, prob = c(0.4, 0.4, 0.2)), nrow = 8)
anc1 <- matrix(sample(0:2, 8 * 6, replace = TRUE, prob = c(0.4, 0.4, 0.2)), nrow = 8)
phased3 <- split_haplotype_multi(hap0, hap1, anc0, anc1,
pop_codes = c(AFR = 0L, EUR = 1L, NAT = 2L))
phased3$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 0Each allele goes to exactly one population, so the three dosages add up to the genotype:
With two populations, split_haplotype_multi() gives the
same result as split_haplotype():
More populations only means more entries in pure_codes
and more rows in mixed_codes, one for every pair of
populations. For four populations:
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)
gt4 <- matrix(sample(0:2, 6 * 16, replace = TRUE, prob = c(0.85, 0.12, 0.03)),
nrow = 6)
anc4 <- matrix(sample(c(pure_4, mixed_4$code), 6 * 16, replace = TRUE), nrow = 6)
out4 <- split_diploid_multi(gt4, anc4, pure_4, mixed_4)
names(out4)
#> [1] "AFR" "EUR" "NAT" "EAS"The proportions in A.1 and B.1 are estimated from the unambiguous carriers of each variant. Rare variants have few carriers, and an estimate from one or two of them can be far from the truth. Consider a variant carried by ten ambiguous AFR/EUR heterozygotes and a single EUR/EUR homozygote:
gt_dom <- matrix(c(rep(1L, 10), 2L), nrow = 1)
anc_dom <- matrix(c(rep(2L, 10), 1L), nrow = 1)
raw <- split_diploid_multi(gt_dom, anc_dom, pure_codes, mixed_codes)
c(AFR = sum(raw$AFR), EUR = sum(raw$EUR), NAT = sum(raw$NAT))
#> AFR EUR NAT
#> 0 12 0The only unambiguous carrier is European, so the estimate is \(P_\text{EUR} = 1\), and all ten ambiguous alleles are assigned to European haplotypes on the strength of one sample.
ancestry_split() therefore pulls the raw split \(\tilde{\pi}_{ij}\) toward 1/2. How far
depends on \(w_{ij}\), the share of the
evidence about the pair that is ambiguous. Only the alleles on
ancestries \(i\) and \(j\) (\(n_i\) and \(n_j\) from B.1), and the pair’s own
ambiguous carriers, contribute to that split. So the raw split and the
pair’s ambiguous share are
\[ \tilde{\pi}_{ij} = \frac{n_i}{n_i + n_j}, \qquad w_{ij} = \frac{N_{ij}^{(1)}}{n_i + n_j + N_{ij}^{(1)}}, \]
and the shrunk split is
\[ \pi_{ij} = (1 - w_{ij})\,\tilde{\pi}_{ij} + 1/2\; w_{ij} = \frac{n_i + 1/2\, N_{ij}^{(1)}}{n_i + n_j + N_{ij}^{(1)}}, \]
i.e. each ambiguous carrier counts as half an allele on each of its two ancestries. When the pair has no unambiguous allele (\(n_i + n_j = 0\), including \(D = 0\)), \(w_{ij} = 1\) and \(\pi_{ij} = 1/2\), the same as the fallback in A.1. Alleles on other ancestries, and other pairs’ ambiguous carriers, do not enter \(w_{ij}\). The carrier gets \(x_i = \pi_{ij}\) and \(x_j = 1 - \pi_{ij}\), and 0 for every other ancestry. When most of the pair’s evidence is unambiguous, \(w_{ij}\) is close to 0 and \(\tilde{\pi}_{ij}\) is used almost unchanged.
The target is 1/2 because an AFR/EUR heterozygote has exactly one African and one European haplotype at the site. Which of the two carries the allele depends on how common the allele is on African versus European haplotypes, not on how much African ancestry the cohort has overall. An earlier version of lantern shrank toward the cohort’s average ancestry (about 0.82 African on chromosome 19 in our data). In simulations based on that data, ambiguous alleles were in fact on the African haplotype half of the time, and shrinking toward 1/2 assigned them more accurately.
split_diploid() and split_diploid_multi()
apply the same shrinkage with shrink = TRUE:
shrunk <- split_diploid_multi(gt_dom, anc_dom, pure_codes, mixed_codes,
shrink = TRUE)
rbind(raw = c(AFR = sum(raw$AFR), EUR = sum(raw$EUR)),
shrunk = c(AFR = sum(shrunk$AFR), EUR = sum(shrunk$EUR)))
#> AFR EUR
#> raw 0.000000 12.000000
#> shrunk 4.166667 7.833333Here the pair’s evidence is the EUR/EUR homozygote’s two alleles (\(n_\text{AFR} = 0\), \(n_\text{EUR} = 2\)) and the ten ambiguous heterozygotes, so \(w_\text{AFR,EUR} = 10/(0 + 2 + 10) = 10/12\). Each heterozygote’s African share becomes \(\pi_\text{AFR,EUR} = (0 + 10/2) / 12 = 5/12 \approx 0.42\) instead of 0:
Adding a hundred Native American carriers changes nothing, because they do not enter \(w_\text{AFR,EUR}\):
gt_nat <- cbind(gt_dom, matrix(1L, nrow = 1, ncol = 100))
anc_nat <- cbind(anc_dom, matrix(pure_codes[["NAT"]], nrow = 1, ncol = 100))
split_diploid_multi(gt_nat, anc_nat, pure_codes, mixed_codes,
shrink = TRUE)$AFR[1, 1]
#> [1] 0.4166667The variant from A.1 has much more unambiguous evidence: \(n_\text{AFR} = 4\) (S1, S2’s two alleles, and one from S6), \(n_\text{EUR} = 2\) (S3 and one from S6) and one ambiguous heterozygote, so \(w_\text{AFR,EUR} = 1/(4 + 2 + 1) = 1/7\). S4’s African share moves only slightly, from \(\tilde{\pi}_\text{AFR,EUR} = 2/3\) to \(\pi_\text{AFR,EUR} = (4 + 1/2) / 7 = 9/14 \approx 0.64\):
You do not need to set this yourself:
ancestry_split(mode = "dosage") always shrinks.
sessionInfo()
#> R version 4.6.2 (2026-10-29)
#> Platform: x86_64-pc-linux-gnu
#> Running under: Ubuntu 26.04.1 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.32.so; LAPACK version 3.12.0
#>
#> locale:
#> [1] LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C
#> [3] LC_TIME=en_US.UTF-8 LC_COLLATE=en_US.UTF-8
#> [5] LC_MONETARY=en_US.UTF-8 LC_MESSAGES=en_US.UTF-8
#> [7] LC_PAPER=en_US.UTF-8 LC_NAME=C
#> [9] LC_ADDRESS=C LC_TELEPHONE=C
#> [11] LC_MEASUREMENT=en_US.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] rmarkdown_2.32 lantern_0.99.1 BiocStyle_2.41.0
#>
#> loaded via a namespace (and not attached):
#> [1] tidyselect_1.2.1 dplyr_1.2.1 Biostrings_2.81.9
#> [4] gdsfmt_1.49.12 fastmap_1.2.0 CompQuadForm_1.4.4
#> [7] digest_0.6.39 rpart_4.1.27 lifecycle_1.0.5
#> [10] survival_3.8-12 magrittr_2.0.5 compiler_4.6.2
#> [13] rlang_1.3.0 sass_0.4.10 tools_4.6.2
#> [16] yaml_2.3.12 data.table_1.18.6.1 knitr_1.52
#> [19] SeqArray_1.53.6 purrr_1.2.2 BiocGenerics_0.59.12
#> [22] sys_3.4.3 GWASExactHW_1.2 nnet_7.3-21
#> [25] grid_4.6.2 stats4_4.6.2 jomo_2.7-6
#> [28] GMMAT_1.5.0 mice_3.19.0 iterators_1.0.14
#> [31] MASS_7.3-66 cli_3.6.6 crayon_1.5.3
#> [34] reformulas_0.4.4 generics_0.1.4 otel_0.2.0
#> [37] minqa_1.2.8 cachem_1.1.0 operator.tools_1.6.3.1
#> [40] splines_4.6.2 parallel_4.6.2 BiocManager_1.30.27
#> [43] XVector_0.53.0 vctrs_0.7.3 boot_1.3-32
#> [46] glmnet_5.2 Matrix_1.7-6 jsonlite_2.0.0
#> [49] IRanges_2.47.5 S4Vectors_0.51.10 mitml_0.4-5
#> [52] maketools_1.3.2 foreach_1.5.2 tidyr_1.3.2
#> [55] jquerylib_0.1.4 glue_1.8.1 nloptr_2.2.1
#> [58] pan_2.0 codetools_0.2-20 shape_1.4.6.1
#> [61] GenomicRanges_1.65.4 lme4_2.1-0 tibble_3.3.1
#> [64] pillar_1.11.1 htmltools_0.5.9 Seqinfo_1.3.2
#> [67] R6_2.6.1 SeqVarTools_1.51.3 Rdpack_2.6.6
#> [70] formula.tools_1.7.1 evaluate_1.0.5 lattice_0.23-1
#> [73] Biobase_2.73.2 rbibutils_2.4.1 backports_1.5.1
#> [76] broom_1.0.13 bslib_0.12.0 Rcpp_1.1.2
#> [79] nlme_3.1-172 mgcv_1.9-5 logistf_1.26.1
#> [82] xfun_0.62 buildtools_1.0.0 pkgconfig_2.0.3