Rabies-Mutation/code/2b_simulate_perfect_clockrate.R at main · RowanDurrant/Rabies-Mutation · GitHub

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
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
library(seqinr)
library(ape)

k = 24
perTimeRates = c(0.1/(12000*26.3), 0.2/(12000*26.3), 0.5/(12000*26.3), 1/(12000*26.3),
                 2/(12000*26.3),0.05/(12000*26.3), 0.15/(12000*26.3), 0.25/(12000*26.3),
                 0.75/(12000*26.3), 1.5/(12000*26.3), 3/(12000*26.3))
source("code/mutations_function_time.R")
set.seed(758)
#noCases = 150
percentCases = c(0.01, 0.05, 0.1, 0.2)
SI_meanlog <- 2.96
SI_sdlog <-  0.82
branches = c(71,83,124,156,166,179,191,269,271)

for(perTimeRate in perTimeRates){

  for(g in branches){
    load(paste0("input/EF_simulations/seed",k,".Rdata"))
    outMatrix = as.data.frame(out$caseMat)
    outMatrix = outMatrix[order(outMatrix$caseID),]
    initialInfections = outMatrix$caseID[outMatrix$parentID == 0]


    previousGen = g
    descendents = c(previousGen)

    offspring = 1
    while(offspring > 0){
      previousGen = outMatrix$caseID[outMatrix$parentID %in% previousGen]
      descendents = append(descendents, previousGen)
      offspring = length(previousGen)
    }

    outMatrix = outMatrix[outMatrix$caseID %in% descendents,]

    timelim = c(0, max(outMatrix$infD))
    for(i in initialInfections){
      outMatrix$sequence[outMatrix$caseID == i] = paste(rep("a",12000), collapse = "")
    }

    allCases = outMatrix$caseID[order(outMatrix$caseID)]

    for(j in allCases){
      if(is.na(outMatrix$sequence[outMatrix$caseID == j]) == T){
        time = outMatrix$infD[outMatrix$caseID == j] -
          outMatrix$infD[outMatrix$caseID == outMatrix$parentID[outMatrix$caseID ==j]]
        outMatrix$sequence[outMatrix$caseID == j] = viralSeq(outMatrix$sequence[outMatrix$caseID == outMatrix$parentID[outMatrix$caseID ==j]],
                                                             perTimeRate, time)
        print(j)

      }
      if(outMatrix$parentID[outMatrix$caseID == j] == 0){
        outMatrix$timeDiff[outMatrix$caseID == j] = rlnorm(1, meanlog = SI_meanlog, sdlog = SI_sdlog)
      }
      else{
        outMatrix$timeDiff[outMatrix$caseID == j] = outMatrix$infD[outMatrix$caseID == j] - outMatrix$infD[outMatrix$caseID == outMatrix$parentID[outMatrix$caseID ==j]]
      }
    }

    tips = outMatrix
    for(noCases in percentCases*nrow(outMatrix)){
      sampledtips = tips[tips$caseID %in% sample(tips$caseID, noCases),]
      write.csv(sampledtips, paste0("output/simulation/simsampledtips_perfectClockRate_",k,"_",perTimeRate,"_",g,"_",noCases,".csv"))
     # write.fasta(as.list(sampledtips$sequence), paste(sampledtips$caseID, sampledtips$infD, sep = "_"),
    #              file.out = paste0("output/simulation/simsampledtips_perfectClockRate_",k,"_",perTimeRate,"_",g,"_",noCases,".fasta"),
     #             open = "w", nbchar = 60, as.string = F)

     # all.fasta=read.dna(paste0("output/simulation/simsampledtips_perfectClockRate_",k,"_",perTimeRate,"_",noCases,".fasta"), format = 'fasta')
    #  alnDist <- dist.dna(all.fasta, model = "raw", as.matrix = TRUE, pairwise.deletion = T)
   #   write.csv(alnDist, paste0("output/simulation/simsampledsnpdists_perfectClockRate_",k,"_",perTimeRate,"_",noCases,".csv"))

    }

    to = as.character(outMatrix$caseID)
    from = as.character(outMatrix$parentID)
    TimeDiff = outMatrix$timeDiff
    allCases= as.data.frame(cbind(from,to,TimeDiff))
    #write.csv(allCases, paste0("output/simulation/simFullCases_perfectClockRate_",k,"_",perTimeRate,"_",g,".csv"))

  }

}