add functions to compute hamming distance pairs for an MSA

src/msa.cpp

@@ -351,3 +351,38 @@ MSA::writeSequenceWeights(std::string output_file)
     output_stream << sequence_weights.at(i) << std::endl;
   }
 };
+
+void
+MSA::computeHammingDistances(void) {
+  hamming_distances = arma::Col<double>(M, arma::fill::zeros);
+  arma::Mat<int> alignment_T = alignment.t();
+
+  int *i_ptr = nullptr;
+  int *j_ptr = nullptr;
+  int count = 0;
+  double id = 0;
+  for (int i = 0; i < M; i++) {
+    i_ptr = alignment_T.colptr(i);
+    for (int j = i + 1; j < M; j++) {
+      count = 0;
+      j_ptr = alignment_T.colptr(j);
+      for (int n = 0; n < N; n++) {
+        if (*(i_ptr+n) == *(j_ptr+n)) {
+          count++;
+        }
+      }
+      id = (double)count / N;
+      if (id > hamming_distances.at(i)) {
+        hamming_distances.at(i) = id;
+      }
+    }
+  }
+};
+
+void
+MSA::writeHammingDistances(std::string output_file) {
+  std::ofstream output_stream(output_file);
+  for (int i = 0; i < M; i++) {
+    output_stream << hamming_distances.at(i) << std::endl;
+  }
+};

src/msa.hpp

@@ -16,11 +16,14 @@ public:
   int N;                              // number of positions
   int Q;                              // number of amino acids
 
+  arma::Col<double> hamming_distances;
+
   MSA(std::string, bool = true, bool = false, double = 0.8);
   MSA(std::string, std::string, bool = false);
   void printAlignment();
   void writeMatrix(std::string);
   void writeSequenceWeights(std::string);
+  void writeHammingDistances(std::string);
 
 private:
   std::vector<SeqRecord> seq_records;
@@ -30,6 +33,7 @@ private:
   void readSequenceWeights(std::string);
   void makeNumericalMatrix(void);
   void computeSequenceWeights(double);
+  void computeHammingDistances(void);
 };
 
 #endif