Working on a bug fix...

diff --git a/htswanalysis/src/GetReadsInSnps/getRegionConsensus.cpp b/htswanalysis/src/GetReadsInSnps/getRegionConsensus.cpp
index 1f0ae1bc712fd332357c818b42cc2ac6eac955b0..4cbd8781ce2b5fe5a677a233d90c18e9b288fe35 100755 (executable)
--- a/htswanalysis/src/GetReadsInSnps/getRegionConsensus.cpp
+++ b/htswanalysis/src/GetReadsInSnps/getRegionConsensus.cpp
@@ -169,6 +169,8 @@ class Nuc {
     }
 
     char consensus() {
+      unsigned int N = size();
+      if(N == 0) { return ' '; }
       unsigned int max = 0; unsigned int max_idx = 0;
       for(unsigned int i = 0; i < 4; i++) { if(n[i] > max) { max = n[i]; max_idx = i; } }
 
@@ -195,8 +197,7 @@ class Nuc {
         }
       }
  
-      unsigned int N = size();
-      if(N == 0) { return ' '; } else if(N < 10) { return tolower(c); } else { return c; }
+      if(N < 10) { return tolower(c); } else { return c; }
     }
 };
 
@@ -222,6 +223,8 @@ class Window : public Loci {
       this->name = name;
       this->length = length; 
       this->sequence = "";
+ 
+      if(length > 10000) { cerr << "ERROR: window of size " << length << endl; exit(1); }
       seq.resize(length);
 
       this->reads = 0;
@@ -240,6 +243,7 @@ class Window : public Loci {
     }
 
     Window& operator=(const Window& r) { 
+     if (this == &r)  return *this;
       Loci::operator=(r);
       this->name = r.name;
       this->length = r.length; 
@@ -266,6 +270,7 @@ class Window : public Loci {
       this->reads++;
       for(unsigned int i = 0; i < r.length(); i++) {
         int seq_idx = offset + i;
+        if(seq_idx > 0 && (unsigned)seq_idx >= seq.size()) { cerr << "Error: offset: " << offset << " seq_idx " << seq_idx << " size: " << seq.size() << " seqlen: " << sequence.length() << endl; }
         if(seq_idx < 0 || (seq_idx >= 0 && (unsigned)seq_idx > this->length) ) { continue; }
         seq[offset + i].add_nuc(r[i]);
       }
@@ -306,13 +311,17 @@ class Window : public Loci {
         unsigned int max_len = sequence.length() - offset;
         unsigned int len = (line_len > max_len)?max_len:line_len;
         for(unsigned int i = offset; i < offset+len; i++) { 
+          if(i >= seq.size()) { 
+            cerr << "Error: offset of " << i << " exceeds read data: " << seq.size() << " in string of " << sequence.length() << endl; 
+            exit(1); 
+          }
           char con = toupper(seq[i].consensus()); 
           // weak consensus if lowercase.
           bool weak_con = seq[i].consensus() != con;
           if(con == ' ' || weak_con || toupper(con) == toupper(sequence[i])) { 
-            output += sequence[i]; 
+            output.append(1,sequence[i]); 
           } else { 
-            output += con; 
+            output.append(1,con); 
             char buff[128];
             sprintf(buff,"%d:%c>%c",i,sequence[i],con);
             string var = buff;
@@ -716,7 +725,7 @@ void count_read_in_features(Windows& windows, Reads& data) {
     //stop if we have run out of features.
     if(wind_it == windows.end()) { break; }
 
-    if(i->pos + i->length > wind_it->pos && i->pos < (wind_it->pos + wind_it->length)) {
+    if(i->pos + i->length() > wind_it->pos && i->pos < (wind_it->pos + wind_it->length)) {
       wind_it->add_read(*i);
     }
   }
@@ -739,7 +748,7 @@ void retrieveSequenceData(ChromList chrom_filenames, Windows& peaks) {
             offset = chrom_file.gcount();
           }
           unsigned int begin = i->pos - 1;
-          unsigned int end   = i->pos+i->length-2;
+          unsigned int end   = i->pos+i->length;
  
           unsigned int begin_pos = offset + (int)begin/50 + begin;      
           unsigned int end_pos = offset + (int)end/50 + end;