---
title: "Getting started with lantern"
author:
  - name: Yu Wang
  - name: Bjoernar Tuftin
  - name: Paul Auer
output: BiocStyle::html_document
vignette: >
  %\VignetteIndexEntry{Getting started with lantern}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r setup, include = FALSE}
knitr::opts_chunk$set(collapse = TRUE, comment = "#>")
```

# Introduction

People with ancestry from more than one population, such as African
Americans, European, are called admixed. In each of their chromosomes, some
segments were inherited from African ancestors and others from European
ancestors. Software such as [RFMix](https://github.com/slowkoni/rfmix)
(Maples et al. 2013) can infer, for every segment in each of the two copies of a chromosome,
which population it came from. This is so called local ancestry.

Admixed people are often left out of genetic association studies, or
analysed as if their whole genome came from one population. Yet the same
allele can relate to a trait differently depending on the ancestral
background it sits on. Linkage disequilibrium differs between
populations, so a variant may tag a nearby causal variant on African
haplotypes but not on European ones. Allele frequencies differ, so a
variant can be common on one background and almost absent on the other.
And an allele's effect can depend on the other variants carried on the
same haplotype. Tractor (Atkinson et al. 2021) uses local ancestry to
estimate a separate effect for each ancestry in genome-wide association
studies of common variants, and showed that this gives accurate
ancestry-specific effect sizes, boosts power, and places association
signals closer to the causal variants.

Rare-variant association tests have the same problem. The rare
variants in a gene are grouped and tested whether carrying them is associated with a
trait. If the variants act on African haplotypes but not on European ones
(or the reverse), pooling all carriers mixes carriers with and without the
effect, the signal is diluted, and the test loses power.

To bridge the gap, lantern aims to use local ancestry to split each
person's count of alternate alleles (their dosage, 0, 1 or 2) into one dosage
per ancestry. A heterozygote whose allele sits on an African haplotype contributes
1 to the African dosage and 0 to the European dosage. lantern then runs the
gene-based SMMAT test from the `r CRANpkg("GMMAT")` package separately for each ancestry,
and combines the per-ancestry p-values for each gene into one p-value with
the Cauchy combination test.

## Which haplotype carries the rare allele?

For most carriers the split is clear. If both copies of the chromosome are
African at the variant (person A below), the allele is African; if both
are European (person B), it is European. The hard case is a heterozygote
with one African and one European copy (person C): the allele is on one of
the two, and only phase, the assignment of each allele to a chromosome
copy, says which.

![Local ancestry at a rare variant in three people. Only person C, with one African and one European haplotype, needs phase to tell which ancestry the allele belongs to.](figures/haplotype-ambiguity.svg){width=100%}

Phasing is least reliable for exactly the variants that rare-variant tests
use. In the UK Biobank, SHAPEIT5 places a singleton allele on the wrong
haplotype about one time in three:

![SHAPEIT5 switch-error rate by minor allele count in the UK Biobank (Hofmeister et al. 2023).](figures/switch-error.svg){width=100%}

lantern therefore offers two ways to split. `mode = "haplotype"` uses the
phased genotypes, like Tractor does for common variants. `mode = "dosage"`
does not use phase: it shares each ambiguous heterozygote's allele between
the two ancestries, with shares estimated from the variant's other carriers
(`vignette("split-intuition")` explains how).

## What the split gains

We simulated traits on chromosome 19 of 3,313 participants of the Jackson
Heart Study, about 80% of whose ancestry is African, with causal effects
carried by alleles on African haplotypes, on European haplotypes, or on
both. The ancestry-split tests did not inflate the type I error, down to
a significance level of \(2.5 \times 10^{-6}\). When the effect was carried on the
European (here the minority-ancestry) haplotypes, testing the European
dosage had about twice the power of SMMAT on unsplit genotypes, and the
unphased split about 1.5 times. When the effect was on African haplotypes
or on both, the African-ancestry test lost almost nothing against SMMAT.

![Power at a significance level of 0.05 by mean effect size, averaged over 100 genes (Wang et al. 2026). AF, EU: tests of the African and European dosages; Phased, Unphased: `mode = "haplotype"` and `mode = "dosage"`.](figures/power.svg){width=100%}

The details are in the preprint: Wang Y, Tuftin B, Raffield LM, Hidalgo B,
Kerns SL, DeWan AT, Leal SM, Auer PL (2026). LANTERN: Leveraging local
ancestry tracts to enhance rare-variant aggregate association testing.
*medRxiv*. doi:
[10.64898/2026.04.24.26351693](https://doi.org/10.64898/2026.04.24.26351693).

## Workflow

The analysis has three steps as follows:

1. `ancestry_split()` reads a phased VCF and an RFMix local-ancestry file
   and returns one dosage matrix per ancestry.
2. `write_ancestry_gds()` saves each matrix as a `r Biocpkg("SeqArray")`
   GDS file, the format GMMAT reads.
3. `ancestry_smmat()` runs SMMAT on each GDS file and combines the
   results.

# Installation

```{r install, eval = FALSE}
if (!requireNamespace("BiocManager", quietly = TRUE))
    install.packages("BiocManager")
BiocManager::install("lantern")
```

```{r load}
library(lantern)
```

# Example data

lantern includes a small simulated dataset: 60 admixed individuals, 200
rare variants on chromosome 19 and eight genes. We simulated a continuous
trait `y` that is raised by alternate alleles in `GENE2`, but only when the
allele sits on an African haplotype. A good analysis should find `GENE2`
in the African-ancestry test and not in the European one.

```{r files}
ext <- system.file("extdata", package = "lantern")
vcf <- file.path(ext, "toy.vcf.gz")
msp <- file.path(ext, "toy.msp.tsv")
```

The VCF holds phased genotypes. A `|` separates the two copies of the
chromosome, so `0|1` means the alternate allele is on the second copy:

```{r vcf}
vcf_lines <- readLines(vcf)
data_lines <- vcf_lines[!startsWith(vcf_lines, "##")]
substr(data_lines[1:3], 1, 60)
```

The MSP file is RFMix's local-ancestry output. Its first line gives the
population codes (here 0 = AFR and 1 = EUR). Each later row is a segment
of the chromosome, with one ancestry call per sample and chromosome copy
(columns `S001.0`, `S001.1`, ...):

```{r msp}
substr(readLines(msp, n = 4), 1, 70)
```

# Step 1: split dosages by local ancestry

`ancestry_split()` matches samples between the two files, finds the
ancestry segment that contains each variant, and splits the dosages.

```{r split}
split <- ancestry_split(vcf, msp, mode = "dosage", chrom = "chr19",
                        verbose = FALSE)
```

The result has one variants x samples dosage matrix per population, named
after the populations in the MSP file:

```{r split-show}
split$AFR[1:4, 1:6]
split$EUR[1:4, 1:6]
```

Most entries are 0 because the variants are rare. `mode = "dosage"`
ignores phase. For a heterozygote with one African and one European
haplotype, the genotype alone does not say which haplotype carries the
allele, so the allele is shared between the two ancestries, with shares
estimated from the other carriers of the same variant. Here is the first
such heterozygote in the data; its two shares add up to its one allele:

```{r split-fraction}
i <- which(split$AFR %% 1 > 0, arr.ind = TRUE)[1, ]
c(AFR = split$AFR[i[1], i[2]], EUR = split$EUR[i[1], i[2]])
```

`vignette("split-intuition")` explains how the shares are calculated. With
phased data you can instead use `mode = "haplotype"`, which assigns every
allele to the ancestry of its own haplotype and never produces fractions.

The result also lists the variants that were kept. Variants with no
alternate allele in these samples are dropped
(`r split$overlap$n_monomorphic_filtered` here):

```{r kept}
head(split$variant_info)
split$overlap[c("n_common", "n_variants_kept", "n_monomorphic_filtered")]
```

# Step 2: write per-ancestry GDS files

GMMAT reads genotypes from GDS files, so each dosage matrix is saved to
its own file. The function returns the file paths, named by population.

```{r gds}
out_dir <- tempfile("lantern-")
gds <- write_ancestry_gds(split, out_dir, verbose = FALSE)
gds
```

# Step 3: test each gene in each ancestry

SMMAT needs three more inputs: a phenotype table, a kinship matrix and a
gene-group file. The phenotype table has one row per sample, with the
trait and covariates:

```{r pheno}
pheno <- read.delim(file.path(ext, "toy_pheno.tsv"))
head(pheno, 3)
```

The kinship matrix describes how related the samples are. Our simulated
samples are unrelated, so we use the identity matrix. In a real study this
would be a genetic relatedness matrix:

```{r kinship}
kinship <- diag(nrow(pheno))
dimnames(kinship) <- list(pheno$id, pheno$id)
```

The gene-group file assigns variants to genes. It has no header and six
columns: gene, chromosome, position, reference allele, alternate allele
and a weight. These are the 25 variants of `GENE2`:

```{r genes}
genes <- read.delim(file.path(ext, "toy_genes.tsv"), header = FALSE,
                    col.names = c("gene", "chr", "pos", "ref", "alt", "weight"))
head(genes[genes$gene == "GENE2", ])
table(genes$gene)
```

`ancestry_smmat()` fits the null model (the trait explained by the
covariates and relatedness, without genotypes) once, then runs SMMAT for
each ancestry. Passing `ancestry_counts` and `variant_info` from Step 1
weights the Cauchy combination by how many samples of each ancestry carry
a gene's variants, so the ancestry with more data counts for more.

```{r smmat}
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
```

The table has one row per gene: `p_AFR` and `p_EUR` are the SMMAT p-values
from each ancestry, `w_AFR` and `w_EUR` the weights, and `p_cauchy` the
combined p-value. `GENE2` stands out:

```{r gene2}
res$results[res$results$gene == "GENE2", ]
```

Its African-ancestry p-value is very small, its European-ancestry p-value
is not, and the combined p-value keeps the signal. For comparison, we test
the same genes without splitting by ancestry. Adding the two dosage
matrices gives back each sample's total dosage:

```{r pooled}
pooled_gds <- file.path(out_dir, "pooled.gds")
write_dosage_gds(split$AFR + split$EUR, split$variant_info,
                 split$sample_ids, pooled_gds)
pooled <- ancestry_smmat(c(POOLED = pooled_gds), pheno, y ~ age + sex,
                         kinship, file.path(ext, "toy_genes.tsv"),
                         verbose = FALSE)
pooled$results[pooled$results$gene == "GENE2", c("gene", "p_POOLED")]
```

The pooled test also detects `GENE2`, with a p-value about three times
larger: carriers whose allele sits on a European haplotype add to the
pooled dosage without adding to the effect.

```{r cleanup, include = FALSE}
unlink(out_dir, recursive = TRUE)
```

# Further reading

- `vignette("split-intuition")` works through the split on small
  matrices, for two and for more populations.
- `vignette("toy-vcf-msp-example")` shows the VCF and MSP formats in more
  detail, including a four-population example.
- `vignette("real-data-chr19")` analyses chromosome 19 in a real cohort
  with binary, continuous and count traits.

# References

Atkinson EG, Maihofer AX, Kanai M, Martin AR, Karczewski KJ, Santoro ML,
Ulirsch JC, Kamatani Y, Okada Y, Finucane HK, Koenen KC, Nievergelt CM,
Daly MJ, Neale BM (2021). Tractor uses local ancestry to enable the
inclusion of admixed individuals in GWAS and to boost power. *Nature
Genetics* 53(2):195--204. doi:
[10.1038/s41588-020-00766-y](https://doi.org/10.1038/s41588-020-00766-y).

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](https://doi.org/10.1016/j.ajhg.2013.06.020).

# Session information

```{r session}
sessionInfo()
```