rtt.Rmd

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

```{r message=FALSE, echo=F, warning=F}
require(MASS, quietly=T)

source("scripts/fit-rtt.R")
rtt1 <- fit.rtt(paste0(params$datadir, "/aligned_allSeqs_sample1.rtt.nwk"))
fit1 <- rtt1$fits

#exclusion list for the RTT lines
rttLineBlacklist <- c("other", "Recombinants")

# extract line segments
segments <- lapply(1:length(fit1), function(i) {
  pango <- names(fit1)[i]
  if (!(pango %in% rttLineBlacklist)){

      p <- predict(fit1[[i]])
      x <- fit1[[i]]$model$x
      list(pango=pango, x1=min(x), x2=max(x), y1=min(p), y2=max(p))
  }
})

#remove null value from segment list
segments <- segments[segments!="NULL"]

# generate interactive plot
pal <- VOCVOI$color
names(pal) <- VOCVOI$name

r2d3(data=toJSON(list(tips=rtt1$tips, palette=as.list(pal), fits=segments), 
                 auto_unbox=TRUE), 
     script="js/rtt.js", container="div", elementId="rtt-element")
```


```{r message=FALSE, echo=F, warning=F}
fit2 <- fit.rtt(paste0(params$datadir,"/aligned_allSeqs_sample2.rtt.nwk"))$fits
fit3 <- fit.rtt(paste0(params$datadir,"/aligned_allSeqs_sample3.rtt.nwk"))$fits
```


### Molecular clock estimates (based on three independent subsamples)

Here we show the estimate of the substitution rate for 3 independent subsamples of different variants of interest (VOI), with their 95% confidence interval. The average rate of substitution within VOI is given by a bamboo colored dashed line. For comparison, the average rate of substitution across all samples is much higher (grey line), indicating that about half of the substitutions occur through normal routes (i.e. not chronic infections) of transmission while the other half occur with the appearance of new, highly divergent VOI, likely due to different evolutionary processes occurring within chronic infections (see [Neher 2022](https://academic.oup.com/ve/article/8/2/veac113/6887176) for details). BA.2 (red) appears to have a higher rate of mutation because it includes highly divergent sub-groups (i.e., potential saltation lineages), including CM.* and sub-types related to BA.2.86.

```{r message=FALSE, echo=F, warning=F}

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

  
  avgGlobalRate <- mean(sec.frame[rownames(sec.frame) %in% c("global", "global1", "global2"),]$est)
  
  sec.frame<-sec.frame[sec.frame$Lineage != "Recombinants",]
  sec.frame<-sec.frame[sec.frame$Lineage != "other",]
    sec.frame<-sec.frame[sec.frame$Lineage != "global",]

  
  avgVOIRate <- mean(sec.frame$est)
  
  pal <- VOCVOI$color
  names(pal) <- VOCVOI$name
  pal["other"] <- "white"
  
  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=.7,
                  position=position_dodge(1)) +
    geom_hline(aes(yintercept = avgGlobalRate), color = "#6F8FAF", linetype = "dashed", size =1) +
    geom_hline(aes(yintercept = avgVOIRate), color = "#D2B04C", linetype = "twodash", size =1) +
    annotate("text", x=12, y=avgGlobalRate*1.05, label = paste0("Global rate (", round(avgGlobalRate,2), ")"), color = "#6F8FAF") +
    annotate("text", x=12, y=avgVOIRate*1.10, label = paste0("VOI rate (", round(avgVOIRate,2), ")"), color = "#D2B04C") +

    labs(y="Substitutions / Genome / Day",
         x="Lineage", fill="Subsample") 
}

```