alg/glseqbrowser.cpp

   1 #include "alg/glseqbrowser.hpp"
   2 #include "mussa_exceptions.hpp"
   3
   4 #include <iostream>
   5 #include <sstream>
   6 #include <stdexcept>
   7
   8 using namespace std;
   9
  10 GlSeqBrowser::GlSeqBrowser()
  11   : border_width(25),
  12     cur_ortho(400.0, 0.0, 600.0, 0.0),
  13     viewport_size(600, 400),
  14     viewport_center((cur_ortho.right-cur_ortho.left)/2+cur_ortho.left),
  15     zoom_level(2),
  16     color_mapper(),
  17     track_container()
  18 {
  19 }
  20
  21 GlSeqBrowser::GlSeqBrowser(const GlSeqBrowser& gt)
  22   : border_width(gt.border_width),
  23     cur_ortho(gt.cur_ortho),
  24     viewport_size(gt.viewport_size),
  25     viewport_center(gt.viewport_center),
  26     zoom_level(gt.zoom_level),
  27     color_mapper(gt.color_mapper),
  28     track_container(gt.track_container)
  29 {
  30 }
  31
  32 void GlSeqBrowser::initializeGL()
  33 {
  34   glEnable(GL_DEPTH_TEST);
  35   glClearColor(1.0, 1.0, 1.0, 0.0);
  36   glBlendFunc(GL_SRC_ALPHA, GL_ONE_MINUS_SRC_ALPHA);
  37   glShadeModel(GL_FLAT);
  38 }
  39
  40 void GlSeqBrowser::resizeGL(int width, int height)
  41 {
  42   viewport_size.x = width;
  43   viewport_size.y = height;
  44   glViewport(0, 0, (GLsizei)width, (GLsizei)height);
  45   update_viewport(viewport_center, zoom_level);
  46 }
  47
  48 void GlSeqBrowser::paintGL() const
  49 {
  50   glClear(GL_COLOR_BUFFER_BIT|GL_DEPTH_BUFFER_BIT);
  51
  52   glPushMatrix();
  53   glMatrixMode(GL_PROJECTION);
  54   glLoadIdentity();
  55   glOrtho(cur_ortho.left, cur_ortho.right,
  56           cur_ortho.bottom, cur_ortho.top,
  57           -50.0, 50);
  58
  59   draw();
  60
  61   glPopMatrix();
  62   glFlush();
  63 }
  64
  65 void GlSeqBrowser::processSelection(GLuint hits, GLuint buffer[], GLuint bufsize, const rect<float>& r)
  66 {
  67   GLuint *ptr;
  68   GLuint names;
  69   GLuint consumed_names = 0;
  70   float z1;
  71   float z2;
  72   GLuint objtype;
  73   GLuint objid;
  74   GLuint path_index = 0;
  75   GLuint pair_key_0 = 0;
  76   GLuint pair_key_1 = 0;
  77   TrackRegion track;
  78
  79   selected_paths.clear();
  80   selected_tracks.clear();
  81
  82   ptr = (GLuint *) buffer;
  83   if (hits > 0)
  84     selectedMode = true;
  85   for (GLuint i=0; i < hits; ++i)
  86   {
  87     if ((i + 5) > bufsize) {
  88       std::clog << "*** selection overflow***" << std::endl;
  89     } else {
  90       consumed_names = 0;
  91       names = *ptr++;
  92       z1 = ((float)*ptr++)/0x7fffffff;
  93       z2 = ((float)*ptr++)/0x7fffffff;
  94       objtype = *ptr++; ++consumed_names;
  95       switch (objtype) {
  96         case MussaSegment:
  97           path_index = *ptr++; ++consumed_names;
  98           pair_key_0 = *ptr++; ++consumed_names;
  99           pair_key_1 = *ptr++; ++consumed_names;
 100           if (path_index < path_segments.size()) {
 101             segment_key k(pair_key_0, pair_key_1);
 102             pair_segment_map::iterator psm_i;
 103             psm_i = path_segments[path_index].find(k);
 104             if (psm_i != path_segments[path_index].end()) {
 105               Segment &seg = psm_i->second;
 106               selected_paths.insert(seg.path_ids.begin(), seg.path_ids.end());
 107             }
 108             // else something else is wrong
 109           } else {
 110             // something wasn't right
 111             clog << "invalid path_index " << path_index
 112                  << " should have been [0,"<<path_segments.size()
 113                  << ") " << endl;
 114           }
 115           break;
 116         case MussaTrack:
 117           objid = *ptr++; ++consumed_names;
 118
 119           int left = track_container[objid].leftbase(r.left);
 120           int right = track_container[objid].rightbase(r.right);
 121           // the static_cast should be ok, since basepairs line up on
 122           // integral values
 123           //TrackRegion track(objid, left, right);
 124           track.set(objid, left, right);
 125           selected_tracks.push_back(track);
 126         break;
 127         default:
 128           cout << "unknown type " << objtype << " ";
 129           for(; consumed_names < names; ++consumed_names) {
 130             cout << consumed_names << "," << *ptr++ << " ";
 131           }
 132           cout << endl;
 133           break;
 134       }
 135     }
 136   }
 137 }
 138
 139 void GlSeqBrowser::selectRegion(int top, int left, int bottom, int right)
 140 {
 141   GLfloat x_scale = cur_ortho.width()/((float)viewport_size.x);
 142   GLfloat y_scale = cur_ortho.height()/((float)viewport_size.y);
 143   GLfloat x_left = cur_ortho.left + (left*x_scale);
 144   GLfloat x_right = cur_ortho.left + (right * x_scale);
 145
 146   if (top > bottom) {
 147     // woah, someone gave us a rectangle with the origin in the lower left
 148     int temp = top;
 149     bottom = top;
 150     top = temp;
 151   }
 152   // swap the orientation of canvas coordinates
 153   GLfloat y_top = cur_ortho.top-(bottom*y_scale);
 154   GLfloat y_bottom = cur_ortho.top - top * y_scale;
 155   selectedRegion = rect<float>(y_top, x_left, y_bottom, x_right);
 156
 157   // hopefully this will make a buffer big enough to receive
 158   // everything being selected
 159   //const size_t pathz_count = mussaAnalysis->paths().refined_pathz.size();
 160   //const GLuint select_buf_size = 1 + 5 * (pathz_count + sequences.size());
 161   const GLuint select_buf_size = 500000;
 162   GLuint selectBuf[select_buf_size];
 163   glSelectBuffer(select_buf_size, selectBuf);
 164   GLint hits;
 165
 166   (void)glRenderMode(GL_SELECT);
 167   glPushMatrix();
 168   glMatrixMode(GL_PROJECTION);
 169   glLoadIdentity();
 170   glOrtho(x_left, x_right, y_top, y_bottom, -50.0, 50.0);
 171   glMatrixMode(GL_MODELVIEW);
 172   glLoadIdentity();
 173
 174   draw();
 175
 176   glFlush();
 177   glPopMatrix();
 178   hits = glRenderMode(GL_RENDER);
 179   processSelection(hits, selectBuf, select_buf_size, selectedRegion);
 180 }
 181
 182 float GlSeqBrowser::border() const
 183 {
 184   return border_width;
 185 }
 186
 187 float GlSeqBrowser::left() const
 188 {
 189   float left;
 190   if (track_container.size() == 0)
 191   {
 192     return cur_ortho.left;
 193   } else {
 194     vector<GlSequence>::const_iterator track_i = track_container.begin();
 195     left = track_i->x();
 196     for( ; track_i != track_container.end(); ++track_i)
 197     {
 198       if (track_i->x() < left) {
 199         left = track_i->x();
 200       }
 201     }
 202     return left-border_width;
 203   }
 204 }
 205
 206 float GlSeqBrowser::right() const
 207 {
 208   float right;
 209   if (track_container.size() == 0) {
 210     return cur_ortho.right;
 211   } else {
 212     vector<GlSequence>::const_iterator track_i = track_container.begin();
 213     right = track_i->right();
 214     for( ; track_i != track_container.end(); ++track_i) {
 215       if (track_i->right() > right)
 216         right = track_i->right();
 217     }
 218     return right+border_width;
 219   }
 220 }
 221
 222 void GlSeqBrowser::setViewportCenter(float x)
 223 {
 224   update_viewport(x, zoom_level);
 225   viewport_center = x;
 226 }
 227
 228 float GlSeqBrowser::viewportLeft() const
 229 {
 230   return cur_ortho.left;
 231 }
 232
 233 float GlSeqBrowser::viewportCenter() const
 234 {
 235   return viewport_center;
 236 }
 237
 238 float GlSeqBrowser::viewportRight() const
 239 {
 240   return cur_ortho.right;
 241 }
 242
 243 float GlSeqBrowser::viewportHeight() const
 244 {
 245   return cur_ortho.top - cur_ortho.bottom;
 246 }
 247
 248 float GlSeqBrowser::viewportWidth() const
 249 {
 250   return cur_ortho.right - cur_ortho.left;
 251 }
 252
 253 double GlSeqBrowser::zoomOut()
 254 {
 255
 256   if (right() - left() > 0) {
 257     cur_ortho.left = left();
 258     cur_ortho.right = right();
 259     zoom_level =  (right() - left()) / (double)viewport_size.x;
 260     return zoom_level;
 261   } else {
 262     // made up number representing 50 bp / pixel
 263     return 50.0;
 264   }
 265 }
 266
 267 double GlSeqBrowser::zoomToSequence()
 268 {
 269   // (experimentally determined zoom level)
 270   const double friendly_zoom = 0.10;
 271   setZoom(friendly_zoom);
 272   return friendly_zoom;
 273 }
 274
 275 void GlSeqBrowser::setZoom(double new_zoom)
 276 {
 277   update_viewport(viewport_center, new_zoom);
 278   zoom_level = new_zoom;
 279 }
 280
 281 double GlSeqBrowser::zoom() const
 282 {
 283   return zoom_level;
 284 }
 285
 286 void GlSeqBrowser::setColorMapper(AnnotationColors& cm)
 287 {
 288   color_mapper = cm;
 289 }
 290
 291 AnnotationColors& GlSeqBrowser::colorMapper()
 292 {
 293   return color_mapper;
 294 }
 295
 296 void GlSeqBrowser::clear()
 297 {
 298   clear_links();
 299   path_segments.clear();
 300   track_container.clear();
 301 }
 302
 303 void GlSeqBrowser::push_sequence(const Sequence &s)
 304 {
 305   GlSequence gs(s, color_mapper);
 306   push_sequence(gs);
 307 }
 308
 309 void GlSeqBrowser::push_sequence(GlSequence &gs)
 310 {
 311   clear_links();
 312   track_container.push_back(gs);
 313   update_layout();
 314   if (track_container.size() > 1)
 315     path_segments.push_back(pair_segment_map());
 316 }
 317
 318 const std::vector<GlSequence>& GlSeqBrowser::sequences() const
 319 {
 320   return track_container;
 321 }
 322
 323 void GlSeqBrowser::clear_links()
 324 {
 325   path_segments.clear();
 326   for (int i = track_container.size()-1; i > 0; --i)
 327   {
 328     path_segments.push_back(pair_segment_map());
 329   }
 330   pathid = 0;
 331 }
 332
 333 void
 334 GlSeqBrowser::link(const vector<int>& path, const vector<bool>& rc, int )
 335 {
 336   if (path.size() < 2) {
 337     // should i throw an error instead?
 338     return;
 339   }
 340   if (path.size() != track_container.size() ) {
 341     stringstream msg;
 342     msg << "Path size [" << path.size() << "] and track size ["
 343         << track_container.size() << "] don't match" << endl;
 344     throw mussa_error(msg.str());
 345   }
 346   if (path.size() != rc.size()) {
 347     throw runtime_error("path and reverse compliment must be the same length");
 348   }
 349   vector<int>::const_iterator path_i = path.begin();
 350   vector<bool>::const_iterator rc_i = rc.begin();
 351   int track_i = 0;
 352   int prev_x = *path_i; ++path_i;
 353   bool prev_rc = *rc_i; ++rc_i;
 354   while (path_i != path.end() and rc_i != rc.end())
 355   {
 356     segment_key p(prev_x, *path_i);
 357     pair_segment_map::iterator found_segment = path_segments[track_i].find(p);
 358     if (found_segment == path_segments[track_i].end()) {
 359       // not already found
 360       float y1 = track_container[track_i].y();
 361             y1 -= track_container[track_i].height()/2;
 362       float y2 = track_container[track_i+1].y();
 363             y2 += track_container[track_i+1].height()/2;
 364
 365       bool rcFlag = (prev_rc or *rc_i) and !(prev_rc and *rc_i);
 366       Segment s(prev_x, y1, *path_i, y2, rcFlag);
 367       s.path_ids.insert(pathid);
 368       path_segments[track_i][p] = s;
 369     } else {
 370       //found
 371       found_segment->second.path_ids.insert(pathid);
 372     }
 373     prev_x = *path_i;
 374     prev_rc = *rc_i;
 375     ++track_i;
 376     ++path_i;
 377     ++rc_i;
 378   }
 379   // pathid is reset by push_sequence
 380   ++pathid;
 381 }
 382
 383 const set<int>& GlSeqBrowser::selectedPaths() const
 384 {
 385   return selected_paths;
 386 }
 387
 388 //! copy sequence from selected track using formating function
 389 void GlSeqBrowser::copySelectedTracks(std::string& copy_buffer,
 390                                       format_track formatter)
 391 {
 392   copy_buffer.clear();
 393
 394   for(selected_track_iterator track_i = selected_tracks.begin();
 395       track_i != selected_tracks.end();
 396       ++track_i)
 397   {
 398     int track_index = track_i->track_id;
 399     if (track_index >= track_container.size()) {
 400       // should this be an exception instead?
 401       clog << "track " << track_index << " > " << track_container.size()
 402            << endl;
 403     } else {
 404       // we should be safe
 405       const Sequence& seq = track_container[track_index].sequence();
 406       copy_buffer += formatter(seq, track_i->left, track_i->right);
 407     }
 408   }
 409 }
 410
 411 //! copy sequence from selected tracks as plain sequences
 412 void GlSeqBrowser::copySelectedTracksAsString(std::string& copy_buffer)
 413 {
 414   struct AsString {
 415     static string formatter(const Sequence& seq, int left, int right)
 416     {
 417       stringstream s;
 418       s << seq.subseq(left, right-left+1) << std::endl;
 419       return s.str();
 420     }
 421   };
 422
 423   copySelectedTracks(copy_buffer, AsString::formatter);
 424 }
 425
 426 //! copy sequence from selected tracks as FASTA sequences
 427 void GlSeqBrowser::copySelectedTracksAsFasta(std::string& copy_buffer)
 428 {
 429   struct AsFasta {
 430     static string formatter(const Sequence& seq, int left, int right)
 431     {
 432       stringstream s;
 433       s << ">" << seq.get_header() << std::endl
 434         << seq.subseq(left, right-left+1) << std::endl;
 435       return s.str();
 436     }
 437   };
 438   copySelectedTracks(copy_buffer, AsFasta::formatter);
 439 }
 440
 441 void GlSeqBrowser::centerOnPath(const vector<int>& paths)
 442 {
 443   if (paths.size() != track_container.size()) {
 444     throw mussa_error("Path length didn't match the number of sequences");
 445   }
 446
 447   for(size_t track_i = 0; track_i != track_container.size(); ++track_i)
 448   {
 449     // -15 = shift more to the left
 450     track_container[track_i].setX((viewport_center-15) - paths[track_i]);
 451   }
 452 }
 453
 454 void GlSeqBrowser::update_viewport(float center, double new_zoom)
 455 {
 456   // limit how close we can get
 457   if (new_zoom < 0.01) {
 458     new_zoom = 0.01;
 459   }
 460   double new_width = (new_zoom * (float)viewport_size.x);
 461   cur_ortho.left = center-new_width/2.0;
 462   cur_ortho.right = center+new_width/2.0;
 463 }
 464
 465 void GlSeqBrowser::update_layout()
 466 {
 467   typedef std::vector<GlSequence>::iterator glseq_itor_type;
 468   float available_height = (float)cur_ortho.top - 2 * (float)border_width;
 469   float max_base_pairs = 0;
 470   size_t track_count = track_container.size();
 471
 472   if (track_count > 1) {
 473     // we have several sequences
 474     float track_spacing = available_height / (track_count-1);
 475     float y = available_height + (float)border_width;
 476     for(glseq_itor_type seq_i = track_container.begin();
 477         seq_i != track_container.end();
 478         ++seq_i, y-=track_spacing)
 479     {
 480       seq_i->setX(0);
 481       seq_i->setY(y);
 482       if (seq_i->length() > max_base_pairs)
 483         max_base_pairs = seq_i->length();
 484     }
 485   } else if (track_count == 1) {
 486     // center the single track
 487     glseq_itor_type seq_i = track_container.begin();
 488     seq_i->setX(0);
 489     seq_i->setY(viewport_size.x /2);
 490     max_base_pairs = seq_i->length();
 491   } else {
 492     // nothing to do as we're empty
 493     return;
 494   }
 495   cur_ortho.right = max_base_pairs + border_width;
 496   cur_ortho.left = -border_width;
 497   cur_ortho.top = viewport_size.x;
 498   cur_ortho.bottom = 0;
 499   viewport_center = (cur_ortho.width()/2) + cur_ortho.left;
 500   zoomOut();
 501 }
 502
 503 void GlSeqBrowser::draw() const
 504 {
 505   glMatrixMode(GL_MODELVIEW);
 506   glInitNames();
 507   glPushName(MussaSegment);
 508   draw_segments();
 509   glLoadName(MussaTrack);
 510   draw_tracks();
 511   glPopName();
 512   // a selection shouldn't have a glName associated with it
 513   draw_selection();
 514 }
 515
 516 void GlSeqBrowser::draw_selection() const
 517 {
 518   // draw selection box
 519   glEnable(GL_BLEND);
 520   glDepthMask(GL_FALSE);
 521   if (selectedMode) {
 522     glColor4f(0.6, 0.6, 0.6, 0.9);
 523     glRectf(selectedRegion.left, selectedRegion.top,
 524             selectedRegion.right, selectedRegion.bottom);
 525   }
 526   glDepthMask(GL_TRUE);
 527   glDisable(GL_BLEND);
 528 }
 529
 530 void GlSeqBrowser::draw_tracks() const
 531 {
 532   for(size_t track_i = 0; track_i != track_container.size(); ++track_i)
 533   {
 534     glPushName(track_i);
 535     track_container[track_i].draw(cur_ortho.left, cur_ortho.right);
 536     glPopName();
 537   }
 538 }
 539
 540 void GlSeqBrowser::draw_segments() const
 541 {
 542   glLineWidth(1);
 543   // each vector contains path_segment_maps of all the connections
 544   // between this track and the next
 545   path_segment_map_vector::const_iterator psmv_i;
 546   for(psmv_i = path_segments.begin();
 547       psmv_i != path_segments.end();
 548       ++psmv_i)
 549   {
 550     path_segment_map_vector::difference_type path_index;
 551     path_index = psmv_i - path_segments.begin();
 552     // these maps contain the pair index (used so we dont keep drawing the
 553     // same segment) and the actual segment structure.
 554     pair_segment_map::const_iterator psm_i;
 555     for(psm_i = psmv_i->begin();
 556         psm_i != psmv_i->end();
 557         ++psm_i)
 558     {
 559       // grab the index into our segment map
 560       const segment_key& key = psm_i->first;
 561       // the second element of our map pair is a segment
 562       const Segment &s = psm_i->second;
 563       // need to do something so we can detect our selection
 564       vector<int> selected;
 565       set_intersection(selected_paths.begin(), selected_paths.end(),
 566                        s.path_ids.begin(), s.path_ids.end(),
 567                        back_inserter(selected));
 568
 569       if (not s.reversed) {
 570         if (selected_paths.size() == 0 or selected.size() > 0) {
 571           glColor3f(1.0, 0.0, 0.0);
 572         } else {
 573           glColor3f(1.0, 0.8, 0.8);
 574         }
 575       } else {
 576         if (selected_paths.size() == 0 or selected.size() > 0) {
 577           glColor3f(0.0, 0.0, 1.0);
 578         } else {
 579           glColor3f(0.8, 0.8, 1.0);
 580         }
 581         /*
 582         if (selected_paths.size() == 0 or selected.size() > 0) {
 583           glColor3f(0.0, 0.0, 1.0);
 584         } else {
 585           glColor3f(0.8, 0.8, 1.0);
 586         }
 587         */
 588       }
 589       // save the multipart name for our segment
 590       glPushName(path_index); glPushName(key.first); glPushName(key.second);
 591       glBegin(GL_LINES);
 592       float seq_start_x = track_container[path_index].x();
 593       float seq_end_x = track_container[path_index+1].x();
 594       glVertex3f(s.start.x + seq_start_x, s.start.y, -1);
 595       glVertex3f(s.end.x   + seq_end_x  , s.end.y, -1);
 596       glEnd();
 597       // clear the names
 598       glPopName(); glPopName(); glPopName();
 599     }
 600   }
 601 }