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

## ----install, eval = FALSE----------------------------------------------------
# if (!requireNamespace("BiocManager", quietly = TRUE))
#     install.packages("BiocManager")
# BiocManager::install("lantern")

## ----load---------------------------------------------------------------------
library(lantern)

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

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

## ----msp----------------------------------------------------------------------
substr(readLines(msp, n = 4), 1, 70)

## ----split--------------------------------------------------------------------
split <- ancestry_split(vcf, msp, mode = "dosage", chrom = "chr19",
                        verbose = FALSE)

## ----split-show---------------------------------------------------------------
split$AFR[1:4, 1:6]
split$EUR[1:4, 1:6]

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

## ----kept---------------------------------------------------------------------
head(split$variant_info)
split$overlap[c("n_common", "n_variants_kept", "n_monomorphic_filtered")]

## ----gds----------------------------------------------------------------------
out_dir <- tempfile("lantern-")
gds <- write_ancestry_gds(split, out_dir, verbose = FALSE)
gds

## ----pheno--------------------------------------------------------------------
pheno <- read.delim(file.path(ext, "toy_pheno.tsv"))
head(pheno, 3)

## ----kinship------------------------------------------------------------------
kinship <- diag(nrow(pheno))
dimnames(kinship) <- list(pheno$id, pheno$id)

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

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

## ----gene2--------------------------------------------------------------------
res$results[res$results$gene == "GENE2", ]

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

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

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