Contributing authors:

Data analysis, code, and maintenance of this notebook: Carmen Lia Murall, Raphaël Poujol, Susanne Kraemer, Arnaud N’Guessan, Sarah Otto, Art Poon, Jesse Shapiro. Input and direction by other members of Pillar 6 (https://covarrnet.ca/our-team/#pillar-6 ), which include: Fiona Brinkman, Caroline Colijn, Jorg Fritz, Morgan Langille, Paul Gordon, Julie Hussin, Jeff Joy, William Hsiao, and Erin Gill.

Sequence collection, generation, release, and feedback on analyses: Canadian laboratories as part of the CPHLN and CanCOGeN are making these data publicly available and contribute feedback on analyses presented here. A complete list of lab authors is in this repository, and more details are below in the Acknowledgement section.

Introduction

This notebook is built to explore Canadian SARS-CoV-2 genomic and epidemiological data with the aim of investigating viral evolution and spread. It is for discussion with pillar 6’s team and for sharing with collaborators, e.g. PH labs. These analyses can spur further research within or across pillars, be used for reports (or data dashboards), support discussions with the science communication pillar for public dissemination, and enable code reuse by public health authorities/laboratories for their internal use.

Canadian genomic and epidemiological data will be regularly pulled from various public sources (see list below) to keep these analyses up-to-date. Only representations of aggregate data will be posted here.

Important caveats and disclaimers:
These analyses represent only a snapshot of SARS-CoV-2 evolution in Canada. Only some infections are detected by PCR testing, only some of those are sent for whole-genome sequencing, and not all sequences are posted to public facing repositories. Sequencing volumes and priorities have changed during the pandemic, and the sequencing strategy is typically a combination of prioritizing outbreaks, travellers, public health investigations, and random sampling for genomic surveillance.

For example, specific variants or populations might be preferentially sequenced at certain times in certain jurisdictions. When possible, these differences in sampling strategies are mentioned but they are not always known. With the arrival of the Omicron wave, many jurisdictions across Canada reached testing and sequencing capacity mid-late December 2021 and thus switched to targeted testing of priority groups (e.g., hospitalized patients, health care workers, and people in high-risk settings). Therefore, from this time onward, case counts are likely underestimated and the sequenced virus diversity is not necessarily representative of the virus circulating in the overall population.

Thus, interpretation of these plots and comparisons between health regions should be made with caution, considering that the data may not be fully representative. These analyses are subject to change given new data and updated lineage designations.

## 1. LOAD processed metadata of Canadian sequences (with latest pangolin, division, and full seq IDs)
#Download metadata from VirusSeq, put the date here:

VirusSeq_date="2022-08-31"
pangoversion="4.1.2 (Viral AI)"


# this can be made more compact for faster loading
meta <- read.csv(gzfile("data_needed/virusseq.metadata.csv.gz"))

# amend our PANGO lineage classifications with those from ViralAI
#vai <- read.csv("data_needed/viralai.csv")
#idx <- match(meta$isolate, vai$isolate)
# all(meta$isolate == vai$isolate[idx], na.rm=T)
#meta$lineage <- vai$lineage[idx]
#write.csv(meta, "data_needed/virusseq.metadata.csv", row.names=F)

meta$sample.collection.date <- as.Date(meta$sample.collection.date)
meta$province <- meta$geo_loc_name..state.province.territory.


# restore pangolearn calls overriden by scorpio - issue #46
#plearn <- gsub(".+assignment ([A-Z]+\\.[0-9.]*).*", "\\1",
#               meta$note[meta$lineage=="None"])
#plearn[which(grepl("_", plearn))] <- "None"
#meta$lineage[meta$lineage=="None"] <- plearn

#make a pango.group column
VOCVOI <- data.frame(
  name=c('Alpha', 'Beta', 'Gamma',
         'Delta',
         'Lambda', 'Omicron BA.1',
         'Omicron BA.2', 'Omicron BA.4', 'Omicron BA.5',
         'A.23.1', 'B.1.438.1','Recombinants'),
         #'Mu', 'A.23.1', 'B.1.438.1'),
  pangodesignation=c('B.1.1.7|B.1.1.7.*', 'B.1.351|B.1.351.*', 'B.1.1.28|B.1.1.28.*',
            'B.1.617.2|B.1.617.2.*',
            'B.1.1.1|B.1.1.1.*', 'B.1.1.529.1|B.1.1.529.1.*',
            'B.1.1.529.2|B.1.1.529.2.*', 'B.1.1.529.4|B.1.1.529.4.*', 'B.1.1.529.5|B.1.1.529.5.*',
            'A.23.1|A.23.1.*', 'B.1.438.1|B.1.438.1.*', 'X*'),
            #'B.1.621|B.1.621.*', 'A.23.1|A.23.1.*', 'B.1.438.1|B.1.438.1.*'),
  color=c('#B29C71', '#F08C3A', '#0000FF',
          '#A6CEE3',
          '#FFFF00', '#8B0000',
          '#FF0000', '#BC544B','#CC66FF',
          '#9AD378', '#3EA534','black')
          #'#FFFF00', '#9AD378', '#3EA534')
)

VOCVOI$pattern <- gsub("\\.","\\\\.",VOCVOI$pangodesignation)
VOCVOI$pattern <- gsub("^","^",VOCVOI$pattern)
VOCVOI$pattern <- gsub("\\|","|^",VOCVOI$pattern)
VOCVOI$pattern <- gsub("([^\\*])\\|","\\1$|",VOCVOI$pattern)
VOCVOI$pattern <- gsub("([^\\*])$","\\1$",VOCVOI$pattern)
VOCVOI$pattern <- gsub("([\\*])$",".*",VOCVOI$pattern)
VOCVOI$pattern <- gsub("([\\*])\\|",".*|",VOCVOI$pattern)

variants <- sapply(VOCVOI$pattern, function(p) 
  grepl(p, meta$rawlineage, perl=T))

# verify that every row matches either 0 or 1 patterns
# table(apply(variants, 1, sum))

meta$pango.group <- 'other'
meta$pango.group[apply(variants, 1, sum)==1] <- VOCVOI$name[unlist(apply(variants, 1, which))]
meta$pango.group <- as.factor(meta$pango.group)

#To detect unassigned lineages due to wrong pattern matching
#print(unique(meta[meta$pango.group=="other",]$lineage))
#pb="B.1.617.2.25.1.2"
#print(variants <- sapply(VOCVOI$pattern, function(p) 
#  grepl(p, pb, perl=T)))
#print(grepl(VOCVOI$pattern[4], pb, perl=T))

pal <- VOCVOI$color
names(pal) <- VOCVOI$name
pal["other"] <- 'grey'  # named character vector

## 2. LOAD epidemiological data (PHAC)


#from: https://health-infobase.canada.ca/covid-19/epidemiological-summary-covid-19-cases.html?stat=num&measure=total&map=pt#a2
epidataCANall <- read.csv(url("https://health-infobase.canada.ca/src/data/covidLive/covid19-download.csv"))
epidataCANall$date <- as.Date(epidataCANall$date)
epidataCANall$prname <- gsub('_', ' ', epidataCANall$prname)
epidate <- tail(epidataCANall,1)$date #download date

epidataCANall$previousvalue=0
#small loop to get the numtoday column from previous versions of this file from the cumulative cases
for(row in 1:nrow(epidataCANall)) {
  p=epidataCANall[row, "prname"]
  subdf=epidataCANall[which( (epidataCANall$date > epidataCANall[row, "date"] & epidataCANall$prname==p) ),]
  if(nrow(subdf)!=0){
    nextrow=which( (epidataCANall$date == min(subdf$date) & epidataCANall$prname==p))
    epidataCANall[nextrow,"previousvalue"]=epidataCANall[row, "totalcases"]
  }
}
epidataCANall$numtoday=epidataCANall$totalcases-epidataCANall$previousvalue

# for barplots
meta$week <- cut(meta$sample.collection.date, 'week')

SARS-CoV-2 in Canada

Variants in Canada over time

# focus on emergence of VoCs in Canada
meta1 <- meta[as.Date(meta$week) > as.Date('2021-01-01'), ]
meta1$week <- as.factor(as.character(meta1$week))
#meta1$pango.group <- meta1$jstrick 

dat <- lapply(unique(meta1$province), function(x) {
  as.data.frame.matrix(table(meta1$week[meta1$province==x],
                             meta1$pango.group[meta1$province==x]))
})
names(dat) <- unique(meta1$province)

tab <- as.data.frame.matrix(table(meta1$week, meta1$pango.group))
dat[['Canada']] <- tab

# pass colour legend to JavaScript layer
dat[['legend']] <- as.list(VOCVOI$color)
names(dat[['legend']]) <- VOCVOI$name
dat[['legend']]$other <- 'grey'

r2d3(data=toJSON(dat), script="js/barplot.js", container="div",
     elementId="barplot-element")

These plots show the changing composition of sequences for all Canadian data posted to the VirusSeq Portal according to PANGO lineage/variant designation (Pango version 4.1.2 (Viral AI)), up to 2022-08-27. Because sampling and sequencing procedures vary by region and time, this does not necessarily reflect the true composition of SARS-CoV-2 viruses in Canada over time.

From the beginning of the pandemic to the fall of 2021, Canadian sequences were mostly of the wildtype lineages (pre-VOCs). By the beginning of summer 2021, the VOCs Alpha and Gamma were the most sequenced lineages overall in Canada. The Delta wave grew late summer with Delta sublineages AY.25 and AY.27 constituting sizable proportions of this wave. Omicron arrived in November of 2021 and currently the three main lineages in Canada are BA.1, BA.1.1, and BA.2, with multiple sub-lineages within each clade, e.g., BA.2.12. See below for jurisdictional differences of these plots.

There are two PANGO lineages that have a Canadian origin and that predominately spread within Canada (with some exportations internationally): B.1.438.1 and B.1.1.176. Other lineages of historical interest in Canada:

Variants of current interest, due to their potential growth advantage, or mutations of potential functional significance: * BA.5.1 * BA.5.2 * BA.5.2.1 * BA.5.6 * Variants with mutation S:R346T (includes some BA.5.2.1 and BA.4.6)

Data notes by province

All analyses draw on the most recent publicly available viral sequence data on ViralSeq and should be interpreted with caution due to lags in reporting and sequencing priorities that can differ across provinces or territories. Note that the NCCID provides a timeline of Canadian events related to each variant: https://nccid.ca/covid-19-variants/.

BC

British Columbia

Provincial sequencing strategy includes a subset of representative positive samples and prioritized cases (outbreaks, long-term care, travel-related, vaccine escape, hospitalized). Additional up-to-date covid data for this province can be found here:
http://www.bccdc.ca/health-info/diseases-conditions/covid-19/data-trends

AB

Alberta

Additional up-to-date COVID data for this province can be found here:
https://www.alberta.ca/stats/covid-19-alberta-statistics.htm#variants-of-concern

MB

Manitoba

Additional up-to-date COVID data for this province can be found here:
https://geoportal.gov.mb.ca/apps/manitoba-covid-19/explore

QC

Quebec

Provincial random sequencing has been temporarily suspended as of Feb 8th, 2021. Quebec provides a list of updates on changes to screening and sequencing strategies, found here (in French): https://www.inspq.qc.ca/covid-19/donnees/variants#methodologie. Additiona up-to-date COVID data for this province can be found here:
https://www.inspq.qc.ca/covid-19/donnees/variants

NS

Nova Scotia

Additional up-to-date COVID data for this province can be found here:
https://experience.arcgis.com/experience/204d6ed723244dfbb763ca3f913c5cad

NB

New Brunswick

Additional up-to-date COVID data for this province can be found here:
https://experience.arcgis.com/experience/8eeb9a2052d641c996dba5de8f25a8aa (NB dashboard)

NL

Newfoundland and Labrador

Additional up-to-date COVID data for this province can be found here:
https://covid-19-newfoundland-and-labrador-gnl.hub.arcgis.com/

Canadian trees

Here we present a subsampled phylogenetic snapshot of SARS-CoV-2 genomes from Canada. 

### metadata and trees
MLtree_sub1 <- read.tree("./data_needed/sample1.rtt.nwk")
ttree_sub1 <- read.tree("./data_needed/sample1.timetree.nwk")
require(lubridate)
dateseq <- seq(ymd('2019-12-01'), ymd('2022-07-01'), by='3 month')

stopifnot(all(sort(MLtree_sub1$tip.label) == sort(ttree_sub1$tip.label)))

get.tipnames <- function(tip.label) {
sapply(tip.label, function(x) {
  tokens <- strsplit(x, "_")[[1]]
  ntok <- length(tokens)
  paste(tokens[1:(ntok-2)], collapse='_')
  })
}
MLtree_sub1$tip.label <- get.tipnames(MLtree_sub1$tip.label)
ttree_sub1$tip.label <- get.tipnames(ttree_sub1$tip.label)
#epis <- sapply(strsplit(tips,"\\|"),`[`,1)
#metasub1<-metaCANall[metaCANall$Accession_ID %in% epis,]
metasub1 <- meta[meta$fasta.header.name %in% ttree_sub1$tip.label, ]
#stopifnot(nrow(metasub1) == Ntip(ttree_sub1))

#scale to number of mutations
MLtree_sub1$edge.length <- MLtree_sub1$edge.length*29903

### needed for plotting
pango_frame <- meta[c("lineage", "pango.group")]
pango_col <- pango_frame %>% distinct()  #make a lookup table for all lineages to pango.group
listpangogp <- unique(metasub1$pango.group)
listpangogp <- sort(listpangogp)
# division colours for heatmap
df_PTs <- data.frame(divisions=metasub1$province) #make df for heatmaps
rownames(df_PTs) <- metasub1$Specimen.Collector.Sample.ID  #needs row names for heatmaps

#division colours
listPTs <- unique(metasub1$province)
listPTs <- listPTs[order(listPTs)]
#listPTs <- listPTs[-11]  # remove hubei
getpal2 <- colorRampPalette(brewer.pal(length(listPTs), "Spectral")) #"Set3
pal2 <- getpal2(length(listPTs))
names(pal2) <- listPTs
#packageVersion("ggtree")
### time tree
mrdate <- max(na.omit(metasub1$sample.collection.date))

p1 <- ggtree(ttree_sub1, mrsd = as.Date(mrdate)) %<+%
  metasub1[ , c("fasta.header.name", "pango.group"), 
            drop=FALSE] + 
  geom_tippoint(color="black", size=1) +
  geom_tippoint(aes(color = pango.group), size=1) + #, shape= lab
  scale_color_manual(breaks=listpangogp, values=pal) +
  #ggtitle("Time tree - Canada - downsample1")+
  coord_cartesian(clip = 'off') + 
  theme_tree2(plot.margin=margin(6, 40, 6, 6)) +
  guides(color = guide_legend(override.aes = list(size = 4) ) ) +
  scale_x_continuous(
    breaks = Date2decimal(dateseq),
    labels = format(dateseq, c("%b", "%Y", "%b", "%b"))) +  
  theme(legend.position = "right", legend.title = element_blank(), 
        legend.text = element_text(size =12))

plot1 <- p1  #don't plot provinces
#plot1 <- gheatmap(p1, df_PTs[, "divisions", drop = FALSE], width= 0.1, color = F, colnames = F, legend_title = F) + 
#  scale_fill_manual(breaks=listPTs, values = pal2)

### diversity ML tree
p2 <- ggtree(MLtree_sub1) %<+% 
  metasub1[ , c("fasta.header.name", "pango.group"), drop=FALSE] + 
  geom_tippoint(color="black", size=1) +
  geom_tippoint(aes(color = pango.group), size=1) + 
  scale_color_manual(breaks=listpangogp, values = pal) +
  #ggtitle("ML tree - Canada - subsample1, n = 8K")+
  coord_cartesian(clip = 'off') + 
  theme_tree2(plot.margin=margin(6, 40, 6, 6)) +
  guides(color = guide_legend(override.aes = list(size = 4) ) ) +
  theme(legend.position = "right", legend.title = element_blank(), 
        legend.text = element_text(size =12))

plot2 <- p2
#plot2 <- gheatmap(p2, df_PTs[, "divisions", drop = FALSE], width= 0.1, color = F, colnames = F, legend_title = F) +
#  scale_fill_manual(breaks=listPTs, values = pal2)

### omicron time tree
metasub_omi <- metasub1[grepl("B\\.1\\.1\\.529\\.", metasub1$rawlineage),]
                        

#omi_epis<-metasub_omi$Accession_ID
#ref<-"EPI_ISL_402124"
#all_omi<-c(omi_epis,ref)
all_omi <- metasub_omi$fasta.header.name
ttree_omi<-keep.tip(ttree_sub1, all_omi)

p3 <- ggtree(ttree_omi, mrsd = as.Date(mrdate)) %<+%
  metasub_omi[ , c("fasta.header.name", "pango.group"), 
               drop=FALSE] + 
  geom_tippoint(color="black", size=1) +
  geom_tippoint(aes(color = pango.group), size=1) + #, shape= lab
  scale_color_manual(breaks=listpangogp,values=pal) +
  #ggtitle("Time tree - Canada - downsample1")+
  coord_cartesian(clip = 'off') + 
  theme_tree2(plot.margin=margin(6, 40, 6, 6)) +
  guides(color = guide_legend(override.aes = list(size = 4) ) ) +
  scale_x_continuous(breaks = Date2decimal(dateseq),
                     labels = format(dateseq, c("%b", "%Y", "%b", "%b"))) +
  theme(legend.position = "right", legend.title = element_blank(), 
        legend.text = element_text(size =12))

plot3 <- gheatmap(p3, df_PTs[, "divisions", drop = FALSE], 
                  width= 0.1, color = F, colnames = F, 
                  legend_title = F) + 
  scale_fill_manual(breaks=listPTs, values = pal2)

### omicron diversity tree
MLtree_omi<-keep.tip(MLtree_sub1, all_omi)

p4 <- ggtree(MLtree_omi) %<+% 
  metasub_omi[ , c("fasta.header.name", "pango.group"), 
               drop=FALSE] + 
  geom_tippoint(color="black", size=1) +
  geom_tippoint(aes(color = pango.group), size=1) + 
  scale_color_manual(breaks=listpangogp, values = pal) +
  #ggtitle("ML tree - Canada - subsample1, n = 8K")+
  coord_cartesian(clip = 'off') + 
  theme_tree2(plot.margin=margin(6, 40, 6, 6)) +
  guides(color = guide_legend(override.aes = list(size = 4) ) ) +
  theme(legend.position = "right", legend.title = element_blank(), 
        legend.text = element_text(size =12))

plot4 <- p4

#plot4 <- gheatmap(p4, df_PTs[, "divisions", drop = FALSE], width= 0.1, color = F, colnames = F, legend_title = F) +
#  scale_fill_manual(breaks=listPTs, values = pal2)

# table(metasub1$pango.group)

Time Tree

Diversity Tree

Omicron diversity tree

Mutational composition of Omicron

Tabulation of the most predominant mutational changes in Omicron, with adjacent rows comparing the composition of Canadian sublineages to that sublineage globally.

Mutational profile of Omicron and its sublineages in Canada and globally for the most prevalent (>75%) point mutations in each category (based on all genomes available on VirusSeq on 2022-08-31).

Major Omicron sublineages in Canada

Here we take a closer look, sub-dividing the three major Omicron sub-lineages currently circulating in Canada.

sublineages_BA1 <- unique(meta$lineage[grepl("B\\.1\\.1\\.529\\.1$|B\\.1\\.1\\.529\\.1\\.", meta$rawlineage)]) # "BA.1"
sublineages_BA2 <- unique(meta$lineage[grepl("B\\.1\\.1\\.529\\.2$|B\\.1\\.1\\.529\\.2\\.", meta$rawlineage)])  # "BA.2" 
#nameBA3 <- unique(meta$lineage[grepl("B\\.1\\.1\\.529\\.3$|B\\.1\\.1\\.529\\.3\\.", meta$rawlineage)])  # "BA.3" 
sublineages_BA4 <- unique(meta$lineage[grepl("B\\.1\\.1\\.529\\.4$|B\\.1\\.1\\.529\\.4\\.", meta$rawlineage)])  # "BA.4" 
sublineages_BA5 <- unique(meta$lineage[grepl("B\\.1\\.1\\.529\\.5$|B\\.1\\.1\\.529\\.5\\.", meta$rawlineage)])  # "BA.5" 


source("scripts/subtype_plotter.R")

BA.1

Sublineage BA.1

The grey zone represent rare lineages (not in top 15) : BA.1.1.12(1), BA.1.1.3(1), BA.1.13.1(1), BA.1.15.3(1), BA.1.17.1(1), BD.1(1), BA.1.12(2), BA.1.8(2), BA.1.1.11(3), BA.1.1.9(3), BA.1.1.4(4), BA.1.1.8(4), BA.1.10(4), BA.1.14.2(4), BA.1.1.17(5), BA.1.14.1(5), BA.1.7(5), BA.1.16(9), BA.1.5(11), BA.1.1.13(12), BA.1.1.7(15), BC.2(18), BA.1.15.2(27), BA.1.1.15(42), BA.1.21(56), BA.1.19(100), BA.1.9(128), BA.1.13(162), BA.1.1.1(164), BA.1.15.1(177), BA.1.1.2(281)

BA.2

Sublineage BA.2

The grey zone represent rare lineages (not in top 15) : BA.2.10.3(1), BA.2.14(1), BA.2.24(1), BA.2.3.8(1), BA.2.38.3(1), BA.2.58(1), BA.2.63(1), BA.2.73(1), BA.2.75.3(1), BA.2.15(2), BA.2.2.1(2), BA.2.27(2), BA.2.3.9(2), BA.2.30(2), BA.2.4(2), BA.2.42(2), BA.2.45(2), BA.2.53(2), BA.2.71(2), BA.2.9.6(2), BA.2.11(3), BA.2.12.2(3), BA.2.29(3), BA.2.3.14(3), BA.2.34(3), BA.2.38.2(3), BA.2.54(3), BA.2.55(3), BA.2.68(3), BA.2.70(3), BA.2.9.4(3), BH.1(3), BA.2.50(4), BA.2.62(4), BA.2.81(4), BA.2.3.12(5), BA.2.3.13(5), BA.2.3.15(5), BA.2.44(5), BA.2.79(5), BA.2.41(6), BA.2.47(6), BA.2.75.2(6), BA.2.78(6), BA.2.16(7), BA.2.23.1(7), BA.2.3.7(7), BA.2.35(7), BA.2.64(7), BG.6(7), BA.2.38.1(8), BA.2.43(8), BA.2.60(8), BA.2.9.1(8), BA.2.61(9), BA.2.3.16(10), BA.2.3.5(10), BA.2.52(10), BA.2.72(10), BA.2.22(11), BA.2.9.2(11), BG.4(11), BA.2.26(12), BA.2.6(12), BA.2.49(13), BA.2.8(14), BA.2.3.1(15), BA.2.40.1(15), BA.2.5(15), BA.2.51(15), BA.2.75.1(15), BA.2.13.1(18), BA.2.17(20), BA.2.32(21), BA.2.48(23), BA.2.2(24), BA.2.74(24), BA.2.9.3(24), BG.5(25), BA.2.7(28), BA.2.31(35), BA.2.13(42), BA.2.9.5(47), BA.2.75(51), BA.2.3.17(55), BA.2.76(57), BA.2.36(63), BA.2.3.2(66), BA.2.56(79), BA.2.18(93), BA.2.1(95), BG.2(105), BA.2.10.1(112), BA.2.37(122), BA.2.3.11(123)

BA.4

Sublineage BA.4

The grey zone represent rare lineages (not in top 15) :

BA.5

Sublineage BA.5

The grey zone represent rare lineages (not in top 15) : BA.5.2.7(1), BA.5.6.1(1), BE.2(1), BF.14(1), BF.15(1), BF.3.1(1), BA.5.1.12(2), BA.5.1.8(2), BA.5.1.9(2), BA.5.2.4(2), BA.5.3.4(2), BF.12(2), BA.5.2.6(3), BA.5.9(3), BA.5.5.1(4), BF.11(4), BA.5.6.2(5), BA.5.3.3(6), BF.3(7), BA.5.1.5(8), BF.6(8), BA.5.1.4(10), BF.1.1(10), BA.5.3(15), BF.2(19), BF.7(24), BA.5.1.6(26), BA.5.10(29), BF.8(29), BA.5.3.2(31), BK.1(39), BA.5.1.3(50), BA.5.1.10(56), BA.5.2.3(65), BA.5.1.2(77), BF.4(77), BA.5.2.2(82), BE.3(93), BA.5.8(110)

Selection on Omicron

Of particular interest when there are newly arriving or emerging lineages is whether they have a selective advantage (and by how much), relative to the lineages already circulating within a population. Here we examine the major sublineages of Omicron currently in Canada and their relative rate of spread.

Currently, the major variants circulating are Omicron sublineages BA.2 (including BG), BA.4, and BA.5 (including BE and BF). BA.4 and BA.5 are spreading in Canada, as illustrated below. Left plot: y-axis is the proportion of sub-lineages BA.4 and BA.5 among Omicron; right plot: y-axis describes the logit function, log(freq(BA.4 or BA.5)/freq(BA.2)), which gives a straight line whose slope is the selection coefficient if selection is constant over time (see methods).

For comparison, Alpha had a selective advantage of 6%-11% per day over preexisting SARS-CoV-2 lineages, and Delta had a selective advantage of about 10% per day over Alpha.

Caveat: Selection coefficients are not estimated for Alberta, which is currently taking a variant-specific sequencing strategy, based on an initial PCR screen, which would skew estimates of selection. Canada-wide estimates, thus, do not include this province. Separate analyses are provided for those provinces with at least 20 BA.2 sequences in the database. These selection analyses must be interpreted with caution due to the potential for non-representative sampling, lags in reporting, and spatial heterogeneity in prevalence of different sublineages across Canada.

startpar2 <- list(p=0.5, s=0.1)
startpar3 <- list(p=c(0.4, 0.1), s=c(0.05, 0.1))
  
#color for each lineage
col <- c(pal["Omicron BA.4"], pal["Omicron BA.5"])

#Set a starting date
#Note that the startdate shouldn't be too much before both alleles become common
#or rare migration events that die off could throw off the estimation procedure 
#(so that the parameter estimates account for the presence of those alleles long in the past).
startdate<-as.Date("2022-05-01") #Using a later date with less sampling noise
source("scripts/plot_selection_estimator.R")


#####################################################
# tabs for displaying in notebook
#each PT tab should have curve plot and breakpoint plot side by side

Canada

Canada without Alberta

BC

British Columbia

SK

Saskatchewan

MB

Manitoba

ON

Ontario

QC

Quebec

East

East provinces

BA.2 sublineage selection

Here we show the trends of the various BA.2.* sublineages over time, relative to the frequency of BA.2 by itself (shown for sublineages with at least 50 (Canada) or 20 (provinces) cases). Proportions shown here are only among BA.2 (stricto) and the lineage illustrated. Note that these plots are not necessarily representative of trends in each province and that mixing of data from different provinces may lead to shifts in frequency that are not due to selection.

startpar2 <- list(p=0.5, s=0.05)
startdate<-as.Date("2022-01-01") 
sublineagedata <- meta %>% filter(
    lineage %in% unique(meta$lineage[grepl("B\\.1\\.1\\.529\\.2$|B\\.1\\.1\\.529\\.2\\.", meta$rawlineage)]), 
    geo_loc_name..state.province.territory. %in% unique(meta$geo_loc_name..state.province.territory.),
    !is.na(sample.collection.date),
    sample.collection.date >= startdate
    ) %>% group_by(sample.collection.date,geo_loc_name..state.province.territory.) %>% dplyr::count(lineage)
source("scripts/plot_selection_estimator_2type.R")
namereference <- "BA.2"


getsublineagelist <- function(mainlineage) {
  togrep=paste(mainlineage,"$|",mainlineage,".",sep="")
  list_of_lineage=unique((meta %>% filter(lineage %in% unique(meta$lineage[grepl(gsub("\\.", "\\\\.",togrep), meta$rawlineage)])))$lineage)
  return(list_of_lineage)
}

getsublineagesiblings <- function(mainlineage) {
  decompose=(strsplit(mainlineage, "\\."))[[1]]
  parent=sapply(list(head(decompose,-1)), paste, collapse = ".")
  grepmainlineage=paste(mainlineage,"$|",mainlineage,".",sep="")
  greparent=paste(parent,"|",parent,".",sep="")
  togrep=paste("^(?!",grepmainlineage,")(",greparent,")",sep="")
  list_of_lineage=unique((meta %>% filter(lineage %in% unique(meta$lineage[grepl(gsub("\\.", "\\\\.",togrep), meta$rawlineage, perl=T)])))$lineage)
  return(list_of_lineage)
}

l="B.1.1.529.5.2"

#getsublineagelist(l)
#getsublineagesiblings(l)
#print(getsublineagelist("B.1.1.529.2.3"))
#plot.selection(region=region, startdate=startdate, reference=namereference, mutants="BA.2.3", #startpar=startpar2) 
#plot.selection(region=region, startdate=startdate, reference=namereference, mutants="BA.2.3", startpar=startpar2) 
#plot.selection(region=region, startdate=startdate, reference=getsublineagesiblings("B.1.1.529.2.3"), mutants="BA.2.3", startpar=startpar2) 
#plot.selection(region=region, startdate=startdate, reference=namereference, mutants=getsublineagelist("B.1.1.529.2.3"), startpar=startpar2)

Canada

Canada (n>50)

BC

British Columbia (n>20)

Ontario

Ontario (n>20)

Quebec

Quebec (n>20)

Saskatchawan

Saskatchawan (n>20)

[1] “Not enough data available”

Manitoba

Manitoba (n>20)

[1] “Not enough data available”

Eastern Provinces

Eastern provinces (n>20)

BA.4 sublineage selection

Here we show the trends of the various BA.4# sublineages over time, relative to the frequency of BA.4 by itself (shown for sublineages with at least 50 (Canada) or 20 (provinces) cases). Proportions shown here are only among BA.2 (stricto) and the lineage illustrated. Note that these plots are not necessarily representative of trends in each province and that mixing of data from different provinces may lead to shifts in frequency that are not due to selection.

Canada

Canada (n>50)

BC

British Columbia (n>20)

[1] “Not enough data available”

Ontario

Ontario (n>20)

Quebec

Quebec (n>20)

Saskatchawan

Saskatchawan (n>20)

[1] “Not enough data available”

Manitoba

Manitoba (n>20)

[1] “Not enough data available”

Eastern Provinces

Eastern provinces (n>20)

[1] “Not enough data available”

BA.5 sublineage selection

Here we show the trends of the various BA.5# sublineages over time, relative to the frequency of BA.5.1 by itself (shown for sublineages with at least 50 (Canada) or 20 (provinces) cases). Proportions shown here are only among BA.5.1 (stricto) and the lineage illustrated. Note that these plots are not necessarily representative of trends in each province and that mixing of data from different provinces may lead to shifts in frequency that are not due to selection.

Canada

Canada (n>50)

BC

British Columbia (n>20)

[1] “Not enough data available”

Ontario

Ontario (n>20)

Quebec

Quebec (n>20)

Saskatchawan

Saskatchawan (n>20)

[1] “Not enough data available”

Manitoba

Manitoba (n>20)

[1] “Not enough data available”

Eastern Provinces

Eastern provinces (n>20)

[1] “Not enough data available”

Root-to-tip analyses

The slope of root-to-tip plots over time provide an estimate of the substitution rate. A lineage with a steeper positive slope than average for SARS-CoV-2 is accumulating mutations at a faster pace, while a lineage that exhibits a jump up (a shift in intercept but not slope) has accumulated more than expected numbers of mutations in a transient period of time (similar to what we saw with Alpha when it first appeared in the UK).

source("scripts/fit-rtt.R")
fit1 <- fit.rtt("./data_needed/sample1.rtt.nwk", plot=TRUE)

fit2 <- fit.rtt("data_needed/sample2.rtt.nwk", plot=FALSE)
fit3 <- fit.rtt("data_needed/sample3.rtt.nwk", plot=FALSE)

Molecular clock estimates (based on three independent subsamples)

if(!is.null(fit1)){
est1 <- get.ci(fit1); est1$rep <- 'Rep1'
est2 <- get.ci(fit2); est2$rep <- 'Rep2'
est3 <- get.ci(fit3); est3$rep <- 'Rep3'
sec.frame <- rbind(est1, est2, est3)
sec.frame$est[sec.frame$est < 0] <- 0
sec.frame$lower.95[sec.frame$lower.95 < 0] <- 0
sec.frame=sec.frame[sec.frame$Lineage != "Recombinants",]
sec.frame=sec.frame[sec.frame$Lineage != "others",]

pal <- VOCVOI$color
names(pal) <- VOCVOI$name

ggplot(sec.frame, aes(x=Lineage, y=est, group=rep)) + 
  geom_bar(stat="identity", color="black", aes(fill=Lineage), position='dodge') + 
  scale_fill_manual(values=pal) + 
  theme(axis.text.x = element_text(size=9, angle=45, hjust=1, vjust=0.95),
        legend.position='none', panel.grid.major=element_line(colour="grey90")) + 
  geom_errorbar(aes(ymin=lower.95, ymax=upper.95), width=.2,
                position=position_dodge(1)) +
  labs(y="Substitutions / Genome / Day",
       x="Lineage", fill="Subsample")
}

Future development

We are in the process of adding or would like to develop code for some of the following analyses:

  • dN/dS (by variant and by gene/domains)
  • Tajima’s D over time
  • clustering analyses
  • genomically inferred epidemiological parameters: R0, serial interval, etc.

With anonymized data on vaccination status, severity/outcome, reason for sequencing (e.g., outbreak, hospitalization, or general sampling), and setting (workplace, school, daycare, LTC, health institution, other), we could analyze genomic characteristics of the virus relative to the epidemiological and immunological conditions in which it is spreading and evolving. Studies on mutational correlations to superspreading events, vaccination status, or comparisons between variants would allow us to better understand transmission and evolution in these environments.

Methodology

Genome data and metadata are sourced from the Canadian VirusSeq Data Portal. Pangolin lineage assignments are generated using the pangoLEARN algorithm. Source code for generating this RMarkdown notebook can be found in [https://github.com/CoVaRR-NET/duotang].

Trees

Phylogenetic trees

Canadian genomes were obtained from the VirusSeq data on the 2022-08-31 and down-sampled to two genomes per lineage, province and month before October 2021, and five genomes per lineage, province and month after October 2021 (about 10,000 genomes in total). We used a Python wrapper of minimap2 (version 2.17) to generate multiple sequence alignments for these genome samples. A maximum likelihood (ML) tree was reconstructed from each alignment using the COVID-19 release of IQ-TREE (version 2.2.0). Outliers were identified in by root-to-tip regression using the R package ape and removed from the dataset. TreeTime was used to reconstruct a time-scaled tree under a strict molecular clock model. The resulting trees were plotted with ggtree.

Mutational composition

Mutation composition graph

We extracted mutation frequencies from unaligned genomes using a custom Python wrapper of minimap2 (version 2.17). These data were supplemented with genomic data and metadata from the NCBI GenNank database, curated by the Nextstrain development team. We used these outputs to generate mutational graphs reporting mutations seen in at least 75% of sequences in the respective variants of concern in Canada. Bars are colored by substitution type, and the corresponding amino acid changes are shown. Genomic position annotations were generated in Python using SnpEFF.

Selection

Selection Coefficents

To estimate selection, we used standard likelihood techniques. In brief, sublineages of interest were prespecified (e.g., BA.1, BA.1.1, BA.2) and counts by day tracked over time. If selection were constant over time, the frequency of sub-type \(i\) at time \(t\) would be expected to rise according to \[p_i(t) = \frac{p_i(0) \exp(s_i t)}{\sum_j p_j(0) \exp(s_j t)},\] where \(s_i\) is the selection coefficient favouring sub-type \(i\). A selection coefficient of \(s_i=0.1\) implies that sub-type \(i\) is expected to rise from 10% to 90% frequency in 44 days (in \(4.4./s_i\) days for other values of \(s_i\)).

At any given time \(t\), the probability of observing \(n_i\) sequences of sublineage \(i\) is multinomially distributed, given the total number of sequences from that day and the frequency of each \(p_i(t)\). Consequently, the likelihood of seeing the observed sequence data over all times \(t\) and over all sublineages \(j\) is proportional to \[L = \prod_t \prod_j p_i(t)^{n_i(t)}.\]

The BBMLE package in R was used to maximize the likelihood of the observed data (using the default optimization method, optim). For each selection coefficient, 95% confidence intervals were obtained by profile likelihood (using uniroot).

Graphs illustrating the rise in frequency of a variant over time are shown (left panels), with the area of each dot proportional to the number of sequences. 95% confidence bands were obtained by randomly drawing 10,000 sets of parameters (\(p_i\) and \(s_i\) for each sub-type) using RandomFromHessianOrMCMC, assuming a multi-normal distribution around the maximum likelihood point (estimated from the Hessian matrix, Pawitan 2001). At each point in time, the 2.5%-97.5% range of values for \(p_i(t)\) are then shown in the confidence bands.

Logit plots (right panels) show \[ln(\frac{p_i(t)}{p_{ref}(t)})\] relative to a given reference genotype (here BA.1), which gives a line whose slope is the strength of selection \(s_i\). Changes in slope indicate changes in selection on a variant (e.g., see Otto et al.).

These estimates of selection ignore heterogeneity within provinces and may be biased by the arrival of travel-related cases while frequencies are very low. Sampling strategies that oversample clustered cases (e.g., sequencing outbreaks) will introduce additional variation beyond the multinomial expectation, but these should lead to one-time shifts in frequency rather than trends over time. Provinces with sampling strategies that are variant specific are removed, unless explicit information about the variant frequencies is available.

Rates

Root-to-tip estimates of substitution rate

Maximum likelihood tree (IQ-TREE) processed with root-to-tip regression and plotting in R.

List of useful tools

Collect a list of bioinformatics, phylogenetic, and modelling tools that are useful for SARS-CoV-2 analyses:

Acknowledgements and sources

We thank all the authors, developers, and contributors to the VirusSeq database for making their SARS-CoV-2 sequences publicly available. We especially thank the Canadian Public Health Laboratory Network, academic sequencing partners, diagnostic hospital labs, and other sequencing partners for the provision of the Canadian sequence data used in this work. Genome sequencing in Canada was supported by a Genome Canada grant to the Canadian COVID-19 Genomic Network (CanCOGeN).

We gratefully acknowledge all the Authors, the Originating laboratories responsible for obtaining the specimens, and the Submitting laboratories for generating the genetic sequence and metadata and sharing via the VirusSeq database, on which this research is based.

  • The Canadian VirusSeq Data Portal (https://virusseq-dataportal.ca) We wish to acknowledge the following organisations/laboratories for contributing data to the Portal: Canadian Public Health Laboratory Network (CPHLN), CanCOGGeN VirusSeq, Saskatchewan - Roy Romanow Provincial Laboratory (RRPL), Nova Scotia Health Authority, Alberta ProvLab North (APLN), Queen’s University / Kingston Health Sciences Centre, National Microbiology Laboratory (NML), Institut National de Sante Publique du Quebec (INSPQ), BCCDC Public Health Laboratory, Public Health Ontario (PHO), Newfoundland and Labrador - Eastern Health, Unity Health Toronto, Ontario Institute for Cancer Research (OICR), Provincial Public Health Laboratory Network of Nova Scotia, Centre Hospitalier Universitaire Georges L. Dumont - New Brunswick, and Manitoba Cadham Provincial Laboratory. Please see the complete list of laboratories included in this repository.

  • Public Health Agency of Canada (PHAC) / National Microbiology Laboratory (NML) - (https://health-infobase.canada.ca/covid-19/epidemiological-summary-covid-19-cases.html)

  • Various provincial public health websites (e.g. INSPQ https://www.inspq.qc.ca/covid-19/donnees/)

Session info

The version numbers of all packages in the current environment as well as information about the R install is reported below.

Hide

Show

sessionInfo()
## R version 4.2.0 (2022-04-22)
## Platform: x86_64-pc-linux-gnu (64-bit)
## Running under: Ubuntu 18.04.6 LTS
## 
## Matrix products: default
## BLAS:   /usr/lib/x86_64-linux-gnu/blas/libblas.so.3.7.1
## LAPACK: /usr/lib/x86_64-linux-gnu/lapack/liblapack.so.3.7.1
## 
## locale:
##  [1] LC_CTYPE=en_CA.UTF-8       LC_NUMERIC=C              
##  [3] LC_TIME=en_CA.UTF-8        LC_COLLATE=en_CA.UTF-8    
##  [5] LC_MONETARY=en_CA.UTF-8    LC_MESSAGES=en_CA.UTF-8   
##  [7] LC_PAPER=en_CA.UTF-8       LC_NAME=C                 
##  [9] LC_ADDRESS=C               LC_TELEPHONE=C            
## [11] LC_MEASUREMENT=en_CA.UTF-8 LC_IDENTIFICATION=C       
## 
## attached base packages:
## [1] stats4    stats     graphics  grDevices utils     datasets  methods  
## [8] base     
## 
## other attached packages:
##  [1] HelpersMG_5.4      rlang_1.0.2        coda_0.19-4        lme4_1.1-29       
##  [5] Matrix_1.4-1       bbmle_1.0.25       jsonlite_1.8.0     r2d3_0.2.6        
##  [9] MASS_7.3-57        viridis_0.6.2      viridisLite_0.4.0  colorspace_2.0-3  
## [13] RColorBrewer_1.1-3 ggtree_3.4.0       phytools_1.0-3     maps_3.4.0        
## [17] tidytree_0.3.9     phylotools_0.2.2   ape_5.6-2          treeio_1.20.0     
## [21] lubridate_1.8.0    knitr_1.39         forcats_0.5.1      stringr_1.4.0     
## [25] dplyr_1.0.9        purrr_0.3.4        readr_2.1.2        tidyr_1.2.0       
## [29] tibble_3.1.7       ggplot2_3.3.6      tidyverse_1.3.1   
## 
## loaded via a namespace (and not attached):
##  [1] minqa_1.2.4             ellipsis_0.3.2          fs_1.5.2               
##  [4] aplot_0.1.6             rstudioapi_0.13         farver_2.1.0           
##  [7] fansi_1.0.3             mvtnorm_1.1-3           xml2_1.3.3             
## [10] codetools_0.2-18        splines_4.2.0           mnormt_2.1.0           
## [13] nloptr_2.0.3            broom_0.8.0             dbplyr_2.2.0           
## [16] compiler_4.2.0          httr_1.4.3              backports_1.4.1        
## [19] assertthat_0.2.1        fastmap_1.1.0           lazyeval_0.2.2         
## [22] cli_3.3.0               htmltools_0.5.2         tools_4.2.0            
## [25] igraph_1.3.2            gtable_0.3.0            glue_1.6.2             
## [28] clusterGeneration_1.3.7 fastmatch_1.1-3         Rcpp_1.0.8.3           
## [31] cellranger_1.1.0        jquerylib_0.1.4         vctrs_0.4.1            
## [34] nlme_3.1-157            xfun_0.31               rvest_1.0.2            
## [37] lifecycle_1.0.1         phangorn_2.8.1          scales_1.2.0           
## [40] hms_1.1.1               parallel_4.2.0          expm_0.999-6           
## [43] yaml_2.3.5              gridExtra_2.3           ggfun_0.0.6            
## [46] yulab.utils_0.0.4       sass_0.4.1              bdsmatrix_1.3-6        
## [49] stringi_1.7.6           highr_0.9               plotrix_3.8-2          
## [52] boot_1.3-28             pkgconfig_2.0.3         evaluate_0.15          
## [55] lattice_0.20-45         labeling_0.4.2          patchwork_1.1.1        
## [58] htmlwidgets_1.5.4       tidyselect_1.1.2        magrittr_2.0.3         
## [61] R6_2.5.1                generics_0.1.2          combinat_0.0-8         
## [64] DBI_1.1.2               pillar_1.7.0            haven_2.5.0            
## [67] withr_2.5.0             scatterplot3d_0.3-41    modelr_0.1.8           
## [70] crayon_1.5.1            utf8_1.2.2              tzdb_0.3.0             
## [73] rmarkdown_2.14          grid_4.2.0              readxl_1.4.0           
## [76] reprex_2.0.1            digest_0.6.29           numDeriv_2016.8-1.1    
## [79] gridGraphics_0.5-1      munsell_0.5.0           ggplotify_0.1.0        
## [82] bslib_0.3.1             quadprog_1.5-8