Split a phased VCF by local ancestry (Step 1)

Parse an RFMix MSP file and a phased VCF, intersect samples, map each variant to its ancestry tract, and split genotypes into per-population dosage matrices using either the proportional ("dosage") or deterministic ("haplotype") split algorithm. This is the package's Step 1 entry point: it only computes in-memory dosage matrices, with no file I/O. Pass the result to write_ancestry_gds (Step 2) to materialize GDS files, and then to ancestry_smmat (Step 3) to run SMMAT + Cauchy combination.

ancestry_split(
  vcf_path,
  msp_path,
  mode = c("dosage", "haplotype"),
  chrom = NULL,
  use_gla = TRUE,
  verbose = TRUE
)

Arguments

vcf_path

Path to phased VCF/BCF file (plain or gzipped). Read with bcftools if it is on PATH (fastest; required for BCF), otherwise with SeqArray. Set options(lantern.vcf_reader = "bcftools") or "seqarray" to force one.

msp_path

Path to RFMix MSP file (plain text or gzipped TSV).

mode

Split algorithm: "dosage" (proportional p1/p2 split via split_diploid_multi, unphased) or "haplotype" (deterministic per-haplotype split via split_haplotype_multi, phased). Either/or – call ancestry_split() twice if both are needed.

chrom

Chromosome to process (e.g., "chr19" or "19"). If NULL, all chromosomes present in the VCF are used. When supplied, vcf_path and msp_path may contain other chromosomes too: MSP tracts outside chrom are dropped up front, and bcftools query is restricted to chrom via -t when the VCF's declared or inferred contigs allow resolving the exact name to use (whichever of chrom, with, or without a leading "chr" matches) – other chromosomes' genotypes are never read into R. If the contig can't be resolved, all variants are queried and filtered in R instead (slower, more memory, same result).

use_gla

Logical, mode = "dosage" only. When TRUE (default), per-arm GLA (global local ancestry) shrinkage is computed and applied automatically (see the "Core Algorithm" section of the package's CLAUDE.md for the formula/rationale). When FALSE, GLA shrinkage is skipped entirely and the split reproduces the original pre-shrinkage p[k] estimator exactly, including its 0.5/0.5 singleton fallback – equivalent to calling split_diploid_multi with gla = NULL. Ignored (no effect) when mode = "haplotype", which is deterministic and has no ambiguous ratio to shrink.

verbose

Print step-by-step progress messages.

Value

Invisibly, a list with one named numeric matrix (variants x samples) per population (named after the MSP population codes, e.g. $AFR, $EUR), plus:

variant_info

data.frame with columns chrom, pos, ref, alt.

sample_ids

Character vector of common sample IDs.

mode

The split mode used.

ancestry_counts

Numeric matrix (variants x populations): for each variant, the number of samples with pure ancestry in that population. Feeds ancestry_smmat's per-gene weight calculation.

tract_info

data.frame of ancestry tracts from the MSP file.

overlap

List of intersection/filtering statistics.

See also

write_ancestry_gds, ancestry_smmat

Examples

vcf <- system.file("extdata", "toy.vcf.gz", package = "lantern")
msp <- system.file("extdata", "toy.msp.tsv", package = "lantern")
split <- ancestry_split(vcf, msp, mode = "dosage", chrom = "chr19",
                        verbose = FALSE)
split$AFR[1:3, 1:5]   # African dosages (variants x samples)
#>      [,1] [,2] [,3] [,4] [,5]
#> [1,]    0    0    0    0    0
#> [2,]    0    0    0    1    0
#> [3,]    0    0    0    0    0
split$EUR[1:3, 1:5]   # European dosages
#>      [,1] [,2] [,3] [,4] [,5]
#> [1,]    0    0    0    0    0
#> [2,]    0    0    0    0    0
#> [3,]    0    0    0    0    0
head(split$variant_info)
#>   chrom     pos ref alt
#> 1 chr19 1168394   T   A
#> 2 chr19 1240171   T   C
#> 3 chr19 2252809   A   G
#> 4 chr19 2330989   C   T
#> 5 chr19 2436455   A   G
#> 6 chr19 2459485   C   G