Quantification of the spread of SARS-CoV-2 variant B.1.1.7 in Switzerland

Author links open overlay panel, , , , , , , , , , , , , , , , , , , …

https://www.sciencedirect.com/science/article/pii/S1755436521000335?via%3Dihub

2.2. Statistical inference

We fit a logistic model to the frequency of B.1.1.7 samples per day to estimate the logistic growth rate a and the sigmoid’s midpoint t0. From that, we derive an estimate of the transmission fitness advantage of B.1.1.7 under a continuous (fc) and a discrete (fd) model. Each model could plausibly describe the actual dynamics, so we present results from both for comparison. Further, we estimate the reproductive number R for the B.1.1.7 and non-B.1.1.7 infections. The mathematical derivations are described in the supplementary materials in the sections A.3 and A.4. Finally, we show the projected number of confirmed infections in the future under the continuous model. We initialize the model on 01 January 2021 with the estimated number of B.1.1.7 and non-B.1.1.7 confirmed infections on that day. We assume a reproductive number for the non-B.1.1.7 infections as estimated on the national level for 01 January-17 January 2021. Further, we assume that the expected generation time is 4.8 days and the fitness advantage is the estimated fc for the region and dataset of interest (Table 1).

Viruses. 2021 Jul; 13(7): 1216.

Published online 2021 Jun 23. doi: 10.3390/v13071216

PMCID: PMC8310006

PMID: 34201862

Fitness Estimation for Viral Variants in the Context of Cellular Coinfection

Huisheng Zhu,1,† Brent E. Allman,2,† and Katia Koelle1,3,*

Amber M. Smith, Academic Editor and Ruian Ke, Academic Editor 

https://www.ncbi.nlm.nih.gov/pmc/articles/PMC8310006/

QIN: This is within-host competition based fitness model, which is different from SARS-CoV-2 genomic surveillance data!

yeast growth fitness analysis, essentiality

applejack:0.yeast.growth.fitness hqin$ ll *Rmd

-rw-r--r--  1 hqin  staff   3.9K Nov 22 22:24 analyze_growth_fitness.Rmd

Analyze growth fitness

H Qin

11/16/2017

Based on 20131028 code and its 20151012 corrections of comma in gene names. For deletion fitness data, download from http://www-deletion.stanford.edu/YDPM/YDPM_index.html

“In the Deletion Project, strains from the deletion collection were monitored under 9 different media conditions selected for the study of mitochondrial function. 5791 heterozygous diploid and 4706 homozygous diploid deletion strains were monitored in parallel using molecular barcodes on fermentable (YPD, YPDGE) and non-fermentable substrates (YPG, YPE, YPL). The YDPM database contains both the raw data and growth rates calculated for each strain in each media condition. Strains can be searched by ORF or Gene name to access growth measurements and data plots for each strain.”

What are the genotypes in the 4 yeast fitness files?

rm(list=ls()); 
setwd("~/github/yeast.growth.fitness"); 
datapath = "data"; 
my.files = c("Regression_Tc1_het.txt", 
             "Regression_Tc1_hom.txt", 
             "Regression_Tc2_het.txt", 
             "Regression_Tc2_hom.txt");

tb1het = read.table(paste('data/',my.files[1],sep=''), header=T, sep='\t', fill=T)
tb2hom = read.table(paste('data/',my.files[2],sep=''), header=T, sep='\t', fill=T)
tb3het = read.table(paste('data/',my.files[3],sep=''), header=T, sep='\t', fill=T)
tb4hom = read.table(paste('data/',my.files[4],sep=''), header=T, sep='\t', fill=T)

Essetial genes are expected to have YPD measurement in tb1het and tb3het, but not in tb2hom and tb4hom.

tb1het$YPDFlag1 = ifelse( is.na(tb1het$YPD), 'noYPDMeasurement', 'HasYPDMeasurement' )
table(tb1het$YPDFlag1)

## 
## HasYPDMeasurement  noYPDMeasurement 
##              5744               174

tb3het$YPDFlag3 = ifelse( is.na(tb3het$YPD), 'noYPDMeasurement', 'HasYPDMeasurement' )
table(tb3het$YPDFlag3)

## 
## HasYPDMeasurement  noYPDMeasurement 
##              5766               152

tb2hom$YPDFlag2 = ifelse( is.na(tb2hom$YPD), 'noYPDMeasurement', 'HasYPDMeasurement' )
table(tb2hom$YPDFlag2)

## 
## HasYPDMeasurement  noYPDMeasurement 
##              4659              1259

tb4hom$YPDFlag4 = ifelse( is.na(tb4hom$YPD), 'noYPDMeasurement', 'HasYPDMeasurement' )
table(tb4hom$YPDFlag4)

## 
## HasYPDMeasurement  noYPDMeasurement 
##              4718              1200

tb12 = merge(tb1het, tb2hom, by.x='orf', by.y='orf')
tb12 = tb12[, c(1,2,grep('Flag',names(tb12)))]

tb34 = merge(tb3het,tb4hom, by.x='orf', by.y='orf')
tb34 = tb34[, c(1,2,grep('Flag',names(tb34)))]

#names(tb12) = c("orf",'gene','Anno1','Growth1')
#names(tb34) = c('orf','gene','Anno2','Growth2')
str(tb12)

## 'data.frame':    5918 obs. of  4 variables:
##  $ orf     : Factor w/ 5918 levels "YAL001C","YAL002W",..: 1 2 3 4 5 6 7 8 9 10 ...
##  $ gene.x  : Factor w/ 5915 levels "37165","AAC1",..: 3133 3395 750 3424 2971 803 906 2912 1470 3425 ...
##  $ YPDFlag1: chr  "HasYPDMeasurement" "HasYPDMeasurement" "HasYPDMeasurement" "HasYPDMeasurement" ...
##  $ YPDFlag2: chr  "HasYPDMeasurement" "HasYPDMeasurement" "noYPDMeasurement" "HasYPDMeasurement" ...

str(tb34)

## 'data.frame':    5918 obs. of  4 variables:
##  $ orf     : Factor w/ 5918 levels "YAL001C","YAL002W",..: 1 2 3 4 5 6 7 8 9 10 ...
##  $ gene.x  : Factor w/ 5915 levels "37165","AAC1",..: 3133 3395 750 3424 2971 803 906 2912 1470 3425 ...
##  $ YPDFlag3: chr  "HasYPDMeasurement" "HasYPDMeasurement" "HasYPDMeasurement" "HasYPDMeasurement" ...
##  $ YPDFlag4: chr  "HasYPDMeasurement" "HasYPDMeasurement" "noYPDMeasurement" "HasYPDMeasurement" ...

tb = merge(tb12,tb34, by.x='orf', by.y='orf')
head(tb)

##       orf gene.x.x          YPDFlag1          YPDFlag2 gene.x.y
## 1 YAL001C     TFC3 HasYPDMeasurement HasYPDMeasurement     TFC3
## 2 YAL002W     VPS8 HasYPDMeasurement HasYPDMeasurement     VPS8
## 3 YAL003W     EFB1 HasYPDMeasurement  noYPDMeasurement     EFB1
## 4 YAL004W  YAL004W HasYPDMeasurement HasYPDMeasurement  YAL004W
## 5 YAL005C     ssa1 HasYPDMeasurement HasYPDMeasurement     ssa1
## 6 YAL007C     ERP2 HasYPDMeasurement HasYPDMeasurement     ERP2
##            YPDFlag3          YPDFlag4
## 1 HasYPDMeasurement HasYPDMeasurement
## 2 HasYPDMeasurement HasYPDMeasurement
## 3 HasYPDMeasurement  noYPDMeasurement
## 4 HasYPDMeasurement HasYPDMeasurement
## 5 HasYPDMeasurement HasYPDMeasurement
## 6 HasYPDMeasurement HasYPDMeasurement

Consistency check for sanity.

tb$name.check= ifelse(tb$gene.x.x == tb$gene.x.y, T, F)
table(tb$name.check)

## 
## TRUE 
## 5918

#tb$essentiality = NA; 
tb$essentiality = 'abnormal'; 

tb$essentiality = ifelse( tb$YPDFlag2=='HasYPDMeasurement' & tb$YPDFlag4=='HasYPDMeasurement', 'nonessential',  tb$essentiality ); #2 and 4 are homozygous deletion
table(tb$essentiality)

## 
##     abnormal nonessential 
##         1264         4654

tb$essentiality = ifelse( tb$YPDFlag2=='noYPDMeasurement' & tb$YPDFlag4=='noYPDMeasurement' & tb$YPDFlag1 == 'HasYPDMeasurement' & tb$YPDFlag3 == 'HasYPDMeasurement' , 'essential',  tb$essentiality ); #2 and 4 are homozygous deletion, 1 and 3 are heterozynous

table(tb$essentiality)

## 
##     abnormal    essential nonessential 
##          122         1142         4654

sum( table(tb$essentiality) )

## [1] 5918

write.csv(tb, "data/SummaryRegressionHetHom20171122.csv", row.names = F );

#what was I doing here? 
#tbReal$essenflag = ifelse( tbReal$Growth1=='growth1' & tbReal$Growth2=='growth2', 'nonessential', 
#                       ifelse(tbReal$Growth1=='nogrowth1' & tbReal$Growth2=='nogrowth2', 'essential', 'abnormal') )
#head(tbReal)
#table(tbReal$essenflag) 
# abnormal    essential nonessential 
# 68         1152         4552 
#Good, consistent results. Use these for further analysis. 2013 Oct 29
# 20151012, found two "TRUE" in orfs
#length(tbReal$orf) #5772
#length(unique(tbReal$orf)) #5772

#write.csv(tbReal, "SummaryRegressionHetHom2015Oct12.csv", row.names = F ) #change 20151012
#write.csv(tbReal, "SummaryRegressionHetHom2015Oct12.csv",  quote=F, row.names=F )

reviews of yeast fitness data analysis

Code '20131028.growth.fitness.R' analyze the 4 yeast fitness measurements and define 'essential', 'nonessential', and 'abnormal' genes.  This codes output 'SummaryRegressionHetHom2015Oct12.csv' which is the basis of many other codes. The importance of these analysis warrants a separate github repo for these effort. 

fix a bug in orf parsing for 20131028.growth.fitness.R

new output file is  "SummaryRegressionHetHom2015Oct12.csv".
The bug is caused by comma in 1 or 2 orfs, such as "DUR1,2".

# This bug is caused by write.csv(, quote=F, row.names=F)
# The problem is caused by "DUR1,2" an ORF that contains a comma.

markdown output of 20131028.growth.fitness.R

Conver 20131028.growth.fitness.R to Rmd. Output is here. Retro added on Nov 20, 2017. 

Analyze growth fitness

H Qin

11/16/2017

Based on 20131028 code and its 20151012 corrections of comma in gene names. For deletion fitness data, download from http://www-deletion.stanford.edu/YDPM/YDPM_index.html

“In the Deletion Project, strains from the deletion collection were monitored under 9 different media conditions selected for the study of mitochondrial function. 5791 heterozygous diploid and 4706 homozygous diploid deletion strains were monitored in parallel using molecular barcodes on fermentable (YPD, YPDGE) and non-fermentable substrates (YPG, YPE, YPL). The YDPM database contains both the raw data and growth rates calculated for each strain in each media condition. Strains can be searched by ORF or Gene name to access growth measurements and data plots for each strain.”

What are the genotypes in the 4 yeast fitness files?

rm(list=ls())
datapath = "data"
my.files = c("Regression_Tc1_het.txt", 
             "Regression_Tc1_hom.txt", 
             "Regression_Tc2_het.txt", 
             "Regression_Tc2_hom.txt")

tb1 = read.table(paste('data/',my.files[1],sep=''), header=T, sep='\t', fill=T)
tb2 = read.table(paste('data/',my.files[2],sep=''), header=T, sep='\t', fill=T)
tb3 = read.table(paste('data/',my.files[3],sep=''), header=T, sep='\t', fill=T)
tb4 = read.table(paste('data/',my.files[4],sep=''), header=T, sep='\t', fill=T)

tb1$AnnotationFlag = ifelse( is.na(tb1$YPD), 'questionable1', 'verified1' )
table(tb1$AnnotationFlag)

## 
## questionable1     verified1 
##           174          5744

# questional verified 
# 174         5744 
x = table(tb1$orf); str(x)

##  'table' int [1:5918(1d)] 1 1 1 1 1 1 1 1 1 1 ...
##  - attr(*, "dimnames")=List of 1
##   ..$ : chr [1:5918] "YAL001C" "YAL002W" "YAL003W" "YAL004W" ...

tb3$AnnotationFlag = ifelse( is.na(tb3$YPD), 'questionable2', 'verified2' )
table(tb3$AnnotationFlag)

## 
## questionable2     verified2 
##           152          5766

# quetionable    verified 
# 152        5766 
hist(tb3$YPD)

x = table(tb3$orf); str(x)

##  'table' int [1:5918(1d)] 1 1 1 1 1 1 1 1 1 1 ...
##  - attr(*, "dimnames")=List of 1
##   ..$ : chr [1:5918] "YAL001C" "YAL002W" "YAL003W" "YAL004W" ...

tb2$GrowthFlag = ifelse( is.na(tb2$YPD), 'nogrowth1', 'growth1' )
table(tb2$GrowthFlag)

## 
##   growth1 nogrowth1 
##      4659      1259

# growth nogrowth 
# 4659     1259 
x = table(tb2$orf); str(x)

##  'table' int [1:5918(1d)] 1 1 1 1 1 1 1 1 1 1 ...
##  - attr(*, "dimnames")=List of 1
##   ..$ : chr [1:5918] "YAL001C" "YAL002W" "YAL003W" "YAL004W" ...

tb4$GrowthFlag = ifelse( is.na(tb4$YPD), 'nogrowth2', 'growth2' )
table(tb4$GrowthFlag)

## 
##   growth2 nogrowth2 
##      4718      1200

# growth nogrowth 
# 4718     1200 
x = table(tb4$orf); str(x)

##  'table' int [1:5918(1d)] 1 1 1 1 1 1 1 1 1 1 ...
##  - attr(*, "dimnames")=List of 1
##   ..$ : chr [1:5918] "YAL001C" "YAL002W" "YAL003W" "YAL004W" ...

tb12 = merge(tb1,tb2, by.x='orf', by.y='orf')
tb12 = tb12[, c(1,2,grep('Flag',names(tb12)))]

tb34 = merge(tb3,tb4, by.x='orf', by.y='orf')
tb34 = tb34[, c(1,2,grep('Flag',names(tb34)))]

names(tb12) = c("orf",'gene','Anno1','Growth1')
names(tb34) = c('orf','gene','Anno2','Growth2')
str(tb12)

## 'data.frame':    5918 obs. of  4 variables:
##  $ orf    : Factor w/ 5918 levels "YAL001C","YAL002W",..: 1 2 3 4 5 6 7 8 9 10 ...
##  $ gene   : Factor w/ 5915 levels "37165","AAC1",..: 3133 3395 750 3424 2971 803 906 2912 1470 3425 ...
##  $ Anno1  : chr  "verified1" "verified1" "verified1" "verified1" ...
##  $ Growth1: chr  "growth1" "growth1" "nogrowth1" "growth1" ...

str(tb34)

## 'data.frame':    5918 obs. of  4 variables:
##  $ orf    : Factor w/ 5918 levels "YAL001C","YAL002W",..: 1 2 3 4 5 6 7 8 9 10 ...
##  $ gene   : Factor w/ 5915 levels "37165","AAC1",..: 3133 3395 750 3424 2971 803 906 2912 1470 3425 ...
##  $ Anno2  : chr  "verified2" "verified2" "verified2" "verified2" ...
##  $ Growth2: chr  "growth2" "growth2" "nogrowth2" "growth2" ...

tb = merge(tb12,tb34, by.x='orf', by.y='orf')
head(tb)

##       orf  gene.x     Anno1   Growth1  gene.y     Anno2   Growth2
## 1 YAL001C    TFC3 verified1   growth1    TFC3 verified2   growth2
## 2 YAL002W    VPS8 verified1   growth1    VPS8 verified2   growth2
## 3 YAL003W    EFB1 verified1 nogrowth1    EFB1 verified2 nogrowth2
## 4 YAL004W YAL004W verified1   growth1 YAL004W verified2   growth2
## 5 YAL005C    ssa1 verified1   growth1    ssa1 verified2   growth2
## 6 YAL007C    ERP2 verified1   growth1    ERP2 verified2   growth2

tb$name.check= ifelse(tb$gene.x == tb$gene.y, T, F)
table(tb$name.check)

## 
## TRUE 
## 5918

# tb[ tb$name.check==F, ]
length( grep('question', tb$Anno1) ) #174

## [1] 174

length( grep('question', tb$Anno2) ) #152

## [1] 152

#tbReal = tb[ tb$Anno1=='verified1' & tb$Anno2=='verified2', ]
tbReal = tb[ tb$Anno1=='verified1' | tb$Anno2=='verified2', ]  #At least one Het growth. The other one might be an error?!
length( grep('question', tbReal$Anno1) ) #28

## [1] 28

length( grep('question', tbReal$Anno2) ) #6

## [1] 6

table(tbReal$Growth1, tbReal$Growth2)

##            
##             growth2 nogrowth2
##   growth1      4552         5
##   nogrowth1      63      1152

#           growth2 nogrowth2
# growth1      4552         5
# nogrowth1      63      1152

# Verify some manually on SGD.
tbReal[tbReal$Growth1=='nogrowth1' & tbReal$Growth2=='nogrowth2',  ] [1:10,]

##          orf    gene.x     Anno1   Growth1    gene.y     Anno2   Growth2
## 3    YAL003W      EFB1 verified1 nogrowth1      EFB1 verified2 nogrowth2
## 15   YAL016W      tpd3 verified1 nogrowth1      tpd3 verified2 nogrowth2
## 24   YAL025C     mak16 verified1 nogrowth1     mak16 verified2 nogrowth2
## 31   YAL032C     PRP45 verified1 nogrowth1     PRP45 verified2 nogrowth2
## 32   YAL033W      POP5 verified1 nogrowth1      POP5 verified2 nogrowth2
## 34 YAL034W-A      MTW1 verified1 nogrowth1      MTW1 verified2 nogrowth2
## 35 YAL035C-A YAL035C-A verified1 nogrowth1 YAL035C-A verified2 nogrowth2
## 39   YAL038W     cdc19 verified1 nogrowth1     cdc19 verified2 nogrowth2
## 42   YAL041W     cdc24 verified1 nogrowth1     cdc24 verified2 nogrowth2
## 44   YAL043C      pta1 verified1 nogrowth1      pta1 verified2 nogrowth2
##    name.check
## 3        TRUE
## 15       TRUE
## 24       TRUE
## 31       TRUE
## 32       TRUE
## 34       TRUE
## 35       TRUE
## 39       TRUE
## 42       TRUE
## 44       TRUE

tbReal[tbReal$Growth1=='growth1' & tbReal$Growth2=='nogrowth2',  ] 

##          orf gene.x     Anno1 Growth1 gene.y     Anno2   Growth2
## 956  YDR050C   TPI1 verified1 growth1   TPI1 verified2 nogrowth2
## 1448 YEL002C   WBP1 verified1 growth1   WBP1 verified2 nogrowth2
## 2685 YIL078W   THS1 verified1 growth1   THS1 verified2 nogrowth2
## 4152 YML126C  ERG13 verified1 growth1  ERG13 verified2 nogrowth2
## 4383 YMR218C TRS130 verified1 growth1 TRS130 verified2 nogrowth2
##      name.check
## 956        TRUE
## 1448       TRUE
## 2685       TRUE
## 4152       TRUE
## 4383       TRUE

#what was I doing here? 
tbReal$essenflag = ifelse( tbReal$Growth1=='growth1' & tbReal$Growth2=='growth2', 'nonessential', 
                       ifelse(tbReal$Growth1=='nogrowth1' & tbReal$Growth2=='nogrowth2', 'essential', 'abnormal') )
head(tbReal)

##       orf  gene.x     Anno1   Growth1  gene.y     Anno2   Growth2
## 1 YAL001C    TFC3 verified1   growth1    TFC3 verified2   growth2
## 2 YAL002W    VPS8 verified1   growth1    VPS8 verified2   growth2
## 3 YAL003W    EFB1 verified1 nogrowth1    EFB1 verified2 nogrowth2
## 4 YAL004W YAL004W verified1   growth1 YAL004W verified2   growth2
## 5 YAL005C    ssa1 verified1   growth1    ssa1 verified2   growth2
## 6 YAL007C    ERP2 verified1   growth1    ERP2 verified2   growth2
##   name.check    essenflag
## 1       TRUE nonessential
## 2       TRUE nonessential
## 3       TRUE    essential
## 4       TRUE nonessential
## 5       TRUE nonessential
## 6       TRUE nonessential

table(tbReal$essenflag) 

## 
##     abnormal    essential nonessential 
##           68         1152         4552

# abnormal    essential nonessential 
# 68         1152         4552 
#Good, consistent results. Use these for further analysis. 2013 Oct 29
# 20151012, found two "TRUE" in orfs
length(tbReal$orf) #5772

## [1] 5772

length(unique(tbReal$orf)) #5772

## [1] 5772

#write.csv(tbReal, "SummaryRegressionHetHom2015Oct12.csv", row.names = F ) #change 20151012
#write.csv(tbReal, "SummaryRegressionHetHom2015Oct12.csv",  quote=F, row.names=F )