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