monclonal single-source seeding

Author: melkebir

CMakeLists.txt

@@ -84,6 +84,19 @@ set( pmh_pr_src
   src/migrationgraph.h
 )
 
+set( ms_hdr
+  src/utils.h
+  src/frequencymatrix.h
+)
+
+set( ms_src
+  src/ms.cpp
+  src/utils.cpp
+  src/frequencymatrix.cpp
+  src/utils.h
+  src/frequencymatrix.h
+)
+
 set( pmh_cti_hdr
   src/utils.h
   src/nonbinaryclonetree.h
@@ -259,6 +272,9 @@ add_executable( pmh_cti ${pmh_cti_src} ${pmh_cti_hdr} )
 target_include_directories( pmh_cti PUBLIC "${LIBLEMON_ROOT}/include" "src" ${Boost_INCLUDE_DIRS} ${GUROBI_INCLUDE_DIR} )
 target_link_libraries( pmh_cti ${CommonLibs} ${GUROBI_LIBRARIES} )
 
+add_executable( ms ${ms_src} ${ms_hdr} )
+target_link_libraries( ms ${CommonLibs} )
+
 add_executable( generatemigrationtrees ${generatemigrationtrees_src} ${generatemigrationtrees_hdr} )
 target_link_libraries( generatemigrationtrees ${CommonLibs} )
 

src/frequencymatrix.cpp

@@ -18,6 +18,30 @@ FrequencyMatrix::FrequencyMatrix()
 {
 }
 
+bool FrequencyMatrix::mS() const
+{
+  bool res = true;
+  for (int s = 0; s < _m; ++s)
+  {
+    for (int t = 0; t < s; ++t)
+    {
+      for (int i = 0; i < _n; ++i)
+      {
+        if (isSurelySubclonal(s, i) && isSurelySubclonal(t, i))
+        {
+          std::cerr << "Mutation '" << _indexToCharacter[i] << "'"
+            << " is surely subclonal in anatomical sites '"
+            << _indexToSample[s] <<  "' and '" << _indexToSample[t] << "'"
+            << std::endl;
+          res = false;
+        }
+      }
+    }
+  }
+  
+  return res;
+}
+
 std::ostream& operator<<(std::ostream& out, const FrequencyMatrix& F)
 {
   out << F._m << " #samples" << std::endl;

src/frequencymatrix.h

@@ -120,6 +120,26 @@ class FrequencyMatrix
     return Sigma;
   }
 
+  /// Return whether a mutation is surely subclonal
+  ///
+  /// @param s Sample index
+  /// @param i Character index
+  bool isSurelySubclonal(int s, int i) const
+  {
+    for (int j = 0; j < _n; ++j)
+    {
+      if (i == j) continue;
+      
+      if (_f[s][i].second < _f[s][j].first && _f[s][i].first > 0)
+        return true;
+    }
+    
+    return false;
+  }
+
+  /// Return whether (F^-, F+^) satisfy a necessary condition for mS
+  bool mS() const;
+  
 private:
   typedef std::pair<double, double> DoublePair;
   typedef std::vector<DoublePair> DoublePairVector;

src/ms.cpp

@@ -0,0 +1,48 @@
+/*
+ * ms.cpp
+ *
+ *  Created on: 3-jul-2017
+ *      Author: M. El-Kebir
+ */
+
+#include "utils.h"
+#include "frequencymatrix.h"
+#include <fstream>
+#include <lemon/arg_parser.h>
+
+int main(int argc, char** argv)
+{
+  lemon::ArgParser ap(argc, argv);
+  ap.other("frequencies", "Frequencies");
+  ap.parse();
+  
+  if (ap.files().size() != 1)
+  {
+    std::cerr << "Error: <frequencies> must be specified" << std::endl;
+    return 1;
+  }
+  
+  std::string filenameFrequencies = ap.files()[0];
+  
+  std::ifstream inFrequencies(filenameFrequencies.c_str());
+  if (!inFrequencies.good())
+  {
+    std::cerr << "Could not open '" << filenameFrequencies << "' for reading" << std::endl;
+    return 1;
+  }
+  
+  FrequencyMatrix F;
+  inFrequencies >> F;
+  inFrequencies.close();
+  
+  if (F.mS())
+  {
+    std::cerr << "Provided frequency matrix is potentially generated by a clone tree with a monoclonal single-source seeding pattern" << std::endl;
+  }
+  else
+  {
+    std::cerr << "Provided frequency matrix is NOT generated by a clone tree with a monoclonal single-source seeding pattern" << std::endl;
+  }
+  
+  return 0;
+}