draw alignment lines

A step toward getting the alignment view... This draws the alignment
lines. And even gets it right occasionally.

Open issues are: reverse complimented sequences get shifted
by the wrong sign. There's an off by one error (ticket:37), and of course
there's a need for a user interface to pick which alignment to use & show.

diff --git a/alg/glseqbrowser.cpp b/alg/glseqbrowser.cpp
index 147f36c04a84bea08009a37f6358da26c8d425ef..ea5dc2e3f4f6475b000ca032241ac84088da09ae 100644 (file)
--- a/alg/glseqbrowser.cpp
+++ b/alg/glseqbrowser.cpp
@@ -514,8 +514,10 @@ void GlSeqBrowser::draw_segments() const
       // save the multipart name for our segment
       glPushName(path_index); glPushName(key.first); glPushName(key.second);
       glBegin(GL_LINES);
-      glVertex3f(s.start.x, s.start.y, -1);
-      glVertex3f(s.end.x, s.end.y, -1);
+      float seq_start_x = track_container[path_index].x();
+      float seq_end_x = track_container[path_index+1].x();
+      glVertex3f(s.start.x + seq_start_x, s.start.y, -1);
+      glVertex3f(s.end.x   + seq_end_x  , s.end.y, -1);
       glEnd();
       // clear the names
       glPopName(); glPopName(); glPopName();

diff --git a/alg/sequence.cpp b/alg/sequence.cpp
index 26215a4d428f15e15657847013f7a3f8ffa0a2b0..d7bf864bb984fc8abe712ed11d582a62d37b8896 100644 (file)
--- a/alg/sequence.cpp
+++ b/alg/sequence.cpp
@@ -83,6 +83,11 @@ Sequence &Sequence::operator=(const std::string& s)
   return *this;
 }
 
+char Sequence::operator[](int index) const
+{
+  return sequence[index];
+}
+
 ostream& operator<<(ostream& out, const Sequence& seq)
 {
   out << "Sequence(" << seq.get_seq() << ")";

diff --git a/alg/sequence.hpp b/alg/sequence.hpp
index eaf65e7dcf2133e89a05ba56cf79c873377f1853..bd2b864df124fd0ce72db5fbfc6f00adf9188017 100644 (file)
--- a/alg/sequence.hpp
+++ b/alg/sequence.hpp
@@ -72,6 +72,7 @@ class Sequence
     //! assignment to constant sequences
     Sequence &operator=(const Sequence&);
     Sequence &operator=(const std::string &);
+    char operator[](int) const;
     friend std::ostream& operator<<(std::ostream& out, const Sequence& seq);
 
     //! set sequence to a (sub)string containing nothing but AGCTN

diff --git a/qui/MussaAlignedWindow.cpp b/qui/MussaAlignedWindow.cpp
index 46e2daad911797a97cb4cf30ea80a54bc0e0e416..e54f099cef78a941739ab46d354ac687c066fda4 100644 (file)
--- a/qui/MussaAlignedWindow.cpp
+++ b/qui/MussaAlignedWindow.cpp
@@ -1,10 +1,11 @@
 #include <list>
+#include <vector>
 #include <sstream>
 
 #include <QVBoxLayout>
 
 #include "qui/MussaAlignedWindow.hpp"
-
+#include "alg/sequence.hpp"
 
 using namespace std;
 
@@ -20,6 +21,7 @@ MussaAlignedWindow::MussaAlignedWindow(Mussa& m,
   setSelectedPaths(m, sel_paths);
   setAlignment(0);
   browser.zoomToSequence();
+  computeMatchLines();
 
   QVBoxLayout *layout = new QVBoxLayout;
   layout->addWidget(&browser);
@@ -68,3 +70,132 @@ void MussaAlignedWindow::setAlignment(size_t alignment_index)
   }
 }
 
+void MussaAlignedWindow::computeMatchLines()
+{
+  const vector<Sequence>& raw_sequence = analysis.sequences();
+  vector<int> aligned_path;
+  int i2, i3;
+  int x_start, x_end;
+  int window_length, win_i;
+  int rc_1 = 0; 
+  int rc_2 = 0;
+  bool rc_color;
+  vector<bool> rc_list;
+  bool full_match;
+  vector<bool> matched;
+  int align_counter;
+
+  align_counter = 0;
+  for(vector<ExtendedConservedPath >::iterator pathz_i=selected_paths.begin(); 
+      pathz_i != selected_paths.end(); 
+      ++pathz_i)
+  {
+    if (true /* show_aligns[align_counter] */)
+    {
+      ExtendedConservedPath& a_path = *pathz_i;
+      window_length = a_path.window_size;
+      // determine which parts of the path are RC relative to first species
+      rc_list = a_path.reverseComplimented();
+      
+      // loop over each bp in the conserved region for all sequences
+      for(win_i = 0; win_i < window_length; win_i++)
+      {
+        aligned_path.clear();
+        // determine which exact base pairs match between the sequences
+        full_match = true;
+        for(i2 = 0; i2 < a_path.size()-1; i2++)
+        {
+          // assume not rc as most likely, adjust below
+          rc_1 = 0;
+          rc_2 = 0;
+          // no matter the case, any RC node needs adjustments
+          if (a_path[i2] < 0)
+            rc_1 = window_length-1;
+          if (a_path[i2+1] < 0)
+            rc_2 = window_length-1;        
+           
+          x_start = (abs(a_path[i2]-rc_1+win_i));
+          x_end =   (abs(a_path[i2+1]-rc_2+win_i));
+          
+          // RC case handling
+          // ugh, and xor...only want rc coloring if just one of the nodes is rc
+          // if both nodes are rc, then they are 'normal' relative to each other
+          if ( ( rc_list[i2] || rc_list[i2+1] ) &&
+              !(rc_list[i2] && rc_list[i2+1] ) )
+          { //the hideous rc matching logic - not complex, but annoying
+            if ( !( ( (raw_sequence[i2][x_start] == 'A') &&
+                    (raw_sequence[i2+1][x_end] == 'T') ) ||
+                  ( (raw_sequence[i2][x_start] == 'T') &&
+                    (raw_sequence[i2+1][x_end] == 'A') ) ||
+                  ( (raw_sequence[i2][x_start] == 'G') &&
+                    (raw_sequence[i2+1][x_end] == 'C') ) ||
+                  ( (raw_sequence[i2][x_start] == 'C') &&
+                    (raw_sequence[i2+1][x_end] == 'G') ) ) )
+              full_match = false;
+          }
+          else
+          {
+            if (!( (raw_sequence[i2][x_start] == raw_sequence[i2+1][x_end]) &&
+                  (raw_sequence[i2][x_start] != 'N') &&
+                  (raw_sequence[i2+1][x_end] != 'N') ) )
+              full_match = false;
+          }
+        }
+        
+        // draw for matches stretching across all sequences
+        if (full_match)
+        {
+          //fl_line_style(FL_SOLID, 1, NULL);
+          
+          // now can draw the line for each bp in this window that matches
+          // grrr, need to ask if anyone cares if I switch the seq 
+          // top-bot order...
+          i3 = 0;
+          //y_loc = y_min + 5;
+          for(i2 = 0; i2 < a_path.size()-1; i2++)
+          {
+            // assume not rc as most likely, adjust below
+            rc_1 = 0;
+            rc_2 = 0;
+            // no matter the case, any RC node needs adjustments
+            if (a_path[i2] < 0)
+            {
+              rc_1 = window_length;        
+            }
+            if (a_path[i2] < 0)
+            {
+              rc_2 = window_length;        
+            }
+            
+            // set offset based on current alignment for which bp to show
+            //offset1 = seq_align_offsets[i2];
+            //offset2 = seq_align_offsets[i2+1];
+            
+            if ( ( rc_list[i2] || rc_list[i2+1] ) &&
+                !(rc_list[i2] && rc_list[i2+1] ) ) {
+              ;
+              //fl_color(FL_BLUE);
+            } else {
+              ;
+              //fl_color(FL_RED);
+            }
+            
+            // maybe shouldn't recalc these, but store values from first loop
+            x_start = (abs((int) (a_path[i2]-rc_1+win_i)));
+            x_end =   (abs((int) (a_path[i2+1]-rc_2+win_i)));
+            aligned_path.push_back(x_start);
+            // if we're on the last time through the loop, save x_end too
+            if (i2 == a_path.size()-2) {
+              aligned_path.push_back(x_end);
+            }
+          }
+        }
+        if (aligned_path.size() > 0) {
+          browser.link(aligned_path, rc_list,1);
+        }
+      }
+    }
+    align_counter++;
+  }
+}
+

diff --git a/qui/MussaAlignedWindow.hpp b/qui/MussaAlignedWindow.hpp
index b828ea587f361a3b8f88dc049f13aaeafff7fa5d..10f7d168ac28dd7dce1270846efc34af31d09a63 100644 (file)
--- a/qui/MussaAlignedWindow.hpp
+++ b/qui/MussaAlignedWindow.hpp
@@ -21,6 +21,7 @@ public slots:
 
 protected:
   void setSelectedPaths(Mussa &m, const std::set<int>& sel_paths);
+  void computeMatchLines();
   
   Mussa& analysis;
   //std::vector<Sequence> sequences;