Overview

ancestry_split() takes two files: a phased VCF with the genotypes and an RFMix (Maples et al. 2013) .msp.tsv file with the local ancestry. This vignette builds both files by hand for a tiny example, small enough that every number in the output can be checked against the input. Part A uses two populations. Part B uses four, to show that nothing about the call changes when there are more populations.

vignette("lantern") runs the full analysis on a larger simulated dataset, and vignette("split-intuition") explains the arithmetic of the split.

Part A: Two populations

The input files

We write the files to a temporary directory. The VCF has three variants and three samples (S1, S2, S3). The genotypes are phased: in 0|1 the first number is the allele on the first copy of the chromosome (haplotype 0) and the second number the allele on the second copy (haplotype 1) per VCF file standard.

tmp_dir  <- tempfile()
dir.create(tmp_dir)
vcf_path <- file.path(tmp_dir, "phased.vcf")

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)

The MSP file divides the chromosome into segments (tracts) and gives, for each sample and haplotype, the population that segment came from. Its first line maps population names to the codes used in the file (AFR = 0, EUR = 1). The ancestry columns are named <sample>.0 and <sample>.1 for the two haplotypes. Our file has two tracts:

Tract	Positions	S1 (hap 0, hap 1)	S2	S3
1	1000 - 2000	AFR, EUR	EUR, AFR	AFR, EUR
2	2000 - 3000	EUR, AFR	AFR, EUR	EUR, AFR

Variants at positions 1100 and 1200 fall in tract 1, and the variant at 2100 in tract 2.

msp_path <- file.path(tmp_dir, "ancestry.msp.tsv")
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)

Splitting with phase (`mode = "haplotype"`)

Because the genotypes are phased, we know which haplotype carries each allele, and the MSP file tells us that haplotype’s population. Each allele is simply counted for that population.

phased <- ancestry_split(vcf_path, msp_path, mode = "haplotype", verbose = FALSE)
names(phased)
#> [1] "AFR"             "EUR"             "variant_info"    "sample_ids"     
#> [5] "mode"            "ancestry_counts" "tract_info"      "overlap"

The result holds one dosage matrix per population, named from the MSP header (AFR, EUR), plus information about the variants and samples. The matrices have one row per variant and one column per sample. To make them easier to read we label them with positions and sample IDs:

label <- function(m, split) {
  dimnames(m) <- list(split$variant_info$pos, split$sample_ids)
  m
}
label(phased$AFR, phased)
#>      S1 S2 S3
#> 1100  0  0  0
#> 1200  1  1  1
#> 2100  0  1  1
label(phased$EUR, phased)
#>      S1 S2 S3
#> 1100  1  1  0
#> 1200  1  0  0
#> 2100  0  0  0

Take S1 at position 1100. Its genotype is 0|1, so the alternate allele is on haplotype 1. In tract 1, S1’s haplotype 1 is European, so the allele counts as 1 European dosage and 0 African. At position 1200, S1 is 1|1 with one African and one European haplotype, giving 1 African and 1 European. Adding the two matrices gives back each sample’s genotype (its total number of alternate alleles):

label(phased$AFR + phased$EUR, phased)
#>      S1 S2 S3
#> 1100  1  1  0
#> 1200  2  1  1
#> 2100  0  1  1

Splitting without phase (`mode = "dosage"`)

Often phase is unreliable for rare variants, or we prefer not to rely on it. mode = "dosage" uses only the unphased genotype (0, 1 or 2) and each sample’s pair of ancestries in the tract.

unphased <- ancestry_split(vcf_path, msp_path, mode = "dosage", verbose = FALSE)
label(unphased$AFR, unphased)
#>       S1  S2  S3
#> 1100 0.5 0.5 0.0
#> 1200 1.0 0.5 0.5
#> 2100 0.0 0.5 0.5
label(unphased$EUR, unphased)
#>       S1  S2  S3
#> 1100 0.5 0.5 0.0
#> 1200 1.0 0.5 0.5
#> 2100 0.0 0.5 0.5

Every sample in this example has one African and one European haplotype in both tracts, so every heterozygote is ambiguous: the genotype does not say which haplotype carries the allele. Homozygotes are not ambiguous. S1 at 1200 has two alternate alleles, one on each haplotype, so it again gets 1 African and 1 European. An ambiguous heterozygote’s allele is shared between its two populations. The shares are estimated from the carriers of the same variant whose ancestry is unambiguous, and pulled toward 1/2 the less such evidence there is. Here no carrier of a single allele is unambiguous, so every ambiguous allele is split exactly in half. vignette("split-intuition") shows the calculation for less extreme cases. The total is still each sample’s genotype:

label(unphased$AFR + unphased$EUR, unphased)
#>      S1 S2 S3
#> 1100  1  1  0
#> 1200  2  1  1
#> 2100  0  1  1

Part B: More than two populations (K = 4)

Some admixed groups, such as many Hispanic/Latino populations, have ancestry from three or more populations. The population names and codes come from the MSP header, so the same call handles any number of populations. Here we simulate 24 samples, ten rare variants and a single tract in which each haplotype is randomly AFR, EUR, NAT or EAS.

First, the VCF, with random phased genotypes:

set.seed(99)
samples4  <- paste0("S", 1:24)
vcf4_path <- file.path(tmp_dir, "phased4.vcf")

gt_lines <- vapply(1:10, function(i) {
  gts <- sample(c("0|0", "0|1", "1|0", "1|1"), length(samples4), replace = TRUE,
                prob = c(0.85, 0.06, 0.06, 0.03))
  paste(c("19", 1000 + i * 100, ".", "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", samples4),
        collapse = "\t"),
  gt_lines
), vcf4_path)

Then the MSP file. The header now declares four populations, and the single tract (positions 1000 to 3000) holds one random population code per haplotype:

msp4_path <- file.path(tmp_dir, "ancestry4.msp.tsv")
hap_cols  <- paste0(rep(samples4, each = 2), c(".0", ".1"))
hap_anc   <- sample(0:3, 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)

The call is the same as in Part A:

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)

There is now one dosage matrix per population, and ancestry_counts, the number of samples with both haplotypes from each population at each variant, has four columns:

names(phased4)[1:4]
#> [1] "AFR" "EUR" "NAT" "EAS"
head(phased4$ancestry_counts, 3)
#>      AFR EUR NAT EAS
#> [1,]   0   3   1   1
#> [2,]   0   3   1   1
#> [3,]   0   3   1   1

As before, the phased split assigns whole alleles, and the unphased split shares an ambiguous allele between the two populations of that sample’s haplotypes. Here is the first sample and variant where the unphased split gives a fraction, with its dosage in each of the four populations under both modes:

i <- which(unphased4$AFR %% 1 > 0, arr.ind = TRUE)[1, ]
c(variant = phased4$variant_info$pos[i[1]], sample = phased4$sample_ids[i[2]])
#> variant  sample 
#>  "1300"   "S11"
pops <- c("AFR", "EUR", "NAT", "EAS")
rbind(haplotype = sapply(phased4[pops], function(m) m[i[1], i[2]]),
      dosage    = sapply(unphased4[pops], function(m) m[i[1], i[2]]))
#>           AFR EUR NAT EAS
#> haplotype 0.0   0   0 1.0
#> dosage    0.5   0   0 0.5

In haplotype mode the allele goes entirely to one population. In dosage mode it is shared between the two populations of the sample’s haplotypes, and the other two populations get 0.

References

Maples BK, Gravel S, Kenny EE, Bustamante CD (2013). RFMix: a discriminative modeling approach for rapid and robust local-ancestry inference. American Journal of Human Genetics 93(2):278–288. doi: 10.1016/j.ajhg.2013.06.020.

Session information

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

Toy example: splitting a VCF + MSP file (two and more populations)