lineageSoftware/simMemoirRandomGuess.R at master · agranado/lineageSoftware · 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
84
85
86
87
88
89
90


#chat with Mark and grace:
#look for optimal edit rate

#April 2nd
#Simulate a tree and then calcualte the distance to a randomized version of itself
#returns a vector with the distances of many comparisions (many randomizations)
#Update March 30th, function works well
simMemoirRandomGuess<-function(nGen,mu,alpha,barcodeLength,nRand){
  #load necessary libraries and functions
  library(phangorn)
  source("simulation2.R")


  os=system("cat ../os.txt",intern = T)
  if(os=="mac"){ #local Alejandro's laptop
    pathName="/Users/alejandrog/MEGA/Caltech/trees/simulation/"
    pathName2="/Users/alejandrog/MEGA/Caltech/trees/simulation"
  }else if(os=="linux"){ #AWS server or any other linux machine (path is for AWS)
    pathName="/home/ubuntu/alejandrog/Caltech/lineage/"
    pathName2="/home/ubuntu/alejandrog/Caltech/lineage"

  }


  #Estimate ML trees:
  #ML trees does not significantly increase the performance of lineage reconstruction
  #It does, actually, performa way worse than normal upgma/distance methods.
  #Moreover, it is computationally expensive, especially when performing parameter analysis.
  #if you want to estimate ML edit the following line:
  estimML=FALSE


  #clear the variable (since it behaves as global)
  if(exists("firstCell")){
    rm(firstCell)
  }
  #create intialize the tree using 1 as the ID
  #Node is a function from data.tree
  firstCell <- Node$new("1")

  #number of letters (this will change for simulation)
  #A goal is to understand how the lenght of the pattern affects the ability to reconstruc lineages

  #intialitze the barcode using "u->unchanged
  firstCell$barcode <-paste(rep("u",barcodeLength),collapse="")

  #all variables of the data.tree structure are global
  #any pointer to a node, points to the real tree.


  for (g in 1:nGen){
    #this function simulates one generation of the tree
    divideCellRecursive2(firstCell,mu,alpha)
  }

  #prints only the barcodes for all leaves
 # print(firstCell,"barcode")
  print("Tree simulation completed")
  #save to file as newick tree
  #save the length of branches plus the ID (which so far is a number)
  newickTree<-ToNewick(firstCell)
  #Generate unique ID for writing file to disc (we'll erase it later)
  fileID = toString(runif(1))
  firstCellFile = paste(pathName,"trees/firstCell",fileID,".nwk",sep="")
  write(newickTree,file=firstCellFile)
  #load the tree from the file as a tree structure.
  trueTree<-read.tree(file=firstCellFile)


  #file is now deleted
  print("True tree read")
  # plot(trueTree,main=paste("True tree ",sep=""))
  randDist = array()
  for (r in 1:nRand){
    randomTrue<- randomTrueTree(trueTree)
    randDist[r]=RF.dist(randomTrue,trueTree)
  }


  system(paste("rm ",firstCellFile,sep=""))

  return(randDist)
}