## ----setup, include=FALSE----------------------------------------------------- knitr::opts_chunk$set(collapse = TRUE, comment = "#>") library(lantern) ## ----a1-data------------------------------------------------------------------ gt <- matrix(c(1, 2, 1, 1, 0, 2), nrow = 1, dimnames = list("var1", paste0("S", 1:6))) anc <- matrix(c(3, 3, 1, 2, 2, 2), nrow = 1, dimnames = list("var1", paste0("S", 1:6))) rbind(genotype = gt[1, ], ancestry = anc[1, ]) ## ----a1-by-hand--------------------------------------------------------------- N <- function(a, g) sum(anc == a & gt == g) # codes: 3 = AFR/AFR, 2 = AFR/EUR, 1 = EUR/EUR n_afr <- 2 * N(3, 2) + N(3, 1) + N(2, 2) n_eur <- 2 * N(1, 2) + N(1, 1) + N(2, 2) D <- n_afr + n_eur c(n_AFR = n_afr, n_EUR = n_eur, D = D, P_AFR = n_afr / D, P_EUR = n_eur / D) ## ----a1-split----------------------------------------------------------------- res <- split_diploid(gt, anc) split_table <- rbind(AFR = res$african[1, ], EUR = res$european[1, ]) colnames(split_table) <- colnames(gt) split_table ## ----a1-edge------------------------------------------------------------------ gt_edge <- matrix(c(1, 0, 1), nrow = 1) anc_edge <- matrix(c(2, 1, 2), nrow = 1) split_diploid(gt_edge, anc_edge) ## ----a2-data------------------------------------------------------------------ ids <- list(paste0("var", 1:3), paste0("S", 1:3)) gt0 <- matrix(c(1, 0, 1, 0, 1, 0, 1, 1, 0), nrow = 3, dimnames = ids) gt1 <- matrix(c(0, 1, 0, 1, 0, 1, 0, 1, 1), nrow = 3, dimnames = ids) a0 <- matrix(c(0, 0, 1, 1, 0, 0, 0, 1, 0), nrow = 3, dimnames = ids) a1 <- matrix(c(1, 1, 0, 0, 1, 1, 1, 0, 1), nrow = 3, dimnames = ids) ## ----a2-split----------------------------------------------------------------- phased <- split_haplotype(gt0, gt1, a0, a1) phased$african phased$european ## ----a2-total----------------------------------------------------------------- all.equal(phased$african + phased$european, unname(gt0 + gt1)) ## ----b1-codes----------------------------------------------------------------- pure_codes <- c(AFR = 3L, EUR = 1L, NAT = 4L) mixed_codes <- data.frame(code = c(2L, 5L, 6L), pop1 = c("AFR", "AFR", "EUR"), pop2 = c("EUR", "NAT", "NAT")) mixed_codes ## ----b1-data------------------------------------------------------------------ set.seed(2025) gt3 <- matrix(sample(0:2, 12 * 30, replace = TRUE, prob = c(0.88, 0.10, 0.02)), nrow = 12, ncol = 30) anc3 <- matrix(sample(c(pure_codes, mixed_codes$code), 12 * 30, replace = TRUE, prob = c(0.30, 0.20, 0.10, 0.20, 0.10, 0.10)), nrow = 12, ncol = 30) out <- split_diploid_multi(gt3, anc3, pure_codes, mixed_codes) names(out) ## ----b1-by-hand--------------------------------------------------------------- v <- which(rowSums(anc3 == 2 & gt3 == 1) > 0)[1] g <- gt3[v, ] a <- anc3[v, ] N <- function(code, genotype) sum(a == code & g == genotype) # codes: 3 = AFR/AFR, 1 = EUR/EUR, 4 = NAT/NAT, 2 = AFR/EUR, 5 = AFR/NAT, 6 = EUR/NAT n_k <- c(AFR = 2 * N(3, 2) + N(3, 1) + N(2, 2) + N(5, 2), EUR = 2 * N(1, 2) + N(1, 1) + N(2, 2) + N(6, 2), NAT = 2 * N(4, 2) + N(4, 1) + N(5, 2) + N(6, 2)) n_k ## ----b1-check----------------------------------------------------------------- s <- which(a == 2 & g == 1)[1] c(by_hand = unname(n_k["AFR"] / (n_k["AFR"] + n_k["EUR"])), lantern = out$AFR[v, s]) ## ----b1-equivalent------------------------------------------------------------ d2 <- split_diploid_multi(gt, anc) all.equal(unname(d2$AFR), unname(res$african)) ## ----b2-split----------------------------------------------------------------- set.seed(2025) hap0 <- matrix(rbinom(8 * 6, 1, 0.1), nrow = 8) hap1 <- matrix(rbinom(8 * 6, 1, 0.1), nrow = 8) anc0 <- matrix(sample(0:2, 8 * 6, replace = TRUE, prob = c(0.4, 0.4, 0.2)), nrow = 8) anc1 <- matrix(sample(0:2, 8 * 6, replace = TRUE, prob = c(0.4, 0.4, 0.2)), nrow = 8) phased3 <- split_haplotype_multi(hap0, hap1, anc0, anc1, pop_codes = c(AFR = 0L, EUR = 1L, NAT = 2L)) phased3$AFR[1:4, ] ## ----b2-total----------------------------------------------------------------- all.equal(phased3$AFR + phased3$EUR + phased3$NAT, hap0 + hap1) ## ----b2-equivalent------------------------------------------------------------ m2 <- split_haplotype_multi(gt0, gt1, a0, a1, pop_codes = c(AFR = 0L, EUR = 1L)) all.equal(unname(m2$AFR), phased$african) ## ----b3-four------------------------------------------------------------------ pure_4 <- c(AFR = 3L, EUR = 1L, NAT = 4L, EAS = 5L) mixed_4 <- data.frame(code = c(2L, 6L, 7L, 8L, 9L, 10L), pop1 = c("AFR", "AFR", "AFR", "EUR", "EUR", "NAT"), pop2 = c("EUR", "NAT", "EAS", "NAT", "EAS", "EAS")) set.seed(11) gt4 <- matrix(sample(0:2, 6 * 16, replace = TRUE, prob = c(0.85, 0.12, 0.03)), nrow = 6) anc4 <- matrix(sample(c(pure_4, mixed_4$code), 6 * 16, replace = TRUE), nrow = 6) out4 <- split_diploid_multi(gt4, anc4, pure_4, mixed_4) names(out4) ## ----c-data------------------------------------------------------------------- gt_dom <- matrix(c(rep(1L, 10), 2L), nrow = 1) anc_dom <- matrix(c(rep(2L, 10), 1L), nrow = 1) raw <- split_diploid_multi(gt_dom, anc_dom, pure_codes, mixed_codes) c(AFR = sum(raw$AFR), EUR = sum(raw$EUR), NAT = sum(raw$NAT)) ## ----c-shrink----------------------------------------------------------------- shrunk <- split_diploid_multi(gt_dom, anc_dom, pure_codes, mixed_codes, shrink = TRUE) rbind(raw = c(AFR = sum(raw$AFR), EUR = sum(raw$EUR)), shrunk = c(AFR = sum(shrunk$AFR), EUR = sum(shrunk$EUR))) ## ----c-share------------------------------------------------------------------ shrunk$AFR[1, 1] ## ----c-other-pop-------------------------------------------------------------- gt_nat <- cbind(gt_dom, matrix(1L, nrow = 1, ncol = 100)) anc_nat <- cbind(anc_dom, matrix(pure_codes[["NAT"]], nrow = 1, ncol = 100)) split_diploid_multi(gt_nat, anc_nat, pure_codes, mixed_codes, shrink = TRUE)$AFR[1, 1] ## ----c-a1--------------------------------------------------------------------- # S4 is the 4th sample split_diploid(gt, anc, shrink = TRUE)$african[1, 4] ## ----session------------------------------------------------------------------ sessionInfo()