Simulations

Author: melkebir

src/simulation/celltree.cpp

@@ -6,7 +6,6 @@
  */
 
 #include "celltree.h"
-#include <fstream>
 
 CellTree::CellTree(double K,
                    double migrationRate,
@@ -140,6 +139,7 @@ void CellTree::migrate()
         {
           assert(t == currentNrAnatomicalSites);
           Node v_t = _G.addNode();
+          _anatomicalSiteMap[v_t] = t;
           _indexToVertexG.push_back(v_t);
           _G.addArc(_indexToVertexG[s], v_t);
 
@@ -152,7 +152,7 @@ void CellTree::migrate()
   }
 }
 
-void CellTree::simulate()
+bool CellTree::simulate()
 {
   std::uniform_real_distribution<> unif(0, 1);
 
@@ -239,7 +239,13 @@ void CellTree::simulate()
     ++_generation;
   }
 
-  constructCloneTree();
+  if (_nrExtantCells > 0)
+  {
+    constructCloneTree();
+    return true;
+  }
+  
+  return false;
 }
 
 MigrationGraph CellTree::getMigrationGraph() const
@@ -357,7 +363,10 @@ void CellTree::constructCloneTree()
   for (ArcIt a(newT); a != lemon::INVALID; ++a)
   {
     int j = arcToCharacter[a];
-    arcToMutationSet[a] = toOriginalColumns[j];
+    if (j != -1)
+    {
+      arcToMutationSet[a] = toOriginalColumns[j];
+    }
   }
 
   char buf[1024];
@@ -371,13 +380,21 @@ void CellTree::constructCloneTree()
     if (v != newRoot)
     {
       label[v] = "";
-      for (int j : toOriginalColumns[arcToCharacter[InArcIt(newT, v)]])
+      if (arcToCharacter[InArcIt(newT, v)] != -1)
       {
-        if (label[v] != "")
-          label[v] += ";";
-        snprintf(buf, 1024, "%d", j);
+        for (int j : toOriginalColumns[arcToCharacter[InArcIt(newT, v)]])
+        {
+          if (label[v] != "")
+            label[v] += ";";
+          snprintf(buf, 1024, "%d", j);
+          label[v] += buf;
+        }
       }
     }
+    else
+    {
+      label[v] = "GL";
+    }
     if (OutArcIt(newT, v) == lemon::INVALID)
     {
       IntSet mutations;
@@ -401,6 +418,14 @@ void CellTree::constructCloneTree()
         snprintf(buf, 1024, "M%d", s);
         anatomicalSiteStr[v] = buf;
       }
+    }
+  }
+  
+  for (NodeIt v(newT); v != lemon::INVALID; ++v)
+  {
+    if (OutArcIt(newT, v) == lemon::INVALID)
+    {
+      label[v] = label[newT.source(InArcIt(newT, v))];
       label[v] += "_" + anatomicalSiteStr[v];
       invMap[label[v]] = v;
     }

src/simulation/celltree.h

@@ -30,8 +30,8 @@ class CellTree
   /// Destructor
   ~CellTree();
 
-  /// Simulate a cell tree
-  void simulate();
+  /// Simulate a cell tree. Return true if simulation was successful and false otherwise.
+  bool simulate();
 
   /// Return migration graph
   MigrationGraph getMigrationGraph() const;

src/simulation/simulationmain.cpp

@@ -8,6 +8,7 @@
 #include "utils.h"
 #include "simulation/celltree.h"
 #include <lemon/arg_parser.h>
+#include <fstream>
 
 int main(int argc, char** argv)
 {
@@ -16,24 +17,60 @@ int main(int argc, char** argv)
   double mutFreqThreshold = 0.05;
   double migrationRate = 0.000001;
   double detectionProb = 7e-8;
+  std::string filenameColorMap;
+  std::string outputDirectory = ".";
 
   lemon::ArgParser ap(argc, argv);
   ap.refOption("s", "Random number generator seed (default: 5)", seed, false)
+    .refOption("c", "Color map file", filenameColorMap, true)
     .refOption("K", "Carrying capacity (default: 10e4)", K)
     .refOption("f", "Mutation frequency threshold (default: 0.05)", mutFreqThreshold)
     .refOption("m", "Migration rate (default: 0.000001)", migrationRate)
-    .refOption("d", "Detection probability (default: 7e-8)", detectionProb);
+    .refOption("d", "Detection probability (default: 7e-8)", detectionProb)
+    .refOption("o", "Output directory (default: '.')", outputDirectory);
+  ap.parse();
+  
+  char buf[1024];
+  
+  // Read color map
+  StringToIntMap colorMap;
+  std::ifstream inColorMap(filenameColorMap.c_str());
+  if (!inColorMap.good())
+  {
+    std::cerr << "Could not open '" << filenameColorMap << "' for reading" << std::endl;
+    return 1;
+  }
+  if (!BaseTree::readColorMap(inColorMap, colorMap))
+  {
+    return 1;
+  }
 
   g_rng = std::mt19937(seed);
-//  for (int i = 0; i < 100; ++i)
+  CellTree cellTree(K, migrationRate, detectionProb, mutFreqThreshold);
+  if (cellTree.simulate())
   {
-    CellTree cellTree(K, migrationRate, detectionProb, mutFreqThreshold);
-    cellTree.simulate();
-    
     const NonBinaryCloneTree& T = cellTree.getCloneTree();
-    StringToIntMap colorMap = T.generateColorMap();
+    MigrationGraph G = cellTree.getMigrationGraph();
+    
+    snprintf(buf, 1024, "%s/T_seed%d.dot", outputDirectory.c_str(), seed);
+    std::ofstream outT_dot(buf);
+    T.writeDOT(outT_dot, colorMap);
+    outT_dot.close();
+    
+    snprintf(buf, 1024, "%s/T_seed%d.tree", outputDirectory.c_str(), seed);
+    std::ofstream outT_tree(buf);
+    T.write(outT_tree);
+    outT_tree.close();
+    
+    snprintf(buf, 1024, "%s/G_seed%d.dot", outputDirectory.c_str(), seed);
+    std::ofstream outG_dot(buf);
+    G.writeDOT(outG_dot, colorMap);
+    outG_dot.close();
 
-    T.writeDOT(std::cout, colorMap);
+    snprintf(buf, 1024, "%s/G_seed%d.tree", outputDirectory.c_str(), seed);
+    std::ofstream outG_tree(buf);
+    G.write(outG_tree);
+    outG_tree.close();
   }
 
   return 0;