Installation

# From Bioconductor (after acceptance):
if (!require("BiocManager", quietly = TRUE))
    install.packages("BiocManager")
BiocManager::install("WithinHostExperiment")

# Development version from GitHub:
BiocManager::install("CuiweiG/WithinHostExperiment")

Introduction

WithinHostExperiment is a Bioconductor infrastructure package for analysing intra-host single nucleotide variants (iSNV) from pathogen deep sequencing data.

The problem

Researchers performing within-host pathogen diversity analysis face a fragmented tooling landscape. Variant callers such as iVar, LoFreq, and Freebayes produce outputs in different formats. Quality control relies on ad hoc scripts. There is no standard Bioconductor container that organises within-host variant data with replicate information, donor-recipient pairing, and QC provenance.

What this package provides

WithinHostExperiment fills the infrastructure gap between upstream variant callers and downstream analysis tools:

A RangedSummarizedExperiment-based container
Replicate-aware, uncertainty-preserving QC (flag, don’t delete)
Within-host diversity quantification
Standardised export compatible with downstream bottleneck estimation tools

Relationship to existing tools

Tool	Role	Relationship
iVar, LoFreq, Freebayes	Variant calling	Upstream – we read output
deepSNV	Low-frequency variant detection	Upstream – complementary
ViralBottleneck	Bottleneck estimation	Downstream – simplified frame
QSutils	Quasispecies diversity	Parallel – different granularity

Why Bioconductor?

This package belongs in Bioconductor rather than CRAN for several reasons:

High-throughput genomic data. Within-host pathogen diversity analysis operates on deep sequencing data (BAM/VCF) – the core domain of Bioconductor.
Shared data structures. WithinHostExperiment extends RangedSummarizedExperiment and interoperates with VRanges, GRanges, and the VariantAnnotation ecosystem.
Complementary to existing packages. It provides an interoperability layer between upstream variant callers and downstream bottleneck estimation workflows through a standardised QC container.
Reproducibility. The flag-don’t-delete QC philosophy and full audit logging align with Bioconductor’s emphasis on reproducible workflows.

Quick start

library(WithinHostExperiment)

# Import iSNV data from VCF
vcf <- system.file("extdata", "test_donor.vcf",
    package = "WithinHostExperiment")
whe <- readWithinHost(vcf,
    colData = S4Vectors::DataFrame(sample_id = "donor"),
    caller = "ivar")
whe
#> WithinHostExperiment with 15 variant sites and 1 samples
#>   assays(5): altFreq, totalDepth, refCount, altCount, qcPass
#>   qcLog: 0 filtering step(s)
#>   transmission pairs: 0

Pipe-friendly workflow

All core functions take a WithinHostExperiment as their first argument, making them compatible with R’s native pipe:

library(GenomicRanges)
primers <- GRanges("test_segment", IRanges::IRanges(
    start = c(40, 880), end = c(65, 950)))

result <- whe |>
    flagISNV(ISNVFilter(minDepth = 200L, minFreq = 0.03)) |>
    flagPrimerSites(primers) |>
    passedISNV()

cat("Input:", nrow(whe), "sites ->",
    "After QC + primer masking:", nrow(result), "sites\n")
#> Input: 15 sites -> After QC + primer masking: 11 sites

Importing data

From VCF files

Each VCF file corresponds to one sample. Sample metadata is provided via colData:

vcf1 <- system.file("extdata", "test_donor.vcf",
    package = "WithinHostExperiment")
vcf2 <- system.file("extdata", "test_recipient.vcf",
    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"
)
whe
#> WithinHostExperiment with 20 variant sites and 2 samples
#>   assays(5): altFreq, totalDepth, refCount, altCount, qcPass
#>   qcLog: 0 filtering step(s)
#>   transmission pairs: 1

From iVar TSV

tsv <- system.file("extdata", "test_donor_ivar.tsv",
    package = "WithinHostExperiment")
whe_tsv <- readWithinHostTable(tsv, format = "ivar")
whe_tsv
#> WithinHostExperiment with 15 variant sites and 1 samples
#>   assays(5): altFreq, totalDepth, refCount, altCount, qcPass
#>   qcLog: 0 filtering step(s)
#>   transmission pairs: 0

The WithinHostExperiment object

The container extends RangedSummarizedExperiment:

# Assays: variant-level matrices (rows = sites, columns = samples)
SummarizedExperiment::assayNames(whe)
#> [1] "altFreq"    "totalDepth" "refCount"   "altCount"   "qcPass"

# rowRanges: genomic coordinates with ref/alt
head(SummarizedExperiment::rowRanges(whe))
#> GRanges object with 6 ranges and 2 metadata columns:
#>           seqnames    ranges strand |         ref         alt
#>              <Rle> <IRanges>  <Rle> | <character> <character>
#>   [1] test_segment        45      * |           C           T
#>   [2] test_segment        78      * |           A           C
#>   [3] test_segment       112      * |           A           T
#>   [4] test_segment       155      * |           A           G
#>   [5] test_segment       198      * |           A           G
#>   [6] test_segment       267      * |           C           A
#>   -------
#>   seqinfo: 1 sequence from an unspecified genome; no seqlengths

# colData: sample metadata
SummarizedExperiment::colData(whe)
#> DataFrame with 2 rows and 4 columns
#>             sample_id     host_id        role     pair_id
#>           <character> <character> <character> <character>
#> donor           donor          P1       donor      pair_1
#> recipient   recipient          P2   recipient      pair_1

# Transmission pairs (from colData role + pair_id)
transmissionPairs(whe)
#> DataFrame with 1 row and 3 columns
#>         donor   recipient     pair_id
#>   <character> <character> <character>
#> 1       donor   recipient      pair_1

Quality control

Accurate iSNV identification requires careful quality control. Roder et al. (2023, mBio) demonstrated that choice of variant caller and use of replicate sequencing have the most significant impact on iSNV discovery. WithinHostExperiment addresses both: it harmonises output from multiple callers (iVar, LoFreq, Freebayes) and provides replicate-aware QC as a first-class feature.

Standard QC with flagISNV

flagISNV() marks variants that fail quality criteria. It does not delete rows – it sets qcPass to FALSE. This preserves all data for potential re-analysis with different thresholds.

# Create a filter
filt <- ISNVFilter(minDepth = 200L, minFreq = 0.03)
filt
#> ISNVFilter object
#>   minDepth:      200 
#>   freq range:   [0.03, 0.5]
#>   minAltReads:   10 
#>   maxStrandBias: 10 
#>   minBaseQual:   20 
#>   minMapQual:    20 
#>   replicateConc: FALSE

# Apply (flag, don't delete)
whe_qc <- flagISNV(whe, filt)

# Check the audit log
qcLog(whe_qc)
#> DataFrame with 4 rows and 6 columns
#>               step   parameter       value n_flagged  n_passed
#>        <character> <character> <character> <integer> <integer>
#> 1     depth_filter    minDepth         200         1        39
#> 2  min_freq_filter     minFreq        0.03         0        39
#> 3  max_freq_filter     maxFreq         0.5         0        39
#> 4 alt_reads_filter minAltReads          10         0        39
#>             timestamp
#>           <character>
#> 1 2026-09-18T19:51:29
#> 2 2026-09-18T19:51:29
#> 3 2026-09-18T19:51:29
#> 4 2026-09-18T19:51:29

# Per-sample QC summary
qcSummary(whe_qc)
#> DataFrame with 2 rows and 7 columns
#>     sample_id total_sites passed_sites flagged_sites pass_rate median_depth
#>   <character>   <integer>    <integer>     <integer> <numeric>    <numeric>
#> 1       donor          15           14             1  0.933333         1600
#> 2   recipient          10           10             0  1.000000         1950
#>   median_freq
#>     <numeric>
#> 1        0.18
#> 2        0.25

# Extract only QC-passed variants
clean <- passedISNV(whe_qc)
cat("Before:", nrow(whe_qc), "sites\n")
#> Before: 20 sites
cat("After passedISNV:", nrow(clean), "sites\n")
#> After passedISNV: 20 sites

The two tables count different things. qcLog() counts cells of the site-by-sample matrix, so its 40 cells include those where a sample has no call at a site; a filter cannot fail a call that is absent, so those cells stay passed. qcSummary() counts only the calls each sample made. passedISNV() keeps a site if any sample passed it: the one flagged call, at position 812 in the donor, is absent from the recipient, so the site is kept and the row count does not change. The pipe example above, which also masks primer sites, is where rows are removed.

Replicate-aware QC

If your data includes technical replicates, flagReplicateDiscordance() identifies sites where replicates disagree:

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

whe_rep <- flagReplicateDiscordance(whe_rep, freqTolerance = 0.05)
qcLog(whe_rep)
#> DataFrame with 2 rows and 6 columns
#>                    step           parameter       value n_flagged  n_passed
#>             <character>         <character> <character> <integer> <integer>
#> 1 replicate_discordance       freqTolerance        0.05         6        26
#> 2 replicate_discordance requireBothDetected        TRUE         0        26
#>             timestamp
#>           <character>
#> 1 2026-09-18T19:51:29
#> 2 2026-09-18T19:51:29

Primer site filtering

For amplicon-based sequencing, variants in primer binding regions are unreliable and should be flagged:

library(GenomicRanges)
primers <- GRanges("test_segment", IRanges::IRanges(
    start = c(40, 880), end = c(65, 950)))
whe_primer <- flagPrimerSites(whe_qc, primers)
qcLog(whe_primer)
#> DataFrame with 5 rows and 6 columns
#>               step        parameter       value n_flagged  n_passed
#>        <character>      <character> <character> <integer> <integer>
#> 1     depth_filter         minDepth         200         1        39
#> 2  min_freq_filter          minFreq        0.03         0        39
#> 3  max_freq_filter          maxFreq         0.5         0        39
#> 4 alt_reads_filter      minAltReads          10         0        39
#> 5    primer_filter n_primer_regions           2         6        33
#>             timestamp
#>           <character>
#> 1 2026-09-18T19:51:29
#> 2 2026-09-18T19:51:29
#> 3 2026-09-18T19:51:29
#> 4 2026-09-18T19:51:29
#> 5 2026-09-18T19:51:30

Diversity analysis

div <- calcDiversity(whe_qc, genomeLength = 1000L)
div
#> DataFrame with 2 rows and 7 columns
#>     sample_id  richness   shannon         pi  watterson   simpson     chao1
#>   <character> <integer> <numeric>  <numeric>  <numeric> <numeric> <numeric>
#> 1       donor        14   6.83667 0.00446839 0.00175997    4.4656        14
#> 2   recipient        10   4.92373 0.00324286 0.00122661    3.2412        10

Transmission analysis

Shared variant statistics

sv <- calcSharedVariants(whe_qc, pairId = "pair_1")
cat("Shared:", sv$n_shared, "\n")
#> Shared: 5
cat("Donor-only:", sv$n_donor_only, "\n")
#> Donor-only: 9
cat("Recipient-only:", sv$n_recipient_only, "\n")
#> Recipient-only: 5

Quick bottleneck estimation

For a lightweight MLE estimate using the beta-binomial method:

nb <- quickBottleneck(whe_qc, pairId = "pair_1", maxNb = 50L)
cat("Nb MLE:", nb$Nb, "\n")
#> Nb MLE: 3
cat("95% CI: [", nb$ci[1], ",", nb$ci[2], "]\n")
#> 95% CI: [ 1 , 7 ]

asViralBottleneckInput() writes the donor and recipient frequencies of a pair into one frame, which is the input a beta-binomial estimator needs:

vb_input <- asViralBottleneckInput(whe_qc, pairId = "pair_1")
head(vb_input)
#>   pos donor_freq recipient_freq
#> 1  45       0.35           0.28
#> 2 112       0.08           0.00
#> 3 198       0.22           0.00
#> 4 267       0.45           0.39
#> 5 334       0.12           0.00
#> 6 401       0.06           0.00

Interoperability

WithinHostExperiment objects can be converted to standard Bioconductor and base R types for use with other packages.

Convert to data.frame

df <- as(whe_qc, "data.frame")
head(df)
#>          chrom position ref alt sample_id altFreq totalDepth refCount altCount
#> 1 test_segment       45   C   T     donor    0.35       1500      975      525
#> 2 test_segment      112   A   T     donor    0.08        800      736       64
#> 3 test_segment      198   A   G     donor    0.22       2200     1716      484
#> 4 test_segment      267   C   A     donor    0.45       3100     1705     1395
#> 5 test_segment      334   G   T     donor    0.12       1200     1056      144
#> 6 test_segment      401   C   T     donor    0.06        500      470       30
#>   qcPass
#> 1   TRUE
#> 2   TRUE
#> 3   TRUE
#> 4   TRUE
#> 5   TRUE
#> 6   TRUE

Convert to VRanges

For integration with VariantAnnotation workflows:

vr <- as(whe_qc, "VRanges")
vr
#> VRanges object with 25 ranges and 0 metadata columns:
#>            seqnames    ranges strand         ref              alt
#>               <Rle> <IRanges>  <Rle> <character> <characterOrRle>
#>    [1] test_segment        45      *           C                T
#>    [2] test_segment       112      *           A                T
#>    [3] test_segment       198      *           A                G
#>    [4] test_segment       267      *           C                A
#>    [5] test_segment       334      *           G                T
#>    ...          ...       ...    ...         ...              ...
#>   [21] test_segment       478      *           G                T
#>   [22] test_segment       523      *           A                G
#>   [23] test_segment       634      *           C                G
#>   [24] test_segment       650      *           T                A
#>   [25] test_segment       889      *           A                C
#>            totalDepth       refDepth       altDepth   sampleNames
#>        <integerOrRle> <integerOrRle> <integerOrRle> <factorOrRle>
#>    [1]           1500            975            525         donor
#>    [2]            800            736             64         donor
#>    [3]           2200           1716            484         donor
#>    [4]           3100           1705           1395         donor
#>    [5]           1200           1056            144         donor
#>    ...            ...            ...            ...           ...
#>   [21]            900            810             90     recipient
#>   [22]           2100           1365            735     recipient
#>   [23]           1900           1121            779     recipient
#>   [24]            400            380             20     recipient
#>   [25]           3200           2240            960     recipient
#>        softFilterMatrix
#>                <matrix>
#>    [1]                 
#>    [2]                 
#>    [3]                 
#>    [4]                 
#>    [5]                 
#>    ...              ...
#>   [21]                 
#>   [22]                 
#>   [23]                 
#>   [24]                 
#>   [25]                 
#>   -------
#>   seqinfo: 1 sequence from an unspecified genome; no seqlengths
#>   hardFilters: NULL

Export for downstream bottleneck tools

The frame shown above under “Quick bottleneck estimation” carries the positions and the two frequencies, which is what a beta-binomial estimator takes. ViralBottleneck itself builds its objects from per-base read counts through CreateTransmissionObject(), so prepare those from the altCount and refCount assays rather than from this frame.

Visualisation

plotFrequencySpectrum(whe_qc)

Allele frequency spectrum of the QC-flagged example pair.

plotPairScatter(whe_qc, pairId = "pair_1")

Donor against recipient allele frequency at each site.

plotQCDashboard(whe_qc)

Four-panel QC dashboard for the example object.

Reporting summaries

calcDiversity() and quickBottleneck() give the figures a report on a transmission pair usually needs:

div <- calcDiversity(whe_qc, genomeLength = 1000L)
cat("Shannon entropy:", round(div$shannon[1], 2), "\n")
#> Shannon entropy: 6.84
cat("Nucleotide diversity (pi):", format(div$pi[1], digits = 3), "\n")
#> Nucleotide diversity (pi): 0.00447
cat("Richness:", div$richness[1], "iSNV sites\n")
#> Richness: 14 iSNV sites

nb <- quickBottleneck(whe_qc, pairId = "pair_1", maxNb = 50L)
cat("Bottleneck Nb:", nb$Nb, "(95% CI:", nb$ci[1], "-", nb$ci[2], ")\n")
#> Bottleneck Nb: 3 (95% CI: 1 - 7 )

The example files describe a synthetic 1,000-base segment, so these numbers show the form of the output and are not evidence that the estimators are accurate. validateDiversity(), used below, compares an estimate with the package’s own plausibility range for a pathogen: a sanity check, not a validation.

Case study: SARS-CoV-2 transmission pair

This section demonstrates the complete workflow on real SARS-CoV-2 iSNV data from Bendall et al. (2023) Nature Communications 14:272. We analyse household HH6 from Bendall et al. (2023). The study’s pair table lists this household in both directions, as HH6_A and HH6_B, so it does not establish who infected whom; the example follows HH6_A and treats MHM1143 as the donor for illustration.

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
#> WithinHostExperiment with 4 variant sites and 2 samples
#>   assays(5): altFreq, totalDepth, refCount, altCount, qcPass
#>   qcLog: 0 filtering step(s)
#>   transmission pairs: 1

QC and diversity

whe_sars2 <- flagISNV(whe_sars2,
    ISNVFilter(minDepth = 100L, minFreq = 0.03))
qcSummary(whe_sars2)
#> DataFrame with 2 rows and 7 columns
#>     sample_id total_sites passed_sites flagged_sites pass_rate median_depth
#>   <character>   <integer>    <integer>     <integer> <numeric>    <numeric>
#> 1     MHM1143           2            2             0       1.0         1714
#> 2     MHM1184           2            1             1       0.5         2667
#>   median_freq
#>     <numeric>
#> 1   0.0470586
#> 2   0.0328049

div_sars2 <- calcDiversity(whe_sars2, genomeLength = 29903L)
div_sars2
#> DataFrame with 2 rows and 7 columns
#>     sample_id  richness   shannon          pi   watterson   simpson     chao1
#>   <character> <integer> <numeric>   <numeric>   <numeric> <numeric> <numeric>
#> 1     MHM1143         2  0.378533 5.99618e-06 8.33587e-06 0.1791992         2
#> 2     MHM1184         1  0.156079 2.34206e-06 3.95021e-06 0.0700084         1

Transmission analysis

sv <- calcSharedVariants(whe_sars2, pairId = "HH6")
cat("Shared:", sv$n_shared, "/ Donor-only:", sv$n_donor_only,
    "/ Recipient-only:", sv$n_recipient_only, "\n")
#> Shared: 0 / Donor-only: 2 / Recipient-only: 1

Other pathogens

To show that the workflow is not specific to SARS-CoV-2, we run it on small synthetic example files for three further pathogens and check the diversity estimates against the plausibility ranges of validateDiversity(). The files are illustrative fixtures, not measurements, and the ranges are the package’s own heuristics, so these sections demonstrate the interface rather than validate the estimates.

Influenza A (H3N2 HA segment)

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"))
whe_flu <- flagISNV(whe_flu, adaptiveFilter("illumina"))
div_flu <- calcDiversity(whe_flu, genomeLength = 1701L,
    indices = c("pi", "shannon"))
cat("Influenza pi:", format(div_flu$pi[1], digits = 3), "\n")
#> Influenza pi: 0.00255
cat("Influenza Shannon:", round(div_flu$shannon[1], 2), "\n")
#> Influenza Shannon: 6.47
validateDiversity(div_flu, pathogen = "influenza")
#> DataFrame with 2 rows and 5 columns
#>        metric   observed expected_min expected_max  in_range
#>   <character>  <numeric>    <numeric>    <numeric> <logical>
#> 1     shannon 6.46662479        0e+00        2e+01      TRUE
#> 2          pi 0.00255148        1e-05        5e-03      TRUE

HIV-1 (pol region)

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"))
whe_hiv <- flagISNV(whe_hiv, adaptiveFilter("illumina"))
div_hiv <- calcDiversity(whe_hiv, genomeLength = 1000L,
    indices = c("pi", "shannon"))
cat("HIV pi:", format(div_hiv$pi[1], digits = 3), "\n")
#> HIV pi: 0.00663
cat("HIV Shannon:", round(div_hiv$shannon[1], 2), "\n")
#> HIV Shannon: 10.01
validateDiversity(div_hiv, pathogen = "hiv")
#> DataFrame with 2 rows and 5 columns
#>        metric    observed expected_min expected_max  in_range
#>   <character>   <numeric>    <numeric>    <numeric> <logical>
#> 1     shannon 10.01379070        0e+00        30.00      TRUE
#> 2          pi  0.00663095        1e-04         0.05      TRUE

Mycobacterium tuberculosis

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"))
whe_tb <- flagISNV(whe_tb, adaptiveFilter("illumina"))
div_tb <- calcDiversity(whe_tb, genomeLength = 4411532L,
    indices = c("pi", "shannon"))
cat("TB pi:", format(div_tb$pi[1], digits = 3), "\n")
#> TB pi: 6.4e-08
cat("TB Shannon:", round(div_tb$shannon[1], 2), "\n")
#> TB Shannon: 0.58
validateDiversity(div_tb, pathogen = "tb")
#> DataFrame with 2 rows and 5 columns
#>        metric    observed expected_min expected_max  in_range
#>   <character>   <numeric>    <numeric>    <numeric> <logical>
#> 1     shannon 5.81456e-01        0e+00        5e+00      TRUE
#> 2          pi 6.40488e-08        1e-09        1e-04      TRUE

Neutrality tests across pathogens

nt_flu <- calcNeutralityTests(whe_flu, genomeLength = 1701L)
nt_hiv <- calcNeutralityTests(whe_hiv, genomeLength = 1000L)
nt_tb  <- calcNeutralityTests(whe_tb,  genomeLength = 4411532L)
results <- data.frame(
    pathogen = c("Influenza", "HIV-1", "M. tuberculosis"),
    n_sites = c(nrow(whe_flu), nrow(whe_hiv), nrow(whe_tb)),
    tajimaD = c(nt_flu$tajimaD, nt_hiv$tajimaD, nt_tb$tajimaD),
    fusFs = c(nt_flu$fusFs, nt_hiv$fusFs, nt_tb$fusFs)
)
results
#>          pathogen n_sites    tajimaD      fusFs
#> 1       Influenza      12  2.3654771  0.9553573
#> 2           HIV-1      20  2.1215611 -0.7177965
#> 3 M. tuberculosis       3 -0.9116534 -1.7188556

These fixtures are synthetic, and the tuberculosis one has only three sites, so the signs above carry no biological meaning. In real within-host data a negative Tajima’s D or Fu’s Fs is often read as purifying selection or expansion after a bottleneck, and a positive value as balancing selection or population structure; either reading needs far more sites than these examples hold.

Advanced features

This section demonstrates additional analytical modules available in WithinHostExperiment. Each subsection provides a minimal working example.

Neutrality tests

Tajima’s D and Fu’s Fs detect departures from neutral evolution. The calcNeutralityTests() method caps the effective sample size at n = 100 to avoid deep-sequencing bias (see ?tajimaD).

nt <- calcNeutralityTests(whe_qc, genomeLength = 1000L)
nt
#> DataFrame with 2 rows and 3 columns
#>     sample_id   tajimaD     fusFs
#>   <character> <numeric> <numeric>
#> 1       donor   1.82676 0.0382499
#> 2   recipient   1.77660 0.6937062

GFF3 gene annotation

Annotate variant sites with gene names from a GFF3 file. This automates the positional overlap needed for dndsWithinHost():

gff_file <- tempfile(fileext = ".gff3")
writeLines(c(
    "##gff-version 3",
    "test_segment\t.\tgene\t1\t500\t.\t+\t.\tName=geneA",
    "test_segment\t.\tgene\t501\t1000\t.\t+\t.\tName=geneB"
), gff_file)
whe_ann <- annotateFromGFF(whe_qc, gff_file)
table(S4Vectors::mcols(
    SummarizedExperiment::rowRanges(whe_ann))$GFF_FEATURE)
#> 
#> geneA geneB 
#>    11     9

For codon-aware annotation (amino acid changes), use annotateCodonChange() with a GFF3 and reference FASTA:

ref_fa <- system.file("extdata", "test_ref.fa",
    package = "WithinHostExperiment")
if (file.exists(ref_fa)) {
    gff_cds <- tempfile(fileext = ".gff3")
    writeLines(c("##gff-version 3",
        "test_segment\t.\tCDS\t1\t999\t.\t+\t0\tgene=testGene"),
        gff_cds)
    whe_codon <- annotateCodonChange(whe_qc, gff_cds, ref_fa)
    mc <- S4Vectors::mcols(
        SummarizedExperiment::rowRanges(whe_codon))
    table(mc$AA_CLASS)
}
#> 
#>      Nonsense Nonsynonymous    Synonymous 
#>             5            12             3

For complex genomes with overlapping reading frames or non-standard codon usage, consider SnpEff or ANNOVAR upstream.

Selection analysis

Within-host dN/dS

dndsWithinHost() counts synonymous and nonsynonymous iSNVs, divides each by the synonymous and nonsynonymous sites of the coding sequences (Nei and Gojobori 1986) and applies the Jukes-Cantor correction. It uses the codon annotation from annotateCodonChange() together with the same GFF3 and reference FASTA:

if (file.exists(ref_fa)) {
    dnds <- dndsWithinHost(whe_codon, gff = gff_cds, refFasta = ref_fa)
    dnds[, c("sample_id", "nS", "nN", "S_sites", "N_sites", "dNdS")]
}
#> DataFrame with 2 rows and 6 columns
#>     sample_id        nS        nN   S_sites   N_sites      dNdS
#>   <character> <integer> <integer> <numeric> <numeric> <numeric>
#> 1       donor         3        11   219.333   725.667   1.10937
#> 2   recipient         1         9   219.333   725.667   2.73465

Selection coefficients from frequency trajectories

For longitudinal data, estimate the selection coefficient from the logit-transformed frequency change:

freqs <- c(0.05, 0.10, 0.25, 0.50)
s_hat <- estimateSelectionCoefficient(freqs,
    times = c(0, 3, 7, 14))
s_hat
#> [1] 0.2490715 0.2746531 0.1569446

Multi-caller consensus

When data from multiple variant callers are available, create a consensus call set:

gr1 <- GenomicRanges::GRanges("seg",
    IRanges::IRanges(c(100, 200, 300), width = 1))
S4Vectors::mcols(gr1)$ref <- c("A", "C", "G")
S4Vectors::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 <- GenomicRanges::GRanges("seg",
    IRanges::IRanges(c(100, 200, 400), width = 1))
S4Vectors::mcols(gr2)$ref <- c("A", "C", "T")
S4Vectors::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 <- consensusISNV(list(ivar = whe1, lofreq = whe2))
cons
#> WithinHostExperiment with 2 variant sites and 1 samples
#>   assays(2): altFreq, qcPass
#>   qcLog: 0 filtering step(s)
#>   transmission pairs: 0

## Pairwise caller concordance
calcCallerConcordance(list(ivar = whe1, lofreq = whe2))
#> DataFrame with 1 row and 6 columns
#>       caller1     caller2  n_shared   n_only1   n_only2   jaccard
#>   <character> <character> <integer> <integer> <integer> <numeric>
#> 1        ivar      lofreq         2         1         1       0.5

Exact bottleneck estimation

exactBottleneck() offers three estimators beyond the approximate quickBottleneck(): "exact_bb", "wright_fisher" and "presence_absence". The presence-absence estimator is shown here.

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

Longitudinal frequency tracking

Track variant frequencies across time points within a host:

gr_l <- GenomicRanges::GRanges("seg1",
    IRanges::IRanges(c(100, 200), width = 1))
S4Vectors::mcols(gr_l)$ref <- c("A", "C")
S4Vectors::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)))

## Track frequencies
trackFrequency(whe_l)
#> DataFrame with 4 rows and 5 columns
#>       host_id variant_key timepoint frequency qc_passed
#>   <character> <character> <numeric> <numeric> <logical>
#> 1          H1  seg1:100:T         1      0.02      TRUE
#> 2          H1  seg1:200:G         1      0.10      TRUE
#> 3          H1  seg1:100:T         2      0.12      TRUE
#> 4          H1  seg1:200:G         2      0.15      TRUE

## Detect emerging variants
detectEmergingVariants(whe_l)
#> DataFrame with 1 row and 7 columns
#>       host_id variant_key time_from   time_to freq_from   freq_to delta_freq
#>   <character> <character> <numeric> <numeric> <numeric> <numeric>  <numeric>
#> 1          H1  seg1:100:T         1         2      0.02      0.12        0.1

SCI strand bias filtering

The SCI score normalises strand bias by the global strand coverage ratio, which suits amplicon data with unequal strand coverage:

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
#> [1]  0.09008968 -3.83136321  0.03744161

Adaptive technology-specific filtering

Pre-configured filter presets for different sequencing platforms:

adaptiveFilter("illumina")
#> ISNVFilter object
#>   minDepth:      100 
#>   freq range:   [0.03, 0.5]
#>   minAltReads:   10 
#>   maxStrandBias: 10 
#>   minBaseQual:   20 
#>   minMapQual:    20 
#>   replicateConc: FALSE
adaptiveFilter("ont")
#> ISNVFilter object
#>   minDepth:      60 
#>   freq range:   [0.1, 0.5]
#>   minAltReads:   10 
#>   maxStrandBias: Inf 
#>   minBaseQual:   20 
#>   minMapQual:    20 
#>   replicateConc: FALSE
adaptiveFilter("pacbio_hifi")
#> ISNVFilter object
#>   minDepth:      50 
#>   freq range:   [0.05, 0.5]
#>   minAltReads:   10 
#>   maxStrandBias: 10 
#>   minBaseQual:   20 
#>   minMapQual:    20 
#>   replicateConc: FALSE

Diversity modelling and statistical utilities

Fit linear models to diversity data and correct for multiple testing:

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
#>          term    estimate std_error    t_value    p_value
#> 1 (Intercept)  1.56666667 0.6085806  2.5742960 0.08218575
#> 2   timepoint -0.20000000 0.2635231 -0.7589466 0.50306833
#> 3      groupB  0.06666667 0.4303315  0.1549193 0.88672102

## Multiple testing correction
raw_p <- c(0.001, 0.04, 0.03, 0.5)
multitestCorrection(raw_p)
#> [1] 0.00400000 0.05333333 0.05333333 0.50000000

Plausibility check

Cross-check observed diversity against the package’s plausibility range for a given pathogen:

div_check <- S4Vectors::DataFrame(pi = 1e-4, shannon = 5.0)
validateDiversity(div_check, pathogen = "sars_cov_2")
#> DataFrame with 2 rows and 5 columns
#>        metric  observed expected_min expected_max  in_range
#>   <character> <numeric>    <numeric>    <numeric> <logical>
#> 1          pi     1e-04        1e-06        0.001      TRUE
#> 2     shannon     5e+00        0e+00       15.000      TRUE

Scope and limitations

This package operates at the level of individual single nucleotide variants (iSNVs). It does not perform:

Haplotype phasing or quasispecies reconstruction, which require read-level (BAM) analysis outside the scope of RangedSummarizedExperiment. For haplotype-aware analysis, consider CliqueSNV (Knyazev et al. 2021 Nucleic Acids Res 49:e102), ShoRAH (Zagordi et al. 2011 BMC Bioinformatics 12:119), or similar tools. Their output can be imported into WithinHostExperiment for downstream diversity and QC analysis.
Co-infection vs diversity distinction, which requires phylogenetic analysis of full haplotypes. The multi-caller consensus and replicate-concordance modules reduce artefacts but do not resolve this.
Stochastic selection coefficient estimation. The estimateSelectionCoefficient() function uses a deterministic logit-frequency model. For drift-aware inference, consider dedicated Wright-Fisher or diffusion-based frameworks.

Session info

sessionInfo()
#> R version 4.6.1 (2026-06-24)
#> Platform: x86_64-pc-linux-gnu
#> Running under: Ubuntu 26.04.1 LTS
#> 
#> Matrix products: default
#> BLAS:   /usr/lib/x86_64-linux-gnu/openblas-pthread/libblas.so.3 
#> LAPACK: /usr/lib/x86_64-linux-gnu/openblas-pthread/libopenblasp-r0.3.32.so;  LAPACK version 3.12.0
#> 
#> locale:
#>  [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C              
#>  [3] LC_TIME=en_US.UTF-8        LC_COLLATE=en_US.UTF-8    
#>  [5] LC_MONETARY=en_US.UTF-8    LC_MESSAGES=en_US.UTF-8   
#>  [7] LC_PAPER=en_US.UTF-8       LC_NAME=C                 
#>  [9] LC_ADDRESS=C               LC_TELEPHONE=C            
#> [11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C       
#> 
#> time zone: Etc/UTC
#> tzcode source: system (glibc)
#> 
#> attached base packages:
#> [1] stats4    stats     graphics  grDevices utils     datasets  methods  
#> [8] base     
#> 
#> other attached packages:
#> [1] GenomicRanges_1.65.4        Seqinfo_1.3.2              
#> [3] IRanges_2.47.5              S4Vectors_0.51.10          
#> [5] BiocGenerics_0.59.12        generics_0.1.4             
#> [7] WithinHostExperiment_0.99.0 rmarkdown_2.32             
#> 
#> loaded via a namespace (and not attached):
#>  [1] tidyselect_1.2.1            dplyr_1.2.1                
#>  [3] farver_2.1.2                blob_1.3.0                 
#>  [5] Biostrings_2.81.9           S7_0.2.2                   
#>  [7] bitops_1.1-0                fastmap_1.2.0              
#>  [9] RCurl_1.98-1.20             VariantAnnotation_1.59.4   
#> [11] GenomicAlignments_1.49.2    XML_3.99-0.24              
#> [13] digest_0.6.39               lifecycle_1.0.5            
#> [15] KEGGREST_1.53.6             RSQLite_3.53.3             
#> [17] magrittr_2.0.5              compiler_4.6.1             
#> [19] rlang_1.3.0                 sass_0.4.10                
#> [21] tools_4.6.1                 yaml_2.3.12                
#> [23] rtracklayer_1.73.0          knitr_1.52                 
#> [25] S4Arrays_1.13.0             labeling_0.4.3             
#> [27] bit_4.6.0                   curl_8.0.0                 
#> [29] DelayedArray_0.39.6         RColorBrewer_1.1-3         
#> [31] abind_1.4-8                 BiocParallel_1.47.0        
#> [33] withr_3.0.3                 sys_3.4.3                  
#> [35] grid_4.6.1                  ggplot2_4.0.3              
#> [37] scales_1.4.0                SummarizedExperiment_1.43.0
#> [39] cli_3.6.6                   crayon_1.5.3               
#> [41] otel_0.2.0                  httr_1.4.9                 
#> [43] rjson_0.2.23                BiocBaseUtils_1.15.1       
#> [45] DBI_1.3.0                   cachem_1.1.0               
#> [47] parallel_4.6.1              AnnotationDbi_1.75.2       
#> [49] XVector_0.53.0              restfulr_0.0.17            
#> [51] matrixStats_1.5.0           vctrs_0.7.3                
#> [53] Matrix_1.7-6                jsonlite_2.0.0             
#> [55] patchwork_1.3.2             bit64_4.8.6                
#> [57] GenomicFeatures_1.65.0      maketools_1.3.2            
#> [59] jquerylib_0.1.4             glue_1.8.1                 
#> [61] codetools_0.2-20            gtable_0.3.6               
#> [63] GenomeInfoDb_1.49.1         BiocIO_1.23.3              
#> [65] UCSC.utils_1.9.0            tibble_3.3.1               
#> [67] pillar_1.11.1               htmltools_0.5.9            
#> [69] BSgenome_1.81.1             R6_2.6.1                   
#> [71] evaluate_1.0.5              lattice_0.23-1             
#> [73] Biobase_2.73.2              png_0.1-9                  
#> [75] Rsamtools_2.29.0            cigarillo_1.3.1            
#> [77] memoise_2.0.1               bslib_0.12.0               
#> [79] SparseArray_1.13.2          xfun_0.61                  
#> [81] MatrixGenerics_1.25.0       buildtools_1.0.0           
#> [83] pkgconfig_2.0.3

Introduction to WithinHostExperiment

Installation

Introduction

The problem

What this package provides

Relationship to existing tools

Why Bioconductor?

Quick start

Pipe-friendly workflow

Importing data

From VCF files

From iVar TSV

The WithinHostExperiment object

Quality control

Standard QC with flagISNV

Replicate-aware QC

Primer site filtering

Diversity analysis

Transmission analysis

Shared variant statistics

Quick bottleneck estimation

Interoperability

Convert to data.frame

Convert to VRanges

Export for downstream bottleneck tools

Visualisation

Reporting summaries

Case study: SARS-CoV-2 transmission pair

QC and diversity

Transmission analysis

Other pathogens

Influenza A (H3N2 HA segment)

HIV-1 (pol region)

Mycobacterium tuberculosis

Neutrality tests across pathogens

Advanced features

Neutrality tests

GFF3 gene annotation

Selection analysis

Within-host dN/dS

Selection coefficients from frequency trajectories

Multi-caller consensus

Exact bottleneck estimation

Longitudinal frequency tracking

SCI strand bias filtering

Adaptive technology-specific filtering

Diversity modelling and statistical utilities

Plausibility check

Scope and limitations

Session info