From: Diane Trout <diane@caltech.edu>
Date: Sat, 18 Mar 2006 08:23:15 +0000 (+0000)
Subject: center on path
X-Git-Url: http://woldlab.caltech.edu/gitweb/?p=mussa.git;a=commitdiff_plain;h=305f9a01baaf833414bfe245ffbfde6dba4c6e19

center on path
This manages to make the various sequences line up along some provided path.
(a vector<int> listing the sequences positions).

This has problems because of the use of the negative flag to indicate
reverse compliment (as it moves the sequence the wrong direction).

Additionally to get this to work I had to get all of the glsequence code
necessary to allow non-zero x values.
---

diff --git a/alg/conserved_path.cpp b/alg/conserved_path.cpp
index 56253c9..af8af7d 100644
--- a/alg/conserved_path.cpp
+++ b/alg/conserved_path.cpp
@@ -58,7 +58,7 @@ std::ostream& operator<<(std::ostream& out, const ConservedPath& path)
   {
     // print the first element
     ConservedPath::const_iterator itor = path.begin();
-    out << ">" << *itor;
+    out << ">" << *itor++;
     // print the remaining elements
     for (;
         itor != path.end();
@@ -186,4 +186,10 @@ ExtendedConservedPath& ExtendedConservedPath::extend(int growth)
   return *this;
 }
 
+std::ostream& operator<<(std::ostream& out, const ExtendedConservedPath& path)
+{
+  out << "<extended_path size=" << path.window_size;
+  out << (ConservedPath&)path;
+  return out;
+}
 
diff --git a/alg/conserved_path.hpp b/alg/conserved_path.hpp
index 17f1bb7..cc2dda1 100644
--- a/alg/conserved_path.hpp
+++ b/alg/conserved_path.hpp
@@ -57,6 +57,9 @@ struct ExtendedConservedPath : public ConservedPath
    */
   ExtendedConservedPath& extend(int growth=1);
 
+  // output some useful information
+  friend std::ostream& operator<<(std::ostream&, const ExtendedConservedPath&);
+
   //! size of extended window
   int window_size;
 };
diff --git a/alg/glseqbrowser.cpp b/alg/glseqbrowser.cpp
index 70637a0..96496a2 100644
--- a/alg/glseqbrowser.cpp
+++ b/alg/glseqbrowser.cpp
@@ -7,11 +7,10 @@
 using namespace std;
 
 GlSeqBrowser::GlSeqBrowser()
-  : border(25),
-    max_ortho(400.0, 0.0, 600.0, 0.0),
-    cur_ortho(max_ortho),
+  : border_width(25),
+    cur_ortho(400.0, 0.0, 600.0, 0.0),
     viewport_size(600, 400),
-    viewport_center(0),
+    viewport_center((cur_ortho.right-cur_ortho.left)/2+cur_ortho.left),
     zoom_level(2),
     color_mapper(),
     track_container()
@@ -19,7 +18,7 @@ GlSeqBrowser::GlSeqBrowser()
 }
 
 GlSeqBrowser::GlSeqBrowser(const GlSeqBrowser& gt)
-  : border(gt.border),
+  : border_width(gt.border_width),
     max_ortho(gt.max_ortho),
     cur_ortho(gt.cur_ortho),
     viewport_size(gt.viewport_size),
@@ -173,14 +172,45 @@ void GlSeqBrowser::selectRegion(int top, int left, int bottom, int right)
   processSelection(hits, selectBuf, select_buf_size);
 }
 
+float GlSeqBrowser::border() const
+{
+  return border_width;
+}
+
 float GlSeqBrowser::left() const
 { 
-  return max_ortho.left; 
+  float left;
+  if (track_container.size() == 0)
+  {
+    return cur_ortho.left;
+  } else {
+    vector<GlSequence>::const_iterator track_i = track_container.begin();    
+    left = track_i->x();
+    for( ; track_i != track_container.end(); ++track_i)
+    {
+      if (track_i->x() < left)
+        left = track_i->x();
+    }
+    return left-border_width;
+  }
 }
 
 float GlSeqBrowser::right() const
 { 
-  return max_ortho.right; 
+  float right;
+  if (track_container.size() == 0)
+  {
+    return cur_ortho.right;
+  } else {
+    vector<GlSequence>::const_iterator track_i = track_container.begin();
+    right = track_i->right();
+    for( ; track_i != track_container.end(); ++track_i)
+    {
+      if (track_i->right() > right)
+        right = track_i->right();
+    }
+    return right+border_width;
+  }
 }
 
 void GlSeqBrowser::setViewportCenter(float x)
@@ -322,6 +352,18 @@ const set<int>& GlSeqBrowser::selectedPaths() const
   return selected_paths;
 }
 
+void GlSeqBrowser::centerOnPath(const vector<int>& paths)
+{
+  if (paths.size() != track_container.size()) {
+    throw mussa_error("Path length didn't match the number of sequences");
+  }
+
+  for(size_t track_i = 0; track_i != track_container.size(); ++track_i)
+  {
+    track_container[track_i].setX(viewport_center - paths[track_i]);
+  }
+}
+
 void GlSeqBrowser::update_viewport(float center, int new_zoom)
 {
   float max_width = max_ortho.width();
@@ -337,14 +379,14 @@ void GlSeqBrowser::update_viewport(float center, int new_zoom)
 void GlSeqBrowser::update_layout()
 {
   typedef std::vector<GlSequence>::iterator glseq_itor_type;
-  float available_height = (float)cur_ortho.top - 2 * (float)border;
+  float available_height = (float)cur_ortho.top - 2 * (float)border_width;
   float max_base_pairs = 0;
   size_t track_count = track_container.size();
 
   if (track_count > 1) {
     // we have several sequences
     float track_spacing = available_height / (track_count-1);
-    float y = available_height + (float)border;
+    float y = available_height + (float)border_width;
     for(glseq_itor_type seq_i = track_container.begin();
         seq_i != track_container.end();
         ++seq_i, y-=track_spacing)
@@ -364,8 +406,8 @@ void GlSeqBrowser::update_layout()
     // nothing to do as we're empty
     return;
   }
-  max_ortho.right = max_base_pairs + border;
-  max_ortho.left = -border;
+  max_ortho.right = max_base_pairs + border_width;
+  max_ortho.left = -border_width;
   max_ortho.top = viewport_size.x;
   max_ortho.bottom = 0;
   cur_ortho = max_ortho;
diff --git a/alg/glseqbrowser.hpp b/alg/glseqbrowser.hpp
index af7efdc..fb9e490 100644
--- a/alg/glseqbrowser.hpp
+++ b/alg/glseqbrowser.hpp
@@ -26,6 +26,9 @@ public:
   //! select a region (using canvas coordinates)
   void GlSeqBrowser::selectRegion(int top, int left, int bottom, int right);
 
+  //! border size
+  float border() const;
+
   //! max world left coordinate
   float left() const;
   //! max world right coordinate
@@ -60,6 +63,8 @@ public:
   void link(const std::vector<int>& path, const std::vector<bool>& isRC, int length);
   //! returns the index of pathids based on order added by link
   const std::set<int>& selectedPaths() const;
+  //! center the provided path in the current viewport
+  void centerOnPath(const std::vector<int>&);
 
   //! Provide a logical name for a type discriminator for our glName stack
   enum FeatureType { MussaTrack, MussaSegment };
@@ -131,7 +136,7 @@ private:
   void draw_selection() const;
 
   //! number of pixels to reserve around the edges of our canvas
-  const int border;
+  const int border_width;
   //! the maximum dimensions that our scene is taking up (world coord)
   rect<float> max_ortho;
   //! the current viewable region (world coord)
diff --git a/alg/glsequence.cpp b/alg/glsequence.cpp
index aaedf18..4a3b989 100644
--- a/alg/glsequence.cpp
+++ b/alg/glsequence.cpp
@@ -60,6 +60,11 @@ GLfloat GlSequence::x() const
   return seq_x;
 }
 
+GLfloat GlSequence::right() const
+{
+  return length()-seq_x;
+}
+
 void GlSequence::setY(GLfloat value)
 {
   seq_y = value;
@@ -82,9 +87,8 @@ GLfloat GlSequence::length() const
 
 Sequence::size_type GlSequence::leftbase(GLfloat left) const
 {
-  assert (seq_x == 0);
-  left = ceil(left);
-  if (left < seq_x)
+  left = ceil(left) - seq_x;
+  if (left < 0)
     return 0;
   else if (left > seq.size() )
     return seq.size();
@@ -94,11 +98,12 @@ Sequence::size_type GlSequence::leftbase(GLfloat left) const
 
 Sequence::size_type GlSequence::rightbase(GLfloat right) const
 {
-  assert (seq_x == 0);
-  right = floor(right);
+  right = floor(right) - seq_x;
   if (right > seq.size())
     return seq.size();
-  else
+  else if ( right < 0) 
+    return 0;
+  else 
     return (Sequence::size_type)right;
 }
 
@@ -115,10 +120,6 @@ Sequence::const_iterator GlSequence::sequence_end() const
 Sequence::const_iterator 
 GlSequence::sequence_begin(GLfloat left, GLfloat right) const
 {
-  // the following code will be wrong when sequences can be slid around
-  // so make sure we break.
-  assert (seq_x == 0);
-
   if ( leftbase(left) > seq.size() or left > right )
     return seq.end();
   else
@@ -128,10 +129,6 @@ GlSequence::sequence_begin(GLfloat left, GLfloat right) const
 Sequence::const_iterator 
 GlSequence::sequence_end(GLfloat left, GLfloat right) const
 {
-  // the following code will be wrong when sequences can be slid around
-  // so make sure we break.
-  assert (seq_x == 0);
-
   if ( rightbase(right) > seq.size() or left > right )
     return seq.end();
   else
diff --git a/alg/glsequence.hpp b/alg/glsequence.hpp
index 3916687..b6ca6eb 100644
--- a/alg/glsequence.hpp
+++ b/alg/glsequence.hpp
@@ -26,6 +26,8 @@ public:
   void setX(GLfloat);
   //! get our starting x (horizontal) coordinate
   GLfloat x() const;
+  //! get our right (horizontal) coordinate (length-x)
+  GLfloat right() const;
   //! set our current y (vertical) position
   void setY(GLfloat);
   //! get our current y (vertical) position
diff --git a/alg/nway_paths.hpp b/alg/nway_paths.hpp
index 23f7698..3e23f15 100644
--- a/alg/nway_paths.hpp
+++ b/alg/nway_paths.hpp
@@ -67,7 +67,11 @@ class NwayPaths
 
     std::list<ConservedPath>::iterator pbegin() { return pathz.begin() ; }
     std::list<ConservedPath>::iterator pend() { return pathz.end() ; }
+    std::list<ConservedPath>::const_iterator pbegin() const { return pathz.begin() ; }
+    std::list<ConservedPath>::const_iterator pend() const { return pathz.end() ; }
     std::list<ExtendedConservedPath>::iterator rpbegin() { return refined_pathz.begin() ; }
+    std::list<ExtendedConservedPath>::const_iterator rpend() const { return refined_pathz.end() ; }
+    std::list<ExtendedConservedPath>::const_iterator rpbegin() const { return refined_pathz.begin() ; }
     std::list<ExtendedConservedPath>::iterator rpend() { return refined_pathz.end() ; }
     // these probably shouldn't be public, but lets start 
     // simple
diff --git a/alg/test/test_glseqbrowser.cpp b/alg/test/test_glseqbrowser.cpp
index ba4fcae..3641d08 100644
--- a/alg/test/test_glseqbrowser.cpp
+++ b/alg/test/test_glseqbrowser.cpp
@@ -45,3 +45,42 @@ BOOST_AUTO_TEST_CASE ( gltracks_connect )
   BOOST_CHECK_EQUAL( gt.path_segments.size(), 0 );
 }
 
+BOOST_AUTO_TEST_CASE( glseqbrowser_center )
+{
+  string s0("AAGGCCTT");
+  string s1("TTGGCCAA");
+  string s2("GATTACAA");
+  Sequence seq0(s0);
+  Sequence seq1(s1);
+  Sequence seq2(s2);
+  AnnotationColors cm;
+  GlSequence glseq0(seq0, cm);
+  GlSequence glseq1(seq1, cm);
+  GlSequence glseq2(seq2, cm);
+
+  GlSeqBrowser gt;
+  gt.push_sequence(glseq0);
+  gt.push_sequence(glseq1);
+  gt.push_sequence(glseq2);
+
+  vector<int> path;
+  path += 0,3,7; 
+
+  BOOST_CHECK_EQUAL( gt.left(), -gt.border() );
+  BOOST_CHECK_EQUAL( gt.right(), s0.size() + gt.border() );
+
+  gt.centerOnPath(path);
+
+  // I didn't bother to compute how the path should shift things
+  // by the magic number 3, i just ran the test and saw how it failed
+  // and hard coded this in.
+  BOOST_CHECK_EQUAL( gt.left(), -3-gt.border() );
+  BOOST_CHECK_EQUAL( gt.right(), 3 + s0.size() + gt.border() );
+
+  // aparently we end up with a different glsequence in the seqbrowser
+  BOOST_CHECK( glseq1.x() != gt.sequences()[1].x() );
+
+  BOOST_CHECK_EQUAL( gt.sequences()[0].x(), gt.viewportCenter()-path[0] );
+  BOOST_CHECK_EQUAL( gt.sequences()[1].x(), gt.viewportCenter()-path[1] );
+  BOOST_CHECK_EQUAL( gt.sequences()[2].x(), gt.viewportCenter()-path[2] );
+}
diff --git a/alg/test/test_glsequence.cpp b/alg/test/test_glsequence.cpp
index 81443a5..9c78530 100644
--- a/alg/test/test_glsequence.cpp
+++ b/alg/test/test_glsequence.cpp
@@ -106,8 +106,53 @@ BOOST_AUTO_TEST_CASE( glsequence_leftright_base )
 
   BOOST_CHECK_EQUAL( glseq.leftbase( -50.0 ), 0 );
   BOOST_CHECK_EQUAL( glseq.leftbase(   0.5 ), 1 );
+  BOOST_CHECK_EQUAL( glseq.leftbase(   5.0 ), 5 );
   BOOST_CHECK_EQUAL( glseq.leftbase( 500.0 ), seq_string.size() );
+  BOOST_CHECK_EQUAL( glseq.rightbase(    0.0 ), 0 );
   BOOST_CHECK_EQUAL( glseq.rightbase( 1000.0 ), seq_string.size() );
   BOOST_CHECK_EQUAL( glseq.rightbase( seq_string.size()-0.5),
                      seq_string.size()-1);
 }
+
+// do our left and right most base computations still work if
+// we move the sequence around?
+BOOST_AUTO_TEST_CASE( glsequence_movex )
+{
+  AnnotationColors cm;
+  std::string seq_string = "AAGGCCTTAAGGCCTT";
+  Sequence seq(seq_string);
+  GlSequence glseq(seq, cm);
+
+  glseq.setX(-5);
+  BOOST_CHECK_EQUAL( glseq.leftbase(-100.0), 0 );
+  BOOST_CHECK_EQUAL( glseq.leftbase(   0.0), 5 );
+  BOOST_CHECK_EQUAL( glseq.leftbase(  -2.0), 3 );
+  BOOST_CHECK_EQUAL( glseq.leftbase( 100.0), seq_string.size() );
+  BOOST_CHECK_EQUAL( glseq.rightbase( 1000.0 ), seq_string.size() );
+  BOOST_CHECK_EQUAL( glseq.rightbase(    8.0 ), 8+5 );
+  BOOST_CHECK_EQUAL( glseq.rightbase(   -7.0 ), 0 );
+}
+
+// Check iterators
+BOOST_AUTO_TEST_CASE( glsequence_check_iterators )
+{
+  AnnotationColors cm;
+  std::string seq_string = "AAGGCCTTAAGGCCTT";
+  Sequence seq(seq_string);
+  GlSequence glseq(seq, cm);
+
+  Sequence::const_iterator seq_begin_i;
+  Sequence::const_iterator seq_end_i;
+
+  BOOST_CHECK(glseq.sequence_begin(5, -5) == seq.end());
+  BOOST_CHECK(glseq.sequence_begin(0, 20) == seq.begin());
+  BOOST_CHECK(glseq.sequence_begin(10,20) == seq.begin()+10);
+
+  BOOST_CHECK(glseq.sequence_end(5, -5) == seq.end());
+  BOOST_CHECK(glseq.sequence_end(0, 20) == seq.end());
+  BOOST_CHECK(glseq.sequence_end(0, 10) == seq.begin()+10);
+
+  glseq.setX(-5);
+  BOOST_CHECK(glseq.sequence_begin(0, 10) == seq.begin()+5);
+  BOOST_CHECK(glseq.sequence_end(0, 15) == seq.end());
+}
diff --git a/qui/MussaAlignedWindow.cpp b/qui/MussaAlignedWindow.cpp
index bd647a5..28e5d1a 100644
--- a/qui/MussaAlignedWindow.cpp
+++ b/qui/MussaAlignedWindow.cpp
@@ -1,21 +1,69 @@
-#include "qui/MussaAlignedWindow.hpp"
+#include <list>
+#include <sstream>
 
 #include <QVBoxLayout>
 
+#include "qui/MussaAlignedWindow.hpp"
+
+
 using namespace std;
 
 MussaAlignedWindow::MussaAlignedWindow(Mussa& m, 
                                        const set<int>& sel_paths, 
                                        QWidget *parent)
   : QWidget(parent),
-    analysis(m),
-    selected_paths(sel_paths)
+    analysis(m)//,
+    //selected_paths(sel_paths)
     
 {
   browser.setSequences(analysis.sequences(), analysis.colorMapper());
+  setSelectedPaths(m, sel_paths);
+  setAlignment(0);
 
   QVBoxLayout *layout = new QVBoxLayout;
   layout->addWidget(&browser);
+  layout->addWidget(&status);
   setLayout(layout);
+
+  ostringstream message;
+  message << "Selected " << selected_paths.size() << " paths";
+  status.showMessage(message.str().c_str(), 5000);
+}
+
+
+void MussaAlignedWindow::setSelectedPaths(Mussa &m, const set<int>& sel_paths)
+{
+  // sets are sorted
+  set<int>::iterator sel_i = sel_paths.begin();
+  list<ExtendedConservedPath>::const_iterator path_i = m.paths().rpbegin();
+  list<ExtendedConservedPath>::const_iterator path_end = m.paths().rpend();
+  size_t path_size = m.paths().refined_pathz.size();
+  size_t pathid=0;
+
+  selected_paths.reserve(sel_paths.size());
+  while (pathid != path_size and sel_i != sel_paths.end())
+  {
+    assert (*sel_i >= 0);
+    size_t sel_pathid = (size_t)(*sel_i);
+    if (pathid == sel_pathid) {
+      cout << "found one " << *path_i << endl;
+      selected_paths.push_back(*path_i);
+      ++pathid;
+      ++path_i;
+      ++sel_i;
+    } else if (pathid < sel_pathid) {
+      ++pathid;
+      ++path_i;
+    } else if (pathid > sel_pathid) {
+      ++sel_i;
+    }
+  }
+}
+
+void MussaAlignedWindow::setAlignment(size_t alignment_index)
+{
+  if (selected_paths.size() > 0) {
+    browser.centerOnPath(selected_paths[alignment_index].track_indexes);
+  }
 }
 
diff --git a/qui/MussaAlignedWindow.hpp b/qui/MussaAlignedWindow.hpp
index a9cbbfc..244a4d5 100644
--- a/qui/MussaAlignedWindow.hpp
+++ b/qui/MussaAlignedWindow.hpp
@@ -3,6 +3,7 @@
 
 #include <set>
 
+#include <QStatusBar>
 #include <QWidget>
 #include "alg/mussa.hpp"
 #include "qui/seqbrowser/SequenceBrowserWidget.hpp"
@@ -11,12 +12,20 @@
 class MussaAlignedWindow : public QWidget
 {
 public:
-  MussaAlignedWindow(Mussa& m, const std::set<int>& sel_paths, QWidget *parent=0);
+  //! construct an aligned window for an analysis and selected paths
+  MussaAlignedWindow(Mussa&, const std::set<int>&, QWidget *parent=0);
+
+public slots:
+  //! use selected_paths[pathid] to set the starting position of our sequence
+  void setAlignment(size_t pathid);
 
 protected:
+  void setSelectedPaths(Mussa &m, const std::set<int>& sel_paths);
+  
   Mussa &analysis;
-  const std::set<int>& selected_paths;
-
+  //const std::set<int>& selected_paths;
+  std::vector<ExtendedConservedPath> selected_paths;
   SequenceBrowserWidget browser;
+  QStatusBar status;
 };
 #endif
diff --git a/qui/seqbrowser/SequenceBrowserWidget.cpp b/qui/seqbrowser/SequenceBrowserWidget.cpp
index 8d4f44e..5571742 100644
--- a/qui/seqbrowser/SequenceBrowserWidget.cpp
+++ b/qui/seqbrowser/SequenceBrowserWidget.cpp
@@ -95,6 +95,11 @@ const std::set<int> SequenceBrowserWidget::selectedPaths() const
   return scrollable_browser.browser().selectedPaths();
 }
 
+void SequenceBrowserWidget::centerOnPath(const vector<int>& paths)
+{
+  scrollable_browser.browser().centerOnPath(paths);
+}
+
 /* This could theoretically be pushed down to some set 
  * of signals and slots connecting SequenceDescriptions and 
  * some signal emitted by the browser's viewportChanged code
diff --git a/qui/seqbrowser/SequenceBrowserWidget.hpp b/qui/seqbrowser/SequenceBrowserWidget.hpp
index 87f5b28..55d55f9 100644
--- a/qui/seqbrowser/SequenceBrowserWidget.hpp
+++ b/qui/seqbrowser/SequenceBrowserWidget.hpp
@@ -39,7 +39,8 @@ public:
   void link(const std::vector<int> &path, const std::vector<bool>& isRC, int length);
   // return set of pathids defined by order of link calls
   const std::set<int> selectedPaths() const;
-  
+  void centerOnPath(const std::vector<int>& paths);
+ 
 public slots:
   //! set the zoom level of our browser
   void setZoom(int);