vignettes/lantern.Rmd
lantern.RmdIn admixed individuals, each stretch of the genome is inherited from
one of several ancestral populations. A rare-variant burden test that
pools all carriers can miss an effect that exists in only one ancestral
background, or wrongly attribute one. lantern uses
local ancestry to split each sample’s alternate-allele dosage into
per-ancestry dosages (e.g. African and European). It then runs
gene-level SMMAT tests (GMMAT) separately
for each ancestry and combines the per-ancestry p-values with the Cauchy
combination test.
The workflow has three steps:
.msp.tsv file with ancestry_split().write_ancestry_gds().ancestry_smmat().
if (!requireNamespace("BiocManager", quietly = TRUE))
install.packages("BiocManager")
BiocManager::install("lantern")lantern reads VCFs with bcftools when it is on
PATH, which is fastest and also reads BCF. Otherwise it
falls back to SeqArray,
so no external tools are needed.
The package ships a small simulated dataset: 60 admixed samples, 200
rare variants on chromosome 19, two ancestral populations (AFR, EUR) and
8 genes. The phenotype carries an effect that is driven only by
African-ancestry alleles in GENE2.
ext <- system.file("extdata", package = "lantern")
vcf <- file.path(ext, "toy.vcf.gz") # phased genotypes
msp <- file.path(ext, "toy.msp.tsv") # RFMix local ancestry
substr(readLines(msp, n = 3), 1, 70)
#> [1] "#Subpopulation order/codes: AFR=0\tEUR=1"
#> [2] "#chm\tspos\tepos\tsgpos\tegpos\tn snps\tS001.0\tS001.1\tS002.0\tS002.1\tS003.0\tS"
#> [3] "chr19\t1e+06\t4999999\t1.00\t5.00\t25\t0\t0\t0\t0\t0\t0\t0\t0\t0\t1\t0\t1\t0\t0\t1\t0\t0\t0\t0"mode = "haplotype" assigns each phased allele to the
ancestry of its own haplotype, so it is exact.
mode = "dosage" ignores phase: it splits ambiguous
mixed-ancestry heterozygotes using population proportions estimated from
the unambiguous carriers, shrunk toward each chromosome arm’s global
ancestry. See vignette("split-intuition") for the
details.
split <- ancestry_split(vcf, msp, mode = "dosage", chrom = "chr19",
verbose = FALSE)
names(split)
#> [1] "AFR" "EUR" "variant_info" "sample_ids"
#> [5] "mode" "ancestry_counts" "tract_info" "overlap"
dim(split$AFR)
#> [1] 184 60
split$AFR[1:4, 1:6]
#> [,1] [,2] [,3] [,4] [,5] [,6]
#> [1,] 0 0 0 0 0 0
#> [2,] 0 0 0 1 0 0
#> [3,] 0 0 0 0 0 0
#> [4,] 0 0 0 0 0 0
split$EUR[1:4, 1:6]
#> [,1] [,2] [,3] [,4] [,5] [,6]
#> [1,] 0 0 0 0 0 0
#> [2,] 0 0 0 0 0 0
#> [3,] 0 0 0 0 0 0
#> [4,] 0 0 0 0 0 0The result also records which variants and samples were kept:
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
split$overlap[c("n_common", "n_variants_kept", "n_monomorphic_filtered")]
#> $n_common
#> [1] 60
#>
#> $n_variants_kept
#> [1] 184
#>
#> $n_monomorphic_filtered
#> [1] 16
out_dir <- tempfile("lantern-")
gds <- write_ancestry_gds(split, out_dir, verbose = FALSE)
gds
#> $AFR
#> [1] "/tmp/Rtmp7eWPFD/lantern-422647e4df14/AFR.gds"
#>
#> $EUR
#> [1] "/tmp/Rtmp7eWPFD/lantern-422647e4df14/EUR.gds"ancestry_smmat() fits the null model once with
GMMAT::glmmkin(), then runs GMMAT::SMMAT() on
each ancestry’s GDS file. Gene-level weights for the Cauchy combination
are the median per-variant carrier counts in each ancestry. The samples
here are unrelated, so the kinship matrix is the identity.
pheno <- read.delim(file.path(ext, "toy_pheno.tsv"))
kinship <- diag(nrow(pheno))
dimnames(kinship) <- list(pheno$id, pheno$id)
res <- ancestry_smmat(gds, pheno, y ~ age + sex, kinship,
gene_group_file = file.path(ext, "toy_genes.tsv"),
ancestry_counts = split$ancestry_counts,
variant_info = split$variant_info,
verbose = FALSE)
res$results
#> gene p_AFR p_EUR w_AFR w_EUR p_cauchy
#> 1 GENE1 8.429637e-01 0.2321365 42 1 8.378715e-01
#> 2 GENE2 1.860597e-08 0.8339380 39 2 1.956013e-08
#> 3 GENE3 5.076218e-01 0.2765924 41 3 4.887590e-01
#> 4 GENE4 3.937408e-01 0.7851541 40 1 4.049334e-01
#> 5 GENE5 2.826767e-01 0.2007261 40 2 2.776638e-01
#> 6 GENE6 6.638040e-01 0.5783176 33 2 6.594389e-01
#> 7 GENE7 1.517208e-01 0.4259828 41 1 1.544780e-01
#> 8 GENE8 6.060560e-01 0.3600061 33 3 5.863210e-01GENE2 is strongly associated in the AFR test but not in
the EUR test, which matches the simulated African-specific effect. A
test that pools both ancestries would dilute that signal.
vignette("split-intuition"): what the split computes,
on small matrices, for two or more populations.vignette("toy-vcf-msp-example"): the VCF and MSP input
formats, with four populations.vignette("real-data-chr19"): a real-cohort analysis
with binary, continuous and count phenotypes.
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 BiocStyle_2.40.0
#>
#> loaded via a namespace (and not attached):
#> [1] tidyselect_1.2.1 dplyr_1.2.1 Biostrings_2.80.2
#> [4] gdsfmt_1.48.2 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-6 magrittr_2.0.5 compiler_4.6.1
#> [13] rlang_1.3.0 sass_0.4.10 tools_4.6.1
#> [16] yaml_2.3.12 data.table_1.18.6.1 knitr_1.52
#> [19] SeqArray_1.52.1 purrr_1.2.2 BiocGenerics_0.58.1
#> [22] GWASExactHW_1.2 desc_1.4.3 nnet_7.3-20
#> [25] grid_4.6.1 stats4_4.6.1 jomo_2.7-6
#> [28] GMMAT_1.5.0 mice_3.19.0 iterators_1.0.14
#> [31] MASS_7.3-65 cli_3.6.6 rmarkdown_2.32
#> [34] crayon_1.5.3 ragg_1.5.2 reformulas_0.4.4
#> [37] generics_0.1.4 otel_0.2.0 minqa_1.2.8
#> [40] cachem_1.1.0 operator.tools_1.6.3.1 splines_4.6.1
#> [43] parallel_4.6.1 BiocManager_1.30.27 XVector_0.52.0
#> [46] vctrs_0.7.3 boot_1.3-32 glmnet_5.1
#> [49] Matrix_1.7-5 jsonlite_2.0.0 bookdown_0.48
#> [52] IRanges_2.46.0 S4Vectors_0.50.3 mitml_0.4-5
#> [55] systemfonts_1.3.2 foreach_1.5.2 tidyr_1.3.2
#> [58] jquerylib_0.1.4 glue_1.8.1 pkgdown_2.2.1
#> [61] nloptr_2.2.1 pan_2.0 codetools_0.2-20
#> [64] shape_1.4.6.1 GenomicRanges_1.64.0 lme4_2.0-6
#> [67] tibble_3.3.1 pillar_1.11.1 htmltools_0.5.9
#> [70] Seqinfo_1.2.0 R6_2.6.1 textshaping_1.0.5
#> [73] SeqVarTools_1.50.1 Rdpack_2.6.6 formula.tools_1.7.1
#> [76] evaluate_1.0.5 lattice_0.22-9 Biobase_2.72.0
#> [79] rbibutils_2.4.1 backports_1.5.1 broom_1.0.13
#> [82] bslib_0.12.0 Rcpp_1.1.2 nlme_3.1-169
#> [85] mgcv_1.9-4 logistf_1.26.1 xfun_0.61
#> [88] fs_2.1.0 pkgconfig_2.0.3