---
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()` takes two files: a phased VCF with the genotypes and an
[RFMix](https://github.com/slowkoni/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.

```{r toy-vcf}
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.

```{r toy-msp}
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.

```{r toy-split-phased}
phased <- ancestry_split(vcf_path, msp_path, mode = "haplotype", verbose = FALSE)
names(phased)
```

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:

```{r toy-label}
label <- function(m, split) {
  dimnames(m) <- list(split$variant_info$pos, split$sample_ids)
  m
}
label(phased$AFR, phased)
label(phased$EUR, phased)
```

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):

```{r toy-total}
label(phased$AFR + phased$EUR, phased)
```

## 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.

```{r toy-split-dosage}
unphased <- ancestry_split(vcf_path, msp_path, mode = "dosage", verbose = FALSE)
label(unphased$AFR, unphased)
label(unphased$EUR, unphased)
```

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:

```{r toy-total-dosage}
label(unphased$AFR + unphased$EUR, unphased)
```

## Related functions

`ancestry_split_dosage()` and `ancestry_split_phased()` are older
interfaces to the same split, kept so existing code keeps working. New
code should use `ancestry_split()`. `ancestry_split()` reads both files
every time it is called, so for a large cohort it is worth saving the
result (for example with `saveRDS()`) if you need it more than once.

***

# 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:

```{r k4-vcf}
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:

```{r k4-msp}
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:

```{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)
```

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:

```{r k4-show}
names(phased4)[1:4]
head(phased4$ancestry_counts, 3)
```

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:

```{r k4-compare}
i <- which(unphased4$AFR %% 1 > 0, arr.ind = TRUE)[1, ]
c(variant = phased4$variant_info$pos[i[1]], sample = phased4$sample_ids[i[2]])
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]]))
```

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.

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

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

## Session information

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