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         {
 118           objid = *ptr++; ++consumed_names;
 119
 120           int left = track_container[objid].leftbase(r.left);
 121           int right = track_container[objid].rightbase(r.right);
 122           // the static_cast should be ok, since basepairs line up on
 123           // integral values
 124           //TrackRegion track(objid, left, right);
 125           track.set(objid, left, right);
 126           selected_tracks.push_back(track);
 127         }
 128         break;
 129         default:
 130           cout << "unknown type " << objtype << " ";
 131           for(; consumed_names < names; ++consumed_names) {
 132             cout << consumed_names << "," << *ptr++ << " ";
 133           }
 134           cout << endl;
 135           break;
 136       }
 137     }
 138   }
 139 }
 140
 141 void GlSeqBrowser::selectRegion(int top, int left, int bottom, int right)
 142 {
 143   GLfloat x_scale = cur_ortho.width()/((float)viewport_size.x);
 144   GLfloat y_scale = cur_ortho.height()/((float)viewport_size.y);
 145   GLfloat x_left = cur_ortho.left + (left*x_scale);
 146   GLfloat x_right = cur_ortho.left + (right * x_scale);
 147
 148   if (top > bottom) {
 149     // woah, someone gave us a rectangle with the origin in the lower left
 150     int temp = top;
 151     bottom = top;
 152     top = temp;
 153   }
 154   // swap the orientation of canvas coordinates
 155   GLfloat y_top = cur_ortho.top-(bottom*y_scale);
 156   GLfloat y_bottom = cur_ortho.top - top * y_scale;
 157   selectedRegion = rect<float>(y_top, x_left, y_bottom, x_right);
 158
 159   // hopefully this will make a buffer big enough to receive
 160   // everything being selected
 161   //const size_t pathz_count = mussaAnalysis->paths().refined_pathz.size();
 162   //const GLuint select_buf_size = 1 + 5 * (pathz_count + sequences.size());
 163   const GLuint select_buf_size = 500000;
 164   GLuint selectBuf[select_buf_size];
 165   glSelectBuffer(select_buf_size, selectBuf);
 166   GLint hits;
 167
 168   (void)glRenderMode(GL_SELECT);
 169   glPushMatrix();
 170   glMatrixMode(GL_PROJECTION);
 171   glLoadIdentity();
 172   glOrtho(x_left, x_right, y_top, y_bottom, -50.0, 50.0);
 173   glMatrixMode(GL_MODELVIEW);
 174   glLoadIdentity();
 175
 176   draw();
 177
 178   glFlush();
 179   glPopMatrix();
 180   hits = glRenderMode(GL_RENDER);
 181   processSelection(hits, selectBuf, select_buf_size, selectedRegion);
 182 }
 183
 184 float GlSeqBrowser::border() const
 185 {
 186   return border_width;
 187 }
 188
 189 float GlSeqBrowser::left() const
 190 {
 191   float left;
 192   if (track_container.size() == 0)
 193   {
 194     return cur_ortho.left;
 195   } else {
 196     vector<GlSequence>::const_iterator track_i = track_container.begin();
 197     left = track_i->x();
 198     for( ; track_i != track_container.end(); ++track_i)
 199     {
 200       if (track_i->x() < left) {
 201         left = track_i->x();
 202       }
 203     }
 204     return left-border_width;
 205   }
 206 }
 207
 208 float GlSeqBrowser::right() const
 209 {
 210   float right;
 211   if (track_container.size() == 0) {
 212     return cur_ortho.right;
 213   } else {
 214     vector<GlSequence>::const_iterator track_i = track_container.begin();
 215     right = track_i->right();
 216     for( ; track_i != track_container.end(); ++track_i) {
 217       if (track_i->right() > right)
 218         right = track_i->right();
 219     }
 220     return right+border_width;
 221   }
 222 }
 223
 224 void GlSeqBrowser::setViewportCenter(float x)
 225 {
 226   update_viewport(x, zoom_level);
 227   viewport_center = x;
 228 }
 229
 230 float GlSeqBrowser::viewportLeft() const
 231 {
 232   return cur_ortho.left;
 233 }
 234
 235 float GlSeqBrowser::viewportCenter() const
 236 {
 237   return viewport_center;
 238 }
 239
 240 float GlSeqBrowser::viewportRight() const
 241 {
 242   return cur_ortho.right;
 243 }
 244
 245 float GlSeqBrowser::viewportHeight() const
 246 {
 247   return cur_ortho.top - cur_ortho.bottom;
 248 }
 249
 250 float GlSeqBrowser::viewportWidth() const
 251 {
 252   return cur_ortho.right - cur_ortho.left;
 253 }
 254
 255 double GlSeqBrowser::zoomOut()
 256 {
 257
 258   if (right() - left() > 0) {
 259     cur_ortho.left = left();
 260     cur_ortho.right = right();
 261     zoom_level =  (right() - left()) / (double)viewport_size.x;
 262     return zoom_level;
 263   } else {
 264     // made up number representing 50 bp / pixel
 265     return 50.0;
 266   }
 267 }
 268
 269 double GlSeqBrowser::zoomToSequence()
 270 {
 271   // (experimentally determined zoom level)
 272   const double friendly_zoom = 0.10;
 273   setZoom(friendly_zoom);
 274   return friendly_zoom;
 275 }
 276
 277 void GlSeqBrowser::setZoom(double new_zoom)
 278 {
 279   update_viewport(viewport_center, new_zoom);
 280   zoom_level = new_zoom;
 281 }
 282
 283 double GlSeqBrowser::zoom() const
 284 {
 285   return zoom_level;
 286 }
 287
 288 void GlSeqBrowser::setColorMapper(boost::shared_ptr<AnnotationColors> cm)
 289 {
 290   color_mapper = cm;
 291 }
 292
 293 const AnnotationColors& GlSeqBrowser::colorMapper()
 294 {
 295   return *color_mapper;
 296 }
 297
 298 void GlSeqBrowser::clear()
 299 {
 300   clear_selection();
 301   clear_links();
 302   path_segments.clear();
 303   track_container.clear();
 304 }
 305
 306 void GlSeqBrowser::clear_selection()
 307 {
 308   selectedMode = false;
 309   selectedRegion.clear();
 310   selected_paths.clear();
 311   selected_tracks.clear();
 312 }
 313
 314 void GlSeqBrowser::push_sequence(const Sequence& s)
 315 {
 316   boost::shared_ptr<Sequence> seq_copy(new Sequence(s));
 317   GlSequence gs(seq_copy, color_mapper);
 318   push_sequence(gs);
 319 }
 320
 321
 322 void GlSeqBrowser::push_sequence(boost::shared_ptr<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 void GlSeqBrowser::setSelectedPaths(std::vector<int> paths)
 403 {
 404   selected_paths.clear();
 405   for(std::vector<int>::iterator itor = paths.begin();
 406       itor != paths.end();
 407       ++itor)
 408   {
 409     selected_paths.insert(*itor);
 410   }
 411 }
 412
 413 const set<int>& GlSeqBrowser::selectedPaths() const
 414 {
 415   return selected_paths;
 416 }
 417
 418 void GlSeqBrowser::appendSelectedTrack(GLuint track, int start, int stop)
 419 {
 420   selected_tracks.push_back(TrackRegion(track, start, stop));
 421 }
 422
 423 list<TrackRegion> GlSeqBrowser::selectedTracks() const
 424 {
 425   return selected_tracks;
 426 }
 427
 428 //! copy sequence from selected track using formating function
 429 template<class Item>
 430 void GlSeqBrowser::copySelectedTracks(std::list<Item>& result,
 431              Item (*formatter)(const Sequence& s, int left, int right))
 432 {
 433   result.clear();
 434
 435   for(selected_track_iterator track_i = selected_tracks.begin();
 436       track_i != selected_tracks.end();
 437       ++track_i)
 438   {
 439     int track_index = track_i->track_id;
 440     if (track_index >= track_container.size()) {
 441       // should this be an exception instead?
 442       clog << "track " << track_index << " > " << track_container.size()
 443            << endl;
 444     } else {
 445       // we should be safe
 446       const Sequence& seq = track_container[track_index].sequence();
 447       result.push_back(formatter(seq, track_i->left, track_i->right));
 448     }
 449   }
 450 }
 451
 452 //! copy sequence from selected tracks as FASTA sequences
 453 void GlSeqBrowser::copySelectedTracksAsFasta(std::string& copy_buffer)
 454 {
 455   std::list<std::string> result;
 456   struct AsFasta {
 457     static string formatter(const Sequence& seq, int left, int right)
 458     {
 459       stringstream s;
 460       s << ">" << seq.get_fasta_header()
 461         << "|" << "subregion=" << left << "-" << right+1
 462         << std::endl
 463         << seq.subseq(left, right-left+1) << std::endl;
 464       return s.str();
 465     }
 466   };
 467   copySelectedTracks(result, AsFasta::formatter);
 468   // I wish there was some way to use for_each and bind here
 469   for (list<string>::iterator result_i = result.begin();
 470        result_i != result.end();
 471        ++result_i)
 472   {
 473     copy_buffer.append(*result_i);
 474   }
 475 }
 476
 477 //! copy sequence from selected tracks as new sequences
 478 void GlSeqBrowser::copySelectedTracksAsSequences(std::list<Sequence>& result)
 479 {
 480   struct AsSequence {
 481     static Sequence formatter(const Sequence& seq, int left, int right)
 482     {
 483       return seq.subseq(left, right-left+1);
 484     }
 485   };
 486   copySelectedTracks(result, AsSequence::formatter);
 487 }
 488
 489 void GlSeqBrowser::copySelectedTracksAsSeqLocation(
 490     std::list<SequenceLocation>& result)
 491 {
 492   struct AsSeqLocation {
 493     static SequenceLocation formatter(const Sequence& seq, int left, int right)
 494     {
 495       return SequenceLocation(seq, left, right);
 496     }
 497   };
 498   copySelectedTracks(result, AsSeqLocation::formatter);
 499 }
 500
 501 //! copy sequence from selected tracks as plain sequences
 502 void GlSeqBrowser::copySelectedTracksAsString(std::string& copy_buffer)
 503 {
 504   std::list<string> result;
 505   struct AsString {
 506     static string formatter(const Sequence& seq, int left, int right)
 507     {
 508       stringstream s;
 509       s << seq.subseq(left, right-left+1) << std::endl;
 510       return s.str();
 511     }
 512   };
 513
 514   copySelectedTracks(result, AsString::formatter);
 515   // I wish there was some way to use for_each and bind here
 516   for (list<string>::iterator result_i = result.begin();
 517        result_i != result.end();
 518        ++result_i)
 519   {
 520     copy_buffer.append(*result_i);
 521   }
 522
 523 }
 524
 525 void GlSeqBrowser::centerOnPath(const vector<int>& paths)
 526 {
 527   if (paths.size() != track_container.size()) {
 528     throw mussa_error("Path length didn't match the number of sequences");
 529   }
 530
 531   for(size_t track_i = 0; track_i != track_container.size(); ++track_i)
 532   {
 533     // -15 = shift more to the left
 534     track_container[track_i].setX((viewport_center-15) - paths[track_i]);
 535   }
 536 }
 537
 538 void GlSeqBrowser::update_viewport(float center, double new_zoom)
 539 {
 540   // limit how close we can get
 541   if (new_zoom < 0.01) {
 542     new_zoom = 0.01;
 543   }
 544   double new_width = (new_zoom * (float)viewport_size.x);
 545   cur_ortho.left = center-new_width/2.0;
 546   cur_ortho.right = center+new_width/2.0;
 547 }
 548
 549 void GlSeqBrowser::update_layout()
 550 {
 551   typedef std::vector<GlSequence>::iterator glseq_itor_type;
 552   float available_height = (float)cur_ortho.top - 2 * (float)border_width;
 553   float max_base_pairs = 0;
 554   size_t track_count = track_container.size();
 555
 556   if (track_count > 1) {
 557     // we have several sequences
 558     float track_spacing = available_height / (track_count-1);
 559     float y = available_height + (float)border_width;
 560     for(glseq_itor_type seq_i = track_container.begin();
 561         seq_i != track_container.end();
 562         ++seq_i, y-=track_spacing)
 563     {
 564       seq_i->setX(0);
 565       seq_i->setY(y);
 566       if (seq_i->size() > max_base_pairs)
 567         max_base_pairs = seq_i->size();
 568     }
 569   } else if (track_count == 1) {
 570     // center the single track
 571     glseq_itor_type seq_i = track_container.begin();
 572     seq_i->setX(0);
 573     seq_i->setY(viewport_size.x /2);
 574     max_base_pairs = seq_i->size();
 575   } else {
 576     // nothing to do as we're empty
 577     return;
 578   }
 579   cur_ortho.right = max_base_pairs + border_width;
 580   cur_ortho.left = -border_width;
 581   cur_ortho.top = viewport_size.x;
 582   cur_ortho.bottom = 0;
 583   viewport_center = (cur_ortho.width()/2) + cur_ortho.left;
 584   zoomOut();
 585 }
 586
 587 void GlSeqBrowser::draw() const
 588 {
 589   glMatrixMode(GL_MODELVIEW);
 590   glInitNames();
 591   glPushName(MussaSegment);
 592   draw_segments();
 593   glLoadName(MussaTrack);
 594   draw_tracks();
 595   glPopName();
 596   // a selection shouldn't have a glName associated with it
 597   draw_selection();
 598 }
 599
 600 void GlSeqBrowser::draw_selection() const
 601 {
 602   // draw selection box
 603   glEnable(GL_BLEND);
 604   glDepthMask(GL_FALSE);
 605   if (selectedMode) {
 606     glColor4f(0.6, 0.6, 0.6, 0.9);
 607     glRectf(selectedRegion.left, selectedRegion.top,
 608             selectedRegion.right, selectedRegion.bottom);
 609   }
 610   glDepthMask(GL_TRUE);
 611   glDisable(GL_BLEND);
 612 }
 613
 614 void GlSeqBrowser::draw_tracks() const
 615 {
 616   for(size_t track_i = 0; track_i != track_container.size(); ++track_i)
 617   {
 618     glPushName(track_i);
 619     track_container[track_i].draw(cur_ortho.left, cur_ortho.right);
 620     glPopName();
 621   }
 622 }
 623
 624 void GlSeqBrowser::draw_segments() const
 625 {
 626   glLineWidth(1);
 627   glEnable(GL_BLEND);
 628   glDepthMask(GL_FALSE);
 629   // each vector contains path_segment_maps of all the connections
 630   // between this track and the next
 631   path_segment_map_vector::const_iterator psmv_i;
 632   for(psmv_i = path_segments.begin();
 633       psmv_i != path_segments.end();
 634       ++psmv_i)
 635   {
 636     path_segment_map_vector::difference_type path_index;
 637     path_index = psmv_i - path_segments.begin();
 638     // these maps contain the pair index (used so we dont keep drawing the
 639     // same segment) and the actual segment structure.
 640     pair_segment_map::const_iterator psm_i;
 641     for(psm_i = psmv_i->begin();
 642         psm_i != psmv_i->end();
 643         ++psm_i)
 644     {
 645       // grab the index into our segment map
 646       const segment_key& key = psm_i->first;
 647       // the second element of our map pair is a segment
 648       const Segment &s = psm_i->second;
 649       // need to do something so we can detect our selection
 650       vector<int> selected;
 651       set_intersection(selected_paths.begin(), selected_paths.end(),
 652                        s.path_ids.begin(), s.path_ids.end(),
 653                        back_inserter(selected));
 654
 655       if (not s.reversed) {
 656         if (selected_paths.size() == 0 or selected.size() > 0) {
 657           glColor4f(1.0, 0.0, 0.0, 1.0);
 658         } else {
 659           glColor4f(1.0, 0.7, 0.7, 0.4);
 660         }
 661       } else {
 662         if (selected_paths.size() == 0 or selected.size() > 0) {
 663           glColor4f(0.0, 0.0, 1.0, 1.0);
 664         } else {
 665           glColor4f(0.7, 0.7, 1.0, 0.4);
 666         }
 667         /*
 668         if (selected_paths.size() == 0 or selected.size() > 0) {
 669           glColor3f(0.0, 0.0, 1.0);
 670         } else {
 671           glColor3f(0.8, 0.8, 1.0);
 672         }
 673         */
 674       }
 675       // save the multipart name for our segment
 676       glPushName(path_index); glPushName(key.first); glPushName(key.second);
 677       glBegin(GL_LINES);
 678       float seq_start_x = track_container[path_index].x();
 679       float seq_end_x = track_container[path_index+1].x();
 680       glVertex3f(s.start.x + seq_start_x, s.start.y, -1);
 681       glVertex3f(s.end.x   + seq_end_x  , s.end.y, -1);
 682       glEnd();
 683       // clear the names
 684       glPopName(); glPopName(); glPopName();
 685     }
 686   }
 687   glDepthMask(GL_TRUE);
 688   glDisable(GL_BLEND);
 689 }