emvolz

require(ape)
require(treedater)
load( 'Untitled.RData' ) 

s <- 2153922 # sequence length 
bounds  <- data.frame( lower = dates, upper = dates + 1 )
(td0 <- dater( unroot(t), dates, s, minblen = 7/365, estimateSampleTimes = bounds, ncpu = 6))
pb0 <- parboot.treedater( td0 , overrideTempConstraint=FALSE, ncpu = 6)

#~ > td0

emvolz

compute.paf <- function(THETA){
  tfgy <- dm( THETA, x0 = X0 , t0 = 1978, t1 = 2014, integrationMethod='lsoda' )
  fs <- tfgy[[2]]
  times <- tfgy[[1]]
  pafs <- sapply( fs, function(f){
    paf <- rowSums( f)[1:3] ; paf <- paf /sum(paf )
  })
  X <- cbind( times, t(pafs) )
  X[ order(times), ]
}

emvolz

#' This version separetely simulates the model producing a trajectory object `tfgy`
#' The likelihood is then computed using this trajectory and the tree 
#' An additional likelihood term compares the simulated trajectory to surveillance data giving the proportion of infeced men who are MSM

#' @param tfgy simuation trajectory produced by demographic.process.model
#' @value proportion of infected men in 2010 that are MSM 
tfgy2prevalenceStat <- function(tfgy)
{
  times <- tfgy[[1]]

emvolz

# |A|
# | |
# | |--B-|
# | |    |
# | |C|  |
# | | |  |D|
# | | |  | |
# a b c  d e

emvolz

#' This script simulates genealogies for a "Binary state speciation and extinction model" (BiSSE)
#' Trees are simulated using a structured coalescent framework implemented in the phydynR package
#' In this model, a wild type mutates at a constant rate to the mutant type which has lower transmission fitness. 
#'
#' @author Erik M Volz <erik.volz@gmail.com>
#' @date 14 January 2020 

library( phydynR )
library( ape )
library( ggtree )

emvolz

invisible('
Comparing likelihood functions for estimating TMRCA assuming star phylogeny
') 

library( lubridate) 

# The data: 
d = data.frame( 
 date = ymd( c("2020/01/16", "2020/01/17", "2019/12/30", "2019/12/30", "2019/12/30", "2019/12/30", "2019/12/30", "2020/01/13", "2020/01/08", "2019/12/30", "2019/12/30", "2019/12/30", "2019/12/24", "2019/12/26", "2019/12/30", "2019/12/30", "2019/12/30", "2020/01/17", "2020/01/15", "2020/01/14") )
 , snps = c(2, 0, 1, 2, 2, 0, 2, 0, 0, 0, 2, 0, 3, 0, 0, 1, 0, 1, 1, 3)

emvolz

invisible('
Extinction probability as a function of 
- Reproduction number R
- Dispersion k
- Number imported infections

Assumes negative binomial offspring distribution
') 

# negative binomial probability generating function 

emvolz

<!--
transmission
-->
<matrixeq id="MatrixEquation.2" spec="phydyn.model.MatrixEquation" destination="E" origin="Il" type="birth">beta*Il*S / ( S + E + Il + Ih +R )</matrixeq>
<matrixeq id="MatrixEquation.3" spec="phydyn.model.MatrixEquation" destination="E" origin="Ih" type="birth">beta*Ih*tau * S / (S + E + Il + Ih +R)</matrixeq>
<!--
stage progression
-->
<matrixeq id="MatrixEquation.0" spec="phydyn.model.MatrixEquation" destination="Il" origin="E" type="migration">if ( t &gt; SEIR_START ) then gamma0*E*(1-p_h) else 0.0</matrixeq>
<matrixeq id="MatrixEquation.1" spec="phydyn.model.MatrixEquation" destination="Ih" origin="E" type="migration">if ( t &gt; SEIR_START ) then gamma0*E*p_h else 0.0</matrixeq>

emvolz

library( sarscov2 )
library( ape )

# define some parameters 
## path to the GISAID alignment: 
fn = 'large_alignment_may5_Erik_Brazil.fasta' 
## A simple name for the region:
region='brazil'
## A path to the file containing small genetic distance pairs
distfn = 'large_alignment_may5_Erik_Brazil.dst'

emvolz

library( sarscov2 )
library( ape )

# define some parameters 
## path to the GISAID alignment, distance file, metadata: 
fn = 'may10_filter1.fas' 
distfn = 'dist.txt'
mdfn = 'may10_dataC.csv'

## A simple name for the region: