-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathMain_Run.R
45 lines (36 loc) · 1.28 KB
/
Main_Run.R
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
#Loading libraries
library(ape)
library(TransPhylo)
library(stringr)
library(reticulate)
library(lubridate)
library(readr)
#Load the ptree from the Neewick tree
ptree<-ptreeFromPhylo(read.tree('MUSCLE_beast_adjusted_tree.nwk'),dateLastSample=2018.222)
brLenChecker(ptree)
#set parameters for AutoSummary
sscalep<-1.5
sshape<-2
gscalep<-1.5
gshape<-2
pop<- 139
c<-autoSummary(ptree,pop,0.99,F,1,gshape,gscalep,sshape,scalep)
#Run TransReport
transReport(c,pop)
#Load TransReport output
py_run_file('predictVSKnown.py')
transmission_report <- read_delim("transmission_report.txt",
"\t", escape_double = FALSE, trim_ws = TRUE)
knowninfec<-data.frame(ID=transmission_report$`Actual Infector`)
knowninfec<-na.omit(knowninfec)
preddata<-data.frame(SAC_ID=transmission_report$`SAC ID`,ID=transmission_report$`Actual Infector`,PID=transmission_report$`Pred. Infector`,Dates=transmission_report$`Pred. Infect. Date`)
preddata<-na.omit(preddata)
#calculate correct transmissions
correct<-preddata[(preddata$ID==preddata$PID),4]
correct<-decimal_date(ymd(correct))
lenc<-length(correct)
#calculate incorrect transmissions
incorrect<-preddata[(preddata$ID!=preddata$PID),4]
incorrect<-decimal_date(ymd(incorrect))
#generate plotCTree
plotCTree(c,TRUE,NA,NA,correct,incorrect)