nyc_skygrowth.R

# Note this depends on some convenience functions in the 'sarscov2' package available at 
# https://github.com/emvolz-phylodynamics/sarscov2Rutils
library( sarscov2 )
library( ape ) 
library( lubridate )

# copy over the ref sequence 
file.copy( system.file( package='sarscov2', 'extdata/ref.fasta' ) , '.' ) 

# align 
system( 'mafft --thread 4 --keeplength --add NYUMC.volz.2020jun03.fa ref.fasta > algn3.fasta')

readline() 

d = read.dna( 'algn3.fasta', 'fasta')
md = read.table( 'md.txt', header=TRUE, sep = '\t' , stringsAs=FALSE)

d1 <- d[ rownames(d)!='NC_045512.2', ]

# get collection dates from metadata
coldates <- as.Date( md$Clinical.Collection.Date [match( rownames(d1), md$GISAID.Virus.Name ) ])

# put tip labels in standard format 
rownames(d1) = paste( sep = '|', rownames(d1), decimal_date( coldates ), 'Il' )
write.dna( d1, file = 'algn4.fasta', format='fasta')

# estimate time trees 
tds = make_starting_trees( 'algn4.fasta', ntres = 20 , ncpu = 6)
saveRDS(tds, file = 'tds.rds' )
# NOTE you can also smooth out node ties using treedater::gibbs_jitter. Skygrowth on these trees gives slighly later peak

# estimate Ne 
sg1 = skygrowth1( tds , res = 30 , tau0 = 1e-4)

# plots
library( cowplot )
p1= ggGR_sarscov2skygrowth( sg1 , date_limits=as.Date(c('2020-03-01','2020-04-21' )) )
p0 = ggNe_sarscov2skygrowth( sg1, date_limits=as.Date(c('2020-03-01','2020-04-21' )) )
p = plot_grid( plotlist = list( p0, p1 )  )
ggsave( p, file='skygrowth-newyork-2020-06-04-a0.png' , width = 7 , height=3.75)
ggsave( p, file='skygrowth-newyork-2020-06-04-a0.svg' , width = 7 , height=3.75)

# time of max Ne 
tmax = apply(sg1$GRmat, MAR=2, FUN=function(x) {
	i <- cumsum( x < 0 )
	sg1$taxis[which(i==3)[1]-3] # 3 consecutive days of decrease
})
tmax = unlist( tmax )
date_decimal( quantile( tmax, c(.5, .025, .975 )) ) 

'
> date_decimal( quantile( tmax, c(.5, .025, .975 )) ) 
                      50%                      2.5%                     97.5% 
"2020-03-29 22:01:58 UTC" "2020-03-19 14:09:50 UTC" "2020-04-05 02:45:14 UTC" 

'