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

## ----library, message=FALSE---------------------------------------------------
library(WithinHostExperiment)

## ----library-gr, message=FALSE------------------------------------------------
library(GenomicRanges)

## ----object-show--------------------------------------------------------------
data(example_whe)
example_whe

## ----object-assays------------------------------------------------------------
SummarizedExperiment::assayNames(example_whe)
head(SummarizedExperiment::assay(example_whe, "altFreq"), 4)

## ----object-rowranges---------------------------------------------------------
head(SummarizedExperiment::rowRanges(example_whe), 3)
GenomeInfoDb::seqlengths(example_whe)

## ----object-coldata-----------------------------------------------------------
SummarizedExperiment::colData(example_whe)

## ----object-qclog-------------------------------------------------------------
qcLog(example_whe)

## ----object-qclog-replace-----------------------------------------------------
example_qc <- flagISNV(example_whe, ISNVFilter(minDepth = 1000L))
reviewed <- example_qc
extra <- S4Vectors::DataFrame(
    step = "manual_review", parameter = "reviewer",
    value = "two sites re-checked by eye", n_flagged = 0L,
    n_passed = sum(SummarizedExperiment::assay(example_qc, "qcPass")),
    timestamp = format(Sys.time(), "%Y-%m-%dT%H:%M:%S"))
qcLog(reviewed) <- rbind(qcLog(example_qc), extra)
qcLog(reviewed)[, c("step", "parameter", "value", "n_flagged", "n_passed")]

## ----object-qcfilters---------------------------------------------------------
qcFilters(example_whe)
length(qcFilters(example_qc))
qcFilters(example_qc)[[1]]

## ----object-slots-subset------------------------------------------------------
part <- example_qc[1:5, ]
qcLog(part)$n_passed
nrow(part) * ncol(part)
length(qcFilters(combineRows(example_qc[1:5, ], example_qc[6:10, ])))

## ----object-validity----------------------------------------------------------
gr <- GRanges("seg", IRanges::IRanges(1:3, width = 1))
try(WithinHostExperiment(
    assays = list(
        altFreq = matrix(c(0.1, 0.2, 0.3), ncol = 1),
        qcPass = matrix(c(1L, 1L, 0L), ncol = 1)),
    rowRanges = gr,
    colData = S4Vectors::DataFrame(sample_id = "S1")))
bad_filters <- example_whe
qcFilters(bad_filters) <- list("minDepth = 100")
try(validObject(bad_filters))
bad_log <- example_whe
qcLog(bad_log) <- qcLog(example_whe)[, c("step", "parameter")]
try(validObject(bad_log))

## ----object-coerce------------------------------------------------------------
head(as(example_whe, "data.frame"), 4)
vr <- as(example_whe, "VRanges")
vr
GenomeInfoDb::seqlengths(vr)

## ----import-vcf---------------------------------------------------------------
vcf1 <- system.file("extdata", "test_donor.vcf",
    package = "WithinHostExperiment")
vcf2 <- system.file("extdata", "test_recipient.vcf",
    package = "WithinHostExperiment")
ref <- Biostrings::readDNAStringSet(
    system.file("extdata", "test_ref.fa", package = "WithinHostExperiment"))
whe <- readWithinHost(
    c(vcf1, vcf2),
    colData = S4Vectors::DataFrame(
        sample_id = c("donor", "recipient"),
        host_id   = c("P1", "P2"),
        role      = c("donor", "recipient"),
        pair_id   = c("pair_1", "pair_1")
    ),
    caller = "ivar",
    genome = ref
)
whe
GenomeInfoDb::seqlengths(whe)

## ----import-vcf-coerce--------------------------------------------------------
v <- VariantAnnotation::readVcf(vcf1)
as(v, "WithinHostExperiment")

## ----import-tsv---------------------------------------------------------------
tsv <- system.file("extdata", "test_donor_ivar.tsv",
    package = "WithinHostExperiment")
whe_tsv <- readWithinHostTable(tsv, format = "ivar")
whe_tsv
GenomeInfoDb::seqlengths(whe_tsv)

## ----qc-filter----------------------------------------------------------------
filt <- ISNVFilter(minDepth = 200L, minFreq = 0.03)
filt

## ----qc-flag------------------------------------------------------------------
whe_qc <- flagISNV(whe, filt)
qcLog(whe_qc)

## ----qc-summary---------------------------------------------------------------
qcSummary(whe_qc)

## ----qc-passed----------------------------------------------------------------
clean <- passedISNV(whe_qc)
cat("Before:", nrow(whe_qc), "sites; after passedISNV:", nrow(clean), "sites\n")

## ----qc-adaptive--------------------------------------------------------------
adaptiveFilter("illumina")
adaptiveFilter("ont")
adaptiveFilter("pacbio_hifi")

## ----qc-replicate-------------------------------------------------------------
rep1 <- system.file("extdata", "test_donor_rep1.vcf",
    package = "WithinHostExperiment")
rep2 <- system.file("extdata", "test_donor_rep2.vcf",
    package = "WithinHostExperiment")
whe_rep <- readWithinHost(
    c(rep1, rep2),
    colData = S4Vectors::DataFrame(
        sample_id       = c("rep1", "rep2"),
        replicate_group = c("donor_reps", "donor_reps"),
        replicate_id    = c("1", "2")
    ),
    caller = "ivar", genome = ref
)
rep_checked <- flagReplicateDiscordance(whe_rep, freqTolerance = 0.05)
qcLog(rep_checked)

## ----qc-replicate-rule--------------------------------------------------------
rule <- ConcordanceRule(freqTolerance = 0.05)
rep_qc <- flagISNV(whe_rep,
    ISNVFilter(minDepth = 200L, replicateConc = rule))
qcLog(rep_qc)
qcFilters(rep_qc)[[1]]

## ----qc-primer----------------------------------------------------------------
primers <- GRanges("test_segment",
    IRanges::IRanges(start = c(40, 880), end = c(65, 950)))
whe_primer <- flagPrimerSites(whe_qc, primers)
qcLog(whe_primer)

## ----qc-sci-------------------------------------------------------------------
sci <- sciStrandBias(
    fwd_alt   = c(50, 2, 30),
    rev_alt   = c(45, 40, 28),
    fwd_total = c(500, 500, 500),
    rev_total = c(480, 480, 480))
sci

## ----qc-contamination---------------------------------------------------------
gr_c <- GRanges("seg", IRanges::IRanges(c(100, 200, 300), width = 1))
mcols(gr_c)$ref <- c("A", "C", "G")
mcols(gr_c)$alt <- c("T", "G", "A")
freq_c <- matrix(c(0.80, 0.05, 0.30, 0.01, 0.40, 0.00), nrow = 3, ncol = 2)
whe_c <- WithinHostExperiment(
    assays = list(altFreq = freq_c), rowRanges = gr_c,
    colData = S4Vectors::DataFrame(sample_id = c("A", "B"),
        run = c("run1", "run1")))
detectCrossContamination(whe_c, runCol = "run")

## ----qc-pipe------------------------------------------------------------------
result <- whe |>
    flagISNV(ISNVFilter(minDepth = 200L, minFreq = 0.03)) |>
    flagPrimerSites(primers) |>
    passedISNV()
cat("Input:", nrow(whe), "sites; after QC and primer masking:",
    nrow(result), "sites\n")

## ----annotate-genes-----------------------------------------------------------
gff_file <- tempfile(fileext = ".gff3")
writeLines(c(
    "##gff-version 3",
    "test_segment\t.\tgene\t1\t500\t.\t+\t.\tID=gene-A;Name=geneA",
    "test_segment\t.\tgene\t501\t1000\t.\t+\t.\tID=gene-B;Name=geneB"
), gff_file)
features <- rtracklayer::import(gff_file)
genes <- features[features$type == "gene"]
whe_ann <- annotateGenes(whe_qc, genes, nameCol = "Name")
table(mcols(SummarizedExperiment::rowRanges(whe_ann))$GENE, useNA = "ifany")

## ----annotate-codon-----------------------------------------------------------
gff_cds <- tempfile(fileext = ".gff3")
writeLines(c("##gff-version 3",
    "test_segment\t.\tCDS\t1\t999\t.\t+\t0\tID=cds-1;gene=testGene"), gff_cds)
cds <- rtracklayer::import(gff_cds)
cds <- cds[cds$type == "CDS"]
whe_codon <- annotateCodonChange(whe_qc, cds, ref)
mc <- mcols(SummarizedExperiment::rowRanges(whe_codon))
table(mc$CONSEQUENCE, useNA = "ifany")
head(as.data.frame(mc[, c("GENEID", "PROTEINLOC", "REFCODON", "VARCODON",
    "REFAA", "VARAA", "CONSEQUENCE")]), 4)

## ----diversity----------------------------------------------------------------
div <- calcDiversity(whe_qc)
div

## ----diversity-plausibility---------------------------------------------------
checkDiversityPlausibility(div, pathogen = "sars_cov_2")

## ----diversity-model----------------------------------------------------------
div_df <- data.frame(
    shannon = c(1.2, 1.5, 0.8, 1.1, 1.9, 0.7),
    timepoint = c(1, 2, 3, 1, 2, 3),
    group = c("A", "A", "A", "B", "B", "B"))
res <- fitDiversityModel(div_df, shannon ~ timepoint, group_col = "group")
res$summary_df
multitestCorrection(c(0.001, 0.04, 0.03, 0.5))

## ----sfs-build----------------------------------------------------------------
sfs_donor <- buildSFS(whe_qc, sampleIdx = "donor", nBins = 10L)
sfs_donor

## ----sfs-correct--------------------------------------------------------------
correctSFSBias(sfs_donor, method = "truncation")
correctSFSBias(sfs_donor, method = "binomial")

## ----sfs-compare--------------------------------------------------------------
sfs_recip <- buildSFS(whe_qc, sampleIdx = "recipient", nBins = 10L)
cmp <- compareSFS(sfs_donor, sfs_recip)
cmp[c("chisq_stat", "df", "chisq_p", "min_expected", "chisq_approx_ok",
    "n1", "n2")]

## ----neutrality-tests---------------------------------------------------------
nt <- calcNeutralityTests(whe_qc)
nt

## ----neutrality-sfs-----------------------------------------------------------
whe_h <- flagISNV(whe,
    ISNVFilter(minDepth = 200L, minFreq = 0.03, maxFreq = 1))
sfs_unf <- buildSFS(whe_h, sampleIdx = "donor", nBins = 10L, fold = FALSE)
neutralityFromSFS(sfs_unf)

## ----selection-pnps-----------------------------------------------------------
pnps <- pNpSWithinHost(whe_qc, cds, ref)
pnps[, c("sample_id", "nS", "nN", "S_sites", "N_sites", "pNpS")]

## ----selection-consensus------------------------------------------------------
cons <- consensusSequences(whe_qc, ref)
S4Vectors::mcols(cons)

## ----selection-kaks, message=FALSE--------------------------------------------
cds_region <- GRanges("test_segment", IRanges::IRanges(1, 999), strand = "+")
cons_cds <- consensusSequences(whe_qc, ref, minFreq = 0.03, region = cds_region)
if (requireNamespace("MSA2dist", quietly = TRUE)) {
    kaks <- cons_cds |> MSA2dist::dnastring2kaks(model = "Li")
    kaks[, c("seq1", "seq2", "Ka", "Ks", "Ka/Ks")]
}

## ----selection-coefficient----------------------------------------------------
freqs <- c(0.05, 0.10, 0.25, 0.50)
s_hat <- estimateSelectionCoefficient(freqs, times = c(0, 3, 7, 14))
s_hat

## ----transmission-pairs-------------------------------------------------------
transmissionPairs(whe_qc)

## ----transmission-shared------------------------------------------------------
sv <- calcSharedVariants(whe_qc, pairId = "pair_1")
cat("Shared:", sv$n_shared, "| donor-only:", sv$n_donor_only,
    "| recipient-only:", sv$n_recipient_only, "\n")
cat("Donor-only at adequate / low / unknown recipient depth:",
    sv$n_donor_only_adequate_depth, "/", sv$n_donor_only_low_depth, "/",
    sv$n_donor_only_unknown_depth, "| recipient call failed QC:",
    sv$n_donor_only_recipient_failed_qc, "\n")

## ----transmission-quick-------------------------------------------------------
nb <- quickBottleneck(whe_qc, pairId = "pair_1", maxNb = 50L)
cat("Nb:", nb$Nb, "with 95% CI [", nb$ci[1], ",", nb$ci[2], "],",
    nb$n_variants, "donor variants, log-likelihood", round(nb$loglik, 2), "\n")

## ----transmission-exact-------------------------------------------------------
nb_pa <- exactBottleneck(whe_qc, pairId = "pair_1", maxNb = 50L,
    method = "presence_absence")
cat("Presence-absence Nb:", nb_pa$Nb, "with 95% CI [", nb_pa$ci[1], ",",
    nb_pa$ci[2], "] from", nb_pa$n_variants, "variants\n")

## ----transmission-export------------------------------------------------------
head(asViralBottleneckInput(whe_qc, pairId = "pair_1"), 4)
head(exportPairFrequencies(whe_qc, pairId = "pair_1"), 4)

## ----longitudinal-track-------------------------------------------------------
gr_l <- GRanges("seg1", IRanges::IRanges(c(100, 200), width = 1))
mcols(gr_l)$ref <- c("A", "C")
mcols(gr_l)$alt <- c("T", "G")
freq_l <- matrix(c(0.02, 0.10, 0.12, 0.15), nrow = 2, ncol = 2)
whe_l <- WithinHostExperiment(
    assays = list(altFreq = freq_l), rowRanges = gr_l,
    colData = S4Vectors::DataFrame(
        sample_id = c("t1", "t2"), host_id = c("H1", "H1"),
        timepoint = c(1, 2)))
trackFrequency(whe_l)

## ----longitudinal-emerging----------------------------------------------------
detectEmergingVariants(whe_l)

## ----longitudinal-temporal----------------------------------------------------
gr_t <- GRanges("seg1", IRanges::IRanges(c(100, 200, 300), width = 1))
mcols(gr_t)$ref <- c("A", "C", "G")
mcols(gr_t)$alt <- c("T", "G", "A")
freq_t <- matrix(c(0.1, 0.05, 0.15, NA, NA, 0.2, 0.2, NA, 0.1),
    nrow = 3, ncol = 3)
whe_t <- WithinHostExperiment(
    assays = list(altFreq = freq_t), rowRanges = gr_t,
    colData = S4Vectors::DataFrame(
        sample_id = c("t1", "t2", "t3"), host_id = c("H1", "H1", "H1"),
        timepoint = c(1, 2, 3)))
whe_t <- flagTemporalInconsistency(whe_t, minTimepoints = 2L)
mcols(SummarizedExperiment::rowRanges(whe_t))$temporal_class
qcLog(whe_t)

## ----callers-consensus--------------------------------------------------------
gr1 <- GRanges("seg", IRanges::IRanges(c(100, 200, 300), width = 1))
mcols(gr1)$ref <- c("A", "C", "G")
mcols(gr1)$alt <- c("T", "G", "A")
whe1 <- WithinHostExperiment(
    assays = list(altFreq = matrix(c(0.10, 0.20, 0.30), ncol = 1)),
    rowRanges = gr1, colData = S4Vectors::DataFrame(sample_id = "S1"))
gr2 <- GRanges("seg", IRanges::IRanges(c(100, 200, 400), width = 1))
mcols(gr2)$ref <- c("A", "C", "T")
mcols(gr2)$alt <- c("T", "G", "C")
whe2 <- WithinHostExperiment(
    assays = list(altFreq = matrix(c(0.12, 0.18, 0.15), ncol = 1)),
    rowRanges = gr2, colData = S4Vectors::DataFrame(sample_id = "S1"))
cons_calls <- consensusISNV(list(ivar = whe1, lofreq = whe2))
cons_calls
mcols(SummarizedExperiment::rowRanges(cons_calls))

## ----callers-concordance------------------------------------------------------
calcCallerConcordance(list(ivar = whe1, lofreq = whe2))

## ----stats-permutation--------------------------------------------------------
set.seed(1)
permutationTest(div_df$shannon, div_df$group, n_perm = 999)

## ----viz-spectrum-------------------------------------------------------------
plotFrequencySpectrum(whe_qc)

## ----viz-pair-----------------------------------------------------------------
plotPairScatter(whe_qc, pairId = "pair_1")

## ----viz-dashboard------------------------------------------------------------
plotQCDashboard(whe_qc)

## ----viz-sfs------------------------------------------------------------------
plotSFS(list(donor = sfs_donor, recipient = sfs_recip))

## ----viz-trajectory-----------------------------------------------------------
plotFrequencyTrajectory(whe_l, hostId = "H1")

## ----case-import--------------------------------------------------------------
donor_file <- system.file("extdata", "MHM1143_rep1.tsv",
    package = "WithinHostExperiment")
recip_file <- system.file("extdata", "MHM1184_rep1.tsv",
    package = "WithinHostExperiment")
whe_sars2 <- readWithinHostTable(
    c(donor_file, recip_file), format = "ivar",
    colData = S4Vectors::DataFrame(
        sample_id = c("MHM1143", "MHM1184"),
        role = c("donor", "recipient"),
        pair_id = c("HH6", "HH6")))
whe_sars2

## ----case-seqlengths----------------------------------------------------------
GenomeInfoDb::seqlengths(whe_sars2) <- c(MN908947.3 = 29903L)
GenomeInfoDb::seqlengths(whe_sars2)

## ----case-qc------------------------------------------------------------------
whe_sars2 <- flagISNV(whe_sars2, ISNVFilter(minDepth = 100L, minFreq = 0.03))
qcSummary(whe_sars2)
div_sars2 <- calcDiversity(whe_sars2)
div_sars2

## ----case-plausibility--------------------------------------------------------
checkDiversityPlausibility(div_sars2, pathogen = "sars_cov_2")

## ----case-transmission--------------------------------------------------------
sv2 <- calcSharedVariants(whe_sars2, pairId = "HH6")
cat("Shared:", sv2$n_shared, "| donor-only:", sv2$n_donor_only,
    "| recipient-only:", sv2$n_recipient_only, "\n")
cat("Donor-only at adequate / low / unknown recipient depth:",
    sv2$n_donor_only_adequate_depth, "/", sv2$n_donor_only_low_depth, "/",
    sv2$n_donor_only_unknown_depth, "| recipient call failed QC:",
    sv2$n_donor_only_recipient_failed_qc, "\n")

## ----pathogen-flu-------------------------------------------------------------
flu_file <- system.file("extdata", "test_influenza_donor.tsv",
    package = "WithinHostExperiment")
whe_flu <- readWithinHostTable(flu_file, format = "ivar",
    colData = S4Vectors::DataFrame(sample_id = "flu_donor"))
GenomeInfoDb::seqlengths(whe_flu) <- c(HA = 1701L)
whe_flu <- flagISNV(whe_flu, adaptiveFilter("illumina"))
div_flu <- calcDiversity(whe_flu, indices = c("pi", "shannon"))
div_flu
checkDiversityPlausibility(div_flu, pathogen = "influenza")

## ----pathogen-hiv-------------------------------------------------------------
hiv_file <- system.file("extdata", "test_hiv_sample.tsv",
    package = "WithinHostExperiment")
whe_hiv <- readWithinHostTable(hiv_file, format = "ivar",
    colData = S4Vectors::DataFrame(sample_id = "hiv_sample"))
GenomeInfoDb::seqlengths(whe_hiv) <- c(pol = 1000L)
whe_hiv <- flagISNV(whe_hiv, adaptiveFilter("illumina"))
div_hiv <- calcDiversity(whe_hiv, indices = c("pi", "shannon"))
div_hiv
checkDiversityPlausibility(div_hiv, pathogen = "hiv")

## ----pathogen-tb--------------------------------------------------------------
tb_file <- system.file("extdata", "test_tb_sample.tsv",
    package = "WithinHostExperiment")
whe_tb <- readWithinHostTable(tb_file, format = "ivar",
    colData = S4Vectors::DataFrame(sample_id = "tb_sample"))
GenomeInfoDb::seqlengths(whe_tb) <- c(Rv0001 = 4411532L)
whe_tb <- flagISNV(whe_tb, adaptiveFilter("illumina"))
div_tb <- calcDiversity(whe_tb, indices = c("pi", "shannon"))
div_tb
checkDiversityPlausibility(div_tb, pathogen = "tb")

## ----pathogen-neutrality------------------------------------------------------
nt_flu <- calcNeutralityTests(whe_flu)
nt_hiv <- calcNeutralityTests(whe_hiv)
nt_tb  <- calcNeutralityTests(whe_tb)
data.frame(
    pathogen = c("Influenza", "HIV-1", "M. tuberculosis"),
    n_sites = c(nt_flu$n_sites, nt_hiv$n_sites, nt_tb$n_sites),
    tajimaD = c(nt_flu$tajimaD, nt_hiv$tajimaD, nt_tb$tajimaD),
    fusFs = c(nt_flu$fusFs, nt_hiv$fusFs, nt_tb$fusFs))

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