## ----setup, include=FALSE-----------------------------------------------------
knitr::opts_chunk$set(collapse = TRUE, comment = "#>")
library(lantern)

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

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

## ----toy-split-phased---------------------------------------------------------
phased <- ancestry_split(vcf_path, msp_path, mode = "haplotype", verbose = FALSE)
names(phased)

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

## ----toy-total----------------------------------------------------------------
label(phased$AFR + phased$EUR, phased)

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

## ----toy-total-dosage---------------------------------------------------------
label(unphased$AFR + unphased$EUR, unphased)

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

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

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

## ----k4-show------------------------------------------------------------------
names(phased4)[1:4]
head(phased4$ancestry_counts, 3)

## ----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]]))

## ----cleanup, include = FALSE-------------------------------------------------
unlink(tmp_dir, recursive = TRUE)

## ----session------------------------------------------------------------------
sessionInfo()