--- title: "Toy example: splitting a VCF + MSP file (two and more populations)" output: rmarkdown::html_vignette vignette: > %\VignetteIndexEntry{Toy example: splitting a VCF + MSP file (two and more populations)} %\VignetteEngine{knitr::rmarkdown} %\VignetteEncoding{UTF-8} --- ```{r setup, include=FALSE} knitr::opts_chunk$set(collapse = TRUE, comment = "#>") library(lantern) ``` ## Overview `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. *** # Part A: Two populations ```{r toy-files} 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=', 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) ``` ```{r toy-split} 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") cat("=== phased (mode = \"haplotype\") ===\n") for (p in pop_names) { cat(p, "dosage:\n"); print(phased[[p]]) } cat("\n=== unphased (mode = \"dosage\") ===\n") for (p in pop_names) { cat(p, "dosage:\n"); print(unphased[[p]]) } cat("\nVariants kept:", nrow(phased$variant_info), "\n") cat("Samples: ", length(phased$sample_ids), "\n") # 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") ``` Each 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`. *** # Part B: More than two populations (K = 4) Nothing about the call changes — same function, same `mode` argument, same return shape, just more populations declared in the MSP header. ```{r k4-setup} 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=', 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) ``` ```{r k4-split} 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") cat("Variants kept:", nrow(phased4$variant_info), "\n") cat("ancestry_counts columns:", paste(colnames(phased4$ancestry_counts), collapse = ", "), "\n\n") cat("=== phased AFR (first 3 variants) ===\n"); print(phased4$AFR[1:3, ]) cat("=== unphased AFR (first 3 variants) ===\n"); print(round(unphased4$AFR[1:3, ], 2)) ``` Next: `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. ## Session information ```{r session} sessionInfo() ```