Beta binomial

Author: melkebir

CMakeLists.txt

@@ -39,6 +39,7 @@ set( pmh_hdr
 set( simulation_hdr
   src/simulation/cell.h
   src/simulation/celltree.h
+  src/simulation/beta_distribution.hpp
   src/nonbinaryclonetree.h
   src/clonetree.h
   src/binarytree.h

src/simulation/celltree.cpp

@@ -6,6 +6,7 @@
  */
 
 #include "celltree.h"
+#include "beta_distribution.hpp"
 
 CellTree::CellTree(double K,
                    double migrationRate,
@@ -173,7 +174,15 @@ void CellTree::simulateReadCounts()
 
         int coverage = poisson(g_rng);
         // divide 2, heterozygous diploid
-        std::binomial_distribution<> binom(coverage, _freq[s][p][i] / 2);
+        double f = _freq[s][p][i] / 2;
+        int alpha = std::min(1, int(f * 100));
+        int beta = 100 - alpha;
+        
+        sftrabbit::beta_distribution<> beta_dist(alpha, beta);
+        double ff = beta_dist(g_rng);
+        
+//        std::binomial_distribution<> binom(coverage, f);
+        std::binomial_distribution<> binom(coverage, ff);
         _var[s][p][i] = binom(g_rng);
         _ref[s][p][i] = coverage - _var[s][p][i];
       }
@@ -187,16 +196,20 @@ void CellTree::writeReadCounts(std::ostream& out) const
 {
   int nrMutations = getNrMutations();
 
-  out << _nrAnatomicalSites << " #anatomical sites" << std::endl;
+  out << _nrActiveAnatomicalSites << " #anatomical sites" << std::endl;
   out << _nrSamplesPerAnatomicalSite * _nrAnatomicalSites << " #samples" << std::endl;
   out << nrMutations << " #mutations" << std::endl;
 
   out << "#sample_index\tsample_label\tanatomical_site_index\tanatomical_site_label"\
   "\tcharacter_index\tcharacter_label\tref\tvar" << std::endl;
 
   char buf[1024];
+  int activeSiteIndex = 0;
   for (int s = 0; s < _nrAnatomicalSites; ++s)
   {
+    if (!_isActiveAnatomicalSite[s])
+      continue;
+    
     std::string label_s;
     if (s == 0)
     {
@@ -210,18 +223,20 @@ void CellTree::writeReadCounts(std::ostream& out) const
 
     for (int p = 0; p < _nrSamplesPerAnatomicalSite; ++p)
     {
-      int pp = s * _nrSamplesPerAnatomicalSite + p;
+      int pp = activeSiteIndex * _nrSamplesPerAnatomicalSite + p;
 
       snprintf(buf, 1024, "%s_%d", label_s.c_str(), p);
       std::string label_p = buf;
       for (int i = 0; i < nrMutations; ++i)
       {
         out << pp << "\t" << label_p << "\t"
-            << s << "\t" << label_s << "\t"
+            << activeSiteIndex << "\t" << label_s << "\t"
             << i << "\t" << i << "\t"
             << _ref[s][p][i] << "\t" << _var[s][p][i] << std::endl;
       }
     }
+    
+    ++activeSiteIndex;
   }
 }
 

src/simulation/simulationmain.cpp

@@ -138,7 +138,7 @@ int main(int argc, char** argv)
                       mutFreqThreshold,
                       maxNrAnatomicalSites,
                       2,
-                      160,
+                      200,
                       migrationPattern);
     if (cellTree.simulate(true))
     {
@@ -158,7 +158,7 @@ int main(int argc, char** argv)
                         mutFreqThreshold,
                         maxNrAnatomicalSites,
                         2,
-                        160,
+                        200,
                         migrationPattern);
       if (cellTree.simulate(false))
       {
@@ -174,6 +174,7 @@ int main(int argc, char** argv)
           std::cout << "FAILURE" << std::endl;
         }
       }
+      else
       {
         std::cout << "FAILURE" << std::endl;
       }