Introduction

Intra-host viral diversity (or quasispecies) drives resilience to the immune responses and therapeutical interventions, while also fostering the emergence of new lineages (Andino and Domingo 2015; Eigen 1996; Jones 2025). The intraHostVariants dataset contains SARS-CoV-2 single nucleotide variants (SNVs) identified from autopsies of individuals with postmortem viral detection and no documented immunocompromise (Manrique et al. 2024).

Although MetaEntropy is tailored primarily for high-density wastewater data, it is also effective for quantifying and analyzing intra-host viral protein diversity, as demonstrated in this vignette. For specialized or comprehensive quasispecies characterization, dedicated tools such as QSutils (Guerrero-Murillo and Gregori i Font 2018) are recommended.

First load the required packages:

lapply(c("MetaEntropy", "ggplot2", "patchwork", "tidyr"),
       library, character.only = TRUE
)

Data Structure and Characteristics

intraHostVariants is a data.frame containing information on 213 SNVs inferred across lung, intestinal, cardiac, renal, and hepatic tissues. Its 16 columns include everything required by MetaEntropy, alongside experimental metadata and standard outputs typical of variant caller tables. The underlying bioinformatic workflow is detailed in the intraHostVariants documentation.

SARS-CoV-2 infections primarily cause respiratory illness, but the virus can spread systemically and replicate independently across different anatomical sites (Stein et al. 2022; Manrique et al. 2024). Consequently, individual tissues may harbor exclusive SNVs or share them at varying frequencies. This within-host diversity is well illustrated by intraHostVariants.

First, let’s look at the distribution of SNVs across anatomical compartments:

addmargins(
       table(intraHostVariants$case,
         intraHostVariants$organ
         ),
       margin = 2
)
#>      
#>       heart intestine kidney liver lung Sum
#>   c18    12        10      4     0    4  30
#>   c24    40        21     18     4    2  85
#>   c28     5        30     36    12    5  88
#>   c5      2         0      1     4    3  10

Two to five SNVs are recorded per respiratory sample, closely aligning with ranges reported in prior swab-based studies (Jones 2025). However, the total number of SNVs recorded per individual ranges from 10 to 88 across all tissues.

Next, let’s examine the distribution of SNV frequencies. First, we pivot the dataset to place the alternative amino acid frequencies into separate columns for each anatomical compartment:

wide_ihv <- tidyr::pivot_wider(intraHostVariants,
         names_from = organ,
         values_from = alt_aa_freq,
         # Use 0 (instead of NA) for undetected
         # mutations.
         values_fill = list(alt_aa_freq = 0)
)

We can now use wide_ihv to visualize these variant frequencies across tissues with parallel coordinates plots:

par(mfcol = c(4,1), mar = c(2, 4, 2, 2) + 0.1)
for(case in sort(unique(wide_ihv$case))){
    thisCase <- t(wide_ihv[wide_ihv$case == case, c(15, 16, 17, 18, 19)])
    matplot(thisCase, type = "b", pch = 19, col = rgb(0, 0, 0, 0.4), lty = 1,
        xaxt = "n", xlab = "",
        main = sub("c", "case ", case),
        ylab = "SNV frequency"
    )
    axis(side = 1, at = 1:5, labels = rownames(thisCase)
    )
}

This reveals that the viral quasispecies varies considerably between anatomical sites. For instance, the lung and heart in case 18 harbor three high-frequency SNVs that are absent from other organs:

intraHostVariants[intraHostVariants$case == "c18" & intraHostVariants$alt_aa_freq > 0.2 ,
                  c(1, 2, 10, 13, 14, 15, 16)]
#>    case organ protein aa_position ref_aa alt_aa alt_aa_freq
#> 12  c18 heart   nsp10         121      G      V    0.350514
#> 18  c18 heart   nsp15         145      E      E    0.985571
#> 38  c18  lung    nsp3         534      A      A    0.988260

Package Integration

Intra-host SNVs can be synonymous or non-synonymous, and the impact of non-synonymous mutations varies depending on whether the resulting amino acid substitution is conservative or radical. From a virological perspective, these phenotypic changes are particularly relevant because they are more likely to alter viral physiology. MetaEntropy is particularly useful here because it calculates amino acid entropy and allows residues to be classified based on their physicochemical properties, offering a nuanced view of functional protein diversity. We will now use MetaEntropy to evaluate amino acid diversity across anatomical sites within each individual.

First, we split the dataset into separate data.frame objects according to the experimental design strata (organ and case), storing them as a list:

strata <- split(intraHostVariants,
                    list(intraHostVariants$organ, intraHostVariants$case),
                    sep = "_", drop = TRUE
)

mn908947.3.ihv <- mn908947.3
# Get the nsp12 3' end and assign it to nsp12_end
nsp12_rows <- mn908947.3.ihv$CDS$protein %in% c("nsp12a", "nsp12b")
nsp12_end <- max(mn908947.3.ihv$CDS$end[nsp12_rows])
# Create an entry for nsp12 from nsp12a
mn908947.3.ihv$CDS$protein[mn908947.3.ihv$CDS$protein == "nsp12a"] <- "nsp12"
# Update the 3' end
mn908947.3.ihv$CDS$end[mn908947.3.ihv$CDS$protein == "nsp12"] <- nsp12_end
# dismiss the old nsp12b annotation
mn908947.3.ihv$CDS <- mn908947.3.ihv$CDS[mn908947.3.ihv$CDS$protein != "nsp12b", ]
# tidy up the environment
rm(nsp12_rows, nsp12_end)

profiles <- lapply(strata, function(df) {
               getEntropySignature(df, position = "POS", ref = "REF", alt = "ALT",
                           genome = mn908947.3.ihv
               )
})

heatmap_entropyProfiles(!!!profiles)

combined_entropy <- lapply(profiles, function(p) p$Entropy) |> do.call(rbind, args = _)
combined_entropy <- cbind(combined_entropy,
                          system = factor(ifelse(grepl("lung", rownames(combined_entropy)),
                                                       "respiratory system", "other"),
                                          levels = c("respiratory system", "other")
                                   )
                    )
# Plot in genomic order
combined_entropy$protein <- factor(
                                   combined_entropy$protein,
                                   levels = mn908947.3.ihv$CDS$protein
)
bp <- ggplot2::ggplot(data = combined_entropy, aes(x = protein, y = entropy)) +
    ggplot2::geom_boxplot(varwidth = T) +
    ggplot2::facet_wrap(~ system, ncol = 1, scales = "free_y") +
    ggplot2::scale_x_discrete(drop = FALSE) +
    ggplot2::theme_bw() +
    ggplot2::theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust = 1))
bp