representative.cpp : compute representative quality if possible

algo/core/representative.cpp

@@ -63,6 +63,10 @@ float KmerRepresentativeComputer::getCoverage() const{
   return coverage;
 }
 
+string KmerRepresentativeComputer::getQuality() const{
+  return quality;
+}
+
 string KmerRepresentativeComputer::getCoverageInfo() const{
   return coverage_info;
 }
@@ -98,6 +102,8 @@ void KmerRepresentativeComputer::compute() {
   size_t length_longest_run = 0;
   size_t seq_index_longest_run = 1;
   Sequence sequence_longest_run;
+  int sequence_used_for_quality = 0;
+  int window_quality_sum [required.length()] = {};
   size_t k = getSeed().length();
 
   for (size_t seq = 1; seq <= sequences.size() && seq <= seq_index_longest_run + stability_limit ; seq++) {
@@ -106,6 +112,7 @@ void KmerRepresentativeComputer::compute() {
     if (sequence.sequence.size() <= length_longest_run) {
       break;
     }
+    
     size_t pos_required = sequence.sequence.find(required);
     if (pos_required == string::npos && revcomp)
       pos_required = sequence.sequence.find(::revcomp(required));
@@ -113,6 +120,14 @@ void KmerRepresentativeComputer::compute() {
     if (pos_required == string::npos) {
       continue;
     }
+    
+    //sum quality
+    if (sequence.quality.length() > 0) {
+      for (int i = 0; i<required.length(); i++) {
+        window_quality_sum[i] += static_cast<int>(sequence.quality[i+pos_required]);
+      }
+      sequence_used_for_quality++;
+    }
 
     size_t pos_end_required = pos_required + required.length();
 
@@ -159,8 +174,26 @@ void KmerRepresentativeComputer::compute() {
       + string_of_int(pos_longest_run + length_longest_run - 1) + "]"
       + "-#" + string_of_int(seq_index_longest_run);
 
-    representative.quality = "";
-
+    quality = "";
+    
+    if (sequence_used_for_quality > 0) {
+      //init default quality
+      for (int i = 0; i < representative.sequence.length(); i++)
+          quality += "!";
+    
+      //add window quality
+      size_t pos_required = representative.sequence.find(required);
+      if (pos_required == string::npos && revcomp) {
+          size_t pos_required_rev = ::revcomp(representative.sequence).find(required);
+          size_t pos_end_required_rev = pos_required_rev+required.length();
+          for (int i = 0; i<required.length(); i++)
+            quality[pos_end_required_rev-i-1] = static_cast<char>(window_quality_sum[i]/sequence_used_for_quality);
+      }else{
+          for (int i = 0; i<required.length(); i++)
+            quality[pos_required+i] = static_cast<char>(window_quality_sum[i]/sequence_used_for_quality);
+      }
+    }
+    
     coverage = (float) length_longest_run / coverage_reference_length;
 
     coverage_info  = string_of_int(length_longest_run) + " bp"

algo/core/representative.h

@@ -24,6 +24,7 @@ protected:
   string required;
   float coverage_reference_length;
   float coverage;
+  string quality;
   string coverage_info;
 public:
   RepresentativeComputer(list<Sequence> &r);
@@ -105,6 +106,7 @@ public:
   // Getters, setters
   string getSeed() const;
   float getCoverage() const;
+  string getQuality() const;
   string getCoverageInfo() const;
 
   /**