ancestry_split() (Step 1) is LANTERN’s entry point once
you have a phased VCF and an RFMix MSP file: it parses both, intersects
samples, maps variants to ancestry tracts, and runs either split via a
single mode argument ("haplotype" or
"dosage") — see vignette("split-intuition")
for what each mode computes on the underlying matrices.
Population count and names come directly from the MSP header, so the same call works unchanged for K = 2 or K > 2. Part A builds a two-population toy VCF + MSP; Part B does the same for four populations.
tmp_dir <- tempfile()
dir.create(tmp_dir)
vcf_path <- file.path(tmp_dir, "phased.vcf")
msp_path <- file.path(tmp_dir, "ancestry.msp.tsv")
# 3 variants × 3 samples, phased GTs
writeLines(c(
"##fileformat=VCFv4.1",
'##FORMAT=<ID=GT,Number=1,Type=String,Description="Genotype">',
paste(c("#CHROM","POS","ID","REF","ALT","QUAL","FILTER","INFO",
"FORMAT","S1","S2","S3"), collapse = "\t"),
paste(c("19","1100",".","A","G","100","PASS",".","GT","0|1","1|0","0|0"), collapse = "\t"),
paste(c("19","1200",".","C","T","100","PASS",".","GT","1|1","0|1","1|0"), collapse = "\t"),
paste(c("19","2100",".","G","A","100","PASS",".","GT","0|0","1|0","0|1"), collapse = "\t")
), vcf_path)
# RFMix MSP: 2 ancestry tracts, pop codes AFR=0 EUR=1
# tract 1 [1000,2000): S1=AFR|EUR S2=EUR|AFR S3=AFR|EUR
# tract 2 [2000,3000): S1=EUR|AFR S2=AFR|EUR S3=EUR|AFR
writeLines(c(
paste(c("#Subpopulation order/codes: AFR=0", "EUR=1"), collapse = "\t"),
paste(c("#chm","spos","epos","sgpos","egpos","n snps",
"S1.0","S1.1","S2.0","S2.1","S3.0","S3.1"), collapse = "\t"),
paste(c("19","1000","2000","0.0","0.1","2","0","1","1","0","0","1"), collapse = "\t"),
paste(c("19","2000","3000","0.1","0.2","1","1","0","0","1","1","0"), collapse = "\t")
), msp_path)phased <- ancestry_split(vcf_path, msp_path, mode = "haplotype", verbose = FALSE)
unphased <- ancestry_split(vcf_path, msp_path, mode = "dosage", verbose = FALSE)
pop_names <- setdiff(names(phased),
c("variant_info", "sample_ids", "mode",
"ancestry_counts", "tract_info", "overlap"))
cat("Populations:", paste(pop_names, collapse = ", "), "\n\n")
#> Populations: AFR, EUR
cat("=== phased (mode = \"haplotype\") ===\n")
#> === phased (mode = "haplotype") ===
for (p in pop_names) { cat(p, "dosage:\n"); print(phased[[p]]) }
#> AFR dosage:
#> [,1] [,2] [,3]
#> [1,] 0 0 0
#> [2,] 1 1 1
#> [3,] 0 1 1
#> EUR dosage:
#> [,1] [,2] [,3]
#> [1,] 1 1 0
#> [2,] 1 0 0
#> [3,] 0 0 0
cat("\n=== unphased (mode = \"dosage\") ===\n")
#>
#> === unphased (mode = "dosage") ===
for (p in pop_names) { cat(p, "dosage:\n"); print(unphased[[p]]) }
#> AFR dosage:
#> [,1] [,2] [,3]
#> [1,] 0.5 0.5 0.0
#> [2,] 1.0 0.5 0.5
#> [3,] 0.0 0.5 0.5
#> EUR dosage:
#> [,1] [,2] [,3]
#> [1,] 0.5 0.5 0.0
#> [2,] 1.0 0.5 0.5
#> [3,] 0.0 0.5 0.5
cat("\nVariants kept:", nrow(phased$variant_info), "\n")
#>
#> Variants kept: 3
cat("Samples: ", length(phased$sample_ids), "\n")
#> Samples: 3
# Sanity check: phased AFR + EUR reconstructs the raw diploid genotype exactly
vcf_lines <- readLines(vcf_path)[-(1:3)]
raw <- do.call(rbind, lapply(vcf_lines, function(l) {
gts <- strsplit(l, "\t")[[1]][10:12]
vapply(gts, function(g) sum(as.integer(strsplit(g, "|", fixed = TRUE)[[1]])), integer(1))
}))
max_err <- max(abs((phased$AFR + phased$EUR) - raw))
cat("Max phased reconstruction error (must be 0):", max_err, "\n")
#> Max phased reconstruction error (must be 0): 0Each call returns a named list: one dosage matrix per population
(named from the MSP header) plus variant_info,
sample_ids, mode,
ancestry_counts, and overlap. Need both splits
on a large cohort routinely? Cache each result (e.g. to RDS) rather than
re-parsing the VCF/MSP on every run — ancestry_split()
always re-parses from scratch.
The backward-compatible ancestry_split_dosage() (matrix
or VCF/MSP input) and ancestry_split_phased() (VCF/MSP
input, optionally also writing GDS) wrap the same logic for existing
code — see ?ancestry_split_dosage and
?ancestry_split_phased.
Nothing about the call changes — same function, same
mode argument, same return shape, just more populations
declared in the MSP header.
tmp_dir4 <- tempfile()
dir.create(tmp_dir4)
vcf4_path <- file.path(tmp_dir4, "phased4.vcf")
msp4_path <- file.path(tmp_dir4, "ancestry4.msp.tsv")
set.seed(99)
n_samp_k4 <- 24
samp_k4 <- paste0("S", seq_len(n_samp_k4))
# 10 variants, phased GT, mostly rare
gt_lines <- vapply(seq_len(10), function(i) {
gts <- sample(c("0|0","0|1","1|0","1|1"), n_samp_k4, replace = TRUE,
prob = c(0.85, 0.06, 0.06, 0.03))
paste(c("19", 1000L + i * 100L, ".", "A", "T", "100", "PASS", ".", "GT", gts),
collapse = "\t")
}, character(1))
writeLines(c(
"##fileformat=VCFv4.1",
'##FORMAT=<ID=GT,Number=1,Type=String,Description="Genotype">',
paste(c("#CHROM","POS","ID","REF","ALT","QUAL","FILTER","INFO","FORMAT", samp_k4),
collapse = "\t"),
gt_lines
), vcf4_path)
# Single tract spanning all variants; each sample gets a random pair of
# the four ancestries per haplotype
pop_codes_k4 <- c(AFR = 0L, EUR = 1L, NAT = 2L, EAS = 3L)
hap_cols <- unlist(lapply(samp_k4, function(s) paste0(s, c(".0", ".1"))))
hap_anc <- sample(pop_codes_k4, length(hap_cols), replace = TRUE)
writeLines(c(
paste(c("#Subpopulation order/codes: AFR=0", "EUR=1", "NAT=2", "EAS=3"),
collapse = "\t"),
paste(c("#chm","spos","epos","sgpos","egpos","n snps", hap_cols), collapse = "\t"),
paste(c("19", "1000", "3000", "0.0", "0.5", "10", hap_anc), collapse = "\t")
), msp4_path)phased4 <- ancestry_split(vcf4_path, msp4_path, mode = "haplotype",
chrom = "19", verbose = FALSE)
unphased4 <- ancestry_split(vcf4_path, msp4_path, mode = "dosage",
chrom = "19", verbose = FALSE)
pop4 <- setdiff(names(phased4),
c("variant_info", "sample_ids", "mode",
"ancestry_counts", "tract_info", "overlap"))
cat("Populations recovered from the MSP header:", paste(pop4, collapse = ", "), "\n")
#> Populations recovered from the MSP header: AFR, EUR, NAT, EAS
cat("Variants kept:", nrow(phased4$variant_info), "\n")
#> Variants kept: 9
cat("ancestry_counts columns:", paste(colnames(phased4$ancestry_counts), collapse = ", "), "\n\n")
#> ancestry_counts columns: AFR, EUR, NAT, EAS
cat("=== phased AFR (first 3 variants) ===\n"); print(phased4$AFR[1:3, ])
#> === phased AFR (first 3 variants) ===
#> [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11] [,12] [,13] [,14]
#> [1,] 0 0 0 0 0 0 0 0 0 0 0 0 0 0
#> [2,] 0 0 0 0 0 0 0 0 0 0 0 0 0 0
#> [3,] 0 0 0 0 0 0 0 0 0 0 0 0 0 0
#> [,15] [,16] [,17] [,18] [,19] [,20] [,21] [,22] [,23] [,24]
#> [1,] 0 0 0 0 0 0 0 0 0 0
#> [2,] 0 1 0 0 0 0 0 1 0 0
#> [3,] 0 0 0 0 0 0 0 0 0 0
cat("=== unphased AFR (first 3 variants) ===\n"); print(round(unphased4$AFR[1:3, ], 2))
#> === unphased AFR (first 3 variants) ===
#> [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11] [,12] [,13] [,14]
#> [1,] 0 0 0 0 0 0 0 0 0 0 0.00 0 0 0
#> [2,] 0 0 0 0 0 0 0 0 0 0 0.00 0 0 0
#> [3,] 0 0 0 0 0 0 0 0 0 0 0.37 0 0 0
#> [,15] [,16] [,17] [,18] [,19] [,20] [,21] [,22] [,23] [,24]
#> [1,] 0 0 0 0 0 0.00 0 0.00 0 0
#> [2,] 0 1 0 0 0 0.86 0 0.88 0 0
#> [3,] 0 0 0 0 0 0.05 0 0.00 0 0Next: vignette("real-data-chr19") takes a split result
like this one through write_ancestry_gds() (Step 2) and
ancestry_smmat() (Step 3) on real genotype data.
sessionInfo()
#> R version 4.6.1 (2026-06-24)
#> 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.0 BiocStyle_2.41.0
#>
#> loaded via a namespace (and not attached):
#> [1] shape_1.4.6.1 xfun_0.61 bslib_0.12.0
#> [4] formula.tools_1.7.1 SeqArray_1.53.4 Biobase_2.73.2
#> [7] CompQuadForm_1.4.4 lattice_0.23-1 Rdpack_2.6.6
#> [10] vctrs_0.7.3 tools_4.6.1 generics_0.1.4
#> [13] stats4_4.6.1 parallel_4.6.1 tibble_3.3.1
#> [16] pan_2.0 jomo_2.7-6 pkgconfig_2.0.3
#> [19] Matrix_1.7-6 data.table_1.18.6.1 S4Vectors_0.51.10
#> [22] lifecycle_1.0.5 compiler_4.6.1 Biostrings_2.81.9
#> [25] Seqinfo_1.3.2 codetools_0.2-20 GWASExactHW_1.2
#> [28] htmltools_0.5.9 sys_3.4.3 buildtools_1.0.0
#> [31] sass_0.4.10 yaml_2.3.12 glmnet_5.1
#> [34] mice_3.19.0 nloptr_2.2.1 pillar_1.11.1
#> [37] crayon_1.5.3 jquerylib_0.1.4 tidyr_1.3.2
#> [40] MASS_7.3-66 cachem_1.1.0 reformulas_0.4.4
#> [43] iterators_1.0.14 rpart_4.1.27 boot_1.3-32
#> [46] mitml_0.4-5 GMMAT_1.5.0 foreach_1.5.2
#> [49] nlme_3.1-171 tidyselect_1.2.1 digest_0.6.39
#> [52] dplyr_1.2.1 purrr_1.2.2 maketools_1.3.2
#> [55] splines_4.6.1 operator.tools_1.6.3.1 fastmap_1.2.0
#> [58] grid_4.6.1 cli_3.6.6 magrittr_2.0.5
#> [61] survival_3.8-12 broom_1.0.13 backports_1.5.1
#> [64] logistf_1.26.1 XVector_0.53.0 SeqVarTools_1.51.3
#> [67] nnet_7.3-21 lme4_2.0-6 otel_0.2.0
#> [70] evaluate_1.0.5 knitr_1.52 rbibutils_2.4.1
#> [73] GenomicRanges_1.65.4 IRanges_2.47.5 mgcv_1.9-4
#> [76] rlang_1.3.0 Rcpp_1.1.2 glue_1.8.1
#> [79] BiocManager_1.30.27 BiocGenerics_0.59.12 gdsfmt_1.49.8
#> [82] minqa_1.2.8 jsonlite_2.0.0 R6_2.6.1