Worked example: SARS-CoV-2 immune escape
Loading packages
First we load some packages used in this example:
lapply(c("ggplot2", "patchwork"), library, character.only = TRUE)Conceptual framework and data
This example uses data from wastewater colected during the first and third COVID-19 waves in Argentina (BioProject PRJNA1183343; Manrique and Jones, n.d.). The first wave was driven by ancestral SARS-CoV-2 lineages, while the third was dominated by Omicron subvariants. Omicron is a highly adapted lineage to humans, including hypervariation in the receptor-binding domain (RBD) in response to immunity.
We will evidence the RBD relevance by an entropy analysis, and will identify key mutations on it.
We begin by building separate datasets for ancestral and Omicron lineages.
firstWave <- wWater[ wWater$wave == "first", ]
thirdWave <- wWater[ wWater$wave == "third", ]A 3% minor variant cutoff was already applied to these data by the variant caller. Following the same logic, here we will also filter out minor variants whose genotype matches that of the reference genome:
firstWave <- firstWave[ firstWave$alt_aa_freq <= 0.97, ]
thirdWave <- thirdWave[ thirdWave$alt_aa_freq <= 0.97, ]Entropy signatures
Now we infer the ancestral and Omicron signatures and compare them graphically.
ancestral <- getEntropySignature(firstWave)
omicron <- getEntropySignature(thirdWave)
# Compare signatures
anc_plot <- plot(ancestral) + ggtitle("Ancestral")
omi_plot <- plot(omicron) + ggtitle("Omicron")
anc_plot / omi_plot
We see that entropy at the S gene is relatively high in the Omicron sublineages.
Omicron sublineages in more detail
Now we look at how entropy is distributed along genome positions in the Omicron sublineages.
plot(omicron, chartType = "entroScan")
There are a handful of S positions that are responsible for the unusually high entropy.
omicron$Entropy$position[ omicron$Entropy$entropy > 0.3 ]
#> [1] 22882 22898 22917 23013 23040 23048 23055 23063Let’s look at the corresponding mutations:
showMutations(omicron, c(22882, 22898, 22917, 23013, 23040, 23048, 23055, 23063))
#> | snv | phenotype | abundance (%) |
#> |:-------:|:---------:|:-------------:|
#> | T22882G | S:N440K | 59 |
#> | G22898A | S:G446S | 60 |
#> | T22917G | S:L452R | 39 |
#> | A23013C | S:E484A | 66 |
#> | A23040G | S:Q493R | 66 |
#> | G23048A | S:G496S | 66 |
#> | A23055G | S:Q498R | 66 |
#> | A23063T | S:N501Y | 66 |The literature tells us that these SNVs affect the RBD (Lan et al. 2020).
Viral RBDs are targeted by the immune system because they are vital for cell binding and entry. There is solid observational and experimental evidence for S:N440K, S:G446S, S:L452R, S:E484A, S:Q493R, S:G496S, S:Q498R and S:N501Y (Guigon et al. 2022; Liang et al. 2022; Liu et al. 2021; Motozono et al. 2022; Rani et al. 2021; Starr, Greaney, Addetia, et al. 2021; Starr, Greaney, Dingens, et al. 2021; Weisblum et al. 2020; Zhang et al. 2022).
# Search PUBMED for more studies on these mutations:
library(rentrez)
#
# create search phrase
mutations <- showMutations(omicron, c(22882, 22898, 22917, 23013, 23040, 23048, 23055, 23063))
searchPhrase <- paste0(gsub("S:", "", mutations$phenotype, fixed = TRUE), "[TIAB]")
searchPhrase <- paste0("SARS-CoV-2[TITLE] AND spike [TIAB] AND (", paste0(searchPhrase, collapse = " OR "), ")")
#
# search
mySearch <- entrez_search(db = "pubmed", term = searchPhrase, retmax = 1000)
#
# Number of studies
mySearch$count
[1] 560
#We can do a formal analysis by assessHotSpot(). The RBD is located between genomic positions 22517 and 23186.
assessHotSpot(omicron, c(22517, 23186))
#>
#> Kruskal-Wallis rank sum test
#>
#> data: entropy by hot_spot
#> Kruskal-Wallis chi-squared = 8.9066, df = 1, p-value = 0.002841assessHotSpot(ancestral, c(22517, 23186))
#>
#> Kruskal-Wallis rank sum test
#>
#> data: entropy by hot_spot
#> Kruskal-Wallis chi-squared = 0.0055779, df = 1, p-value = 0.9405