vignettes/toy-vcf-msp-example.Rmd
toy-vcf-msp-example.Rmdancestry_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 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
#> [4] fastmap_1.2.0 xfun_0.61 cachem_1.1.0
#> [7] knitr_1.52 htmltools_0.5.9 rmarkdown_2.32
#> [10] lifecycle_1.0.5 cli_3.6.6 sass_0.4.10
#> [13] pkgdown_2.2.1 data.table_1.18.6.1 textshaping_1.0.5
#> [16] jquerylib_0.1.4 systemfonts_1.3.2 compiler_4.6.1
#> [19] tools_4.6.1 ragg_1.5.2 bslib_0.12.0
#> [22] evaluate_1.0.5 yaml_2.3.12 otel_0.2.0
#> [25] jsonlite_2.0.0 rlang_1.3.0 fs_2.1.0