make Sequence a subclass of std::string
[mussa.git] / 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(AnnotationColors& cm)
289 {
290   color_mapper = cm;
291 }
292
293 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   GlSequence gs(s, color_mapper);
317   push_sequence(gs);
318 }
319
320 void GlSeqBrowser::push_sequence(GlSequence &gs)
321 {
322   clear_links();
323   track_container.push_back(gs);
324   update_layout();
325   if (track_container.size() > 1)
326     path_segments.push_back(pair_segment_map());
327 }
328
329 const std::vector<GlSequence>& GlSeqBrowser::sequences() const
330 {
331   return track_container;
332 }
333
334 void GlSeqBrowser::clear_links()
335 {
336   path_segments.clear();
337   for (int i = track_container.size()-1; i > 0; --i)
338   {
339     path_segments.push_back(pair_segment_map());
340   }
341   pathid = 0;
342 }
343
344 void 
345 GlSeqBrowser::link(const vector<int>& path, const vector<bool>& rc, int )
346 {
347   if (path.size() < 2) {
348     // should i throw an error instead?
349     return;
350   }
351   if (path.size() != track_container.size() ) {
352     stringstream msg;
353     msg << "Path size [" << path.size() << "] and track size [" 
354         << track_container.size() << "] don't match" << endl;
355     throw mussa_error(msg.str());
356   }
357   if (path.size() != rc.size()) {
358     throw runtime_error("path and reverse compliment must be the same length");
359   }
360   vector<int>::const_iterator path_i = path.begin();
361   vector<bool>::const_iterator rc_i = rc.begin();
362   int track_i = 0;
363   int prev_x = *path_i; ++path_i;
364   bool prev_rc = *rc_i; ++rc_i;
365   while (path_i != path.end() and rc_i != rc.end())
366   {
367     segment_key p(prev_x, *path_i);
368     pair_segment_map::iterator found_segment = path_segments[track_i].find(p);
369     if (found_segment == path_segments[track_i].end()) {
370       // not already found
371       float y1 = track_container[track_i].y();
372             y1 -= track_container[track_i].height()/2;
373       float y2 = track_container[track_i+1].y();
374             y2 += track_container[track_i+1].height()/2;
375       
376       bool rcFlag = (prev_rc or *rc_i) and !(prev_rc and *rc_i);
377       Segment s(prev_x, y1, *path_i, y2, rcFlag);
378       s.path_ids.insert(pathid);
379       path_segments[track_i][p] = s;
380     } else {
381       //found
382       found_segment->second.path_ids.insert(pathid);
383     }
384     prev_x = *path_i;
385     prev_rc = *rc_i;
386     ++track_i;
387     ++path_i;
388     ++rc_i;
389   }
390   // pathid is reset by push_sequence
391   ++pathid;
392 }
393
394 const set<int>& GlSeqBrowser::selectedPaths() const
395 {
396   return selected_paths;
397 }
398
399 //! copy sequence from selected track using formating function
400 void GlSeqBrowser::copySelectedTracks(std::string& copy_buffer, 
401                                       format_track formatter)
402 {
403   copy_buffer.clear();
404
405   for(selected_track_iterator track_i = selected_tracks.begin();
406       track_i != selected_tracks.end();
407       ++track_i)
408   {
409     int track_index = track_i->track_id;
410     if (track_index >= track_container.size()) {
411       // should this be an exception instead?
412       clog << "track " << track_index << " > " << track_container.size() 
413            << endl;
414     } else {
415       // we should be safe
416       const Sequence& seq = track_container[track_index].sequence();
417       copy_buffer += formatter(seq, track_i->left, track_i->right);
418     }
419   }
420 }
421
422 //! copy sequence from selected tracks as FASTA sequences
423 void GlSeqBrowser::copySelectedTracksAsFasta(std::string& copy_buffer)
424 {
425   struct AsFasta {
426     static string formatter(const Sequence& seq, int left, int right)
427     {
428       stringstream s;
429       s << ">" << seq.get_header() 
430         << "|" << "subregion=" << left << "-" << right+1
431         << std::endl
432         << seq.subseq(left, right-left+1) << std::endl;
433       return s.str();
434     }
435   };
436   copySelectedTracks(copy_buffer, AsFasta::formatter);
437 }
438
439 //! copy sequence from selected tracks as new sequences
440 /*
441 void GlSeqBrowser::copySelectedTracksAsSequence(std::string& copy_buffer)
442 {
443   struct AsSequence {
444     static string formatter(const Sequence& seq, int left, int right)
445     {
446       stringstream s;
447       s << ">" << seq.get_header() 
448         << "|" << "subregion=" << left << "-" << right+1
449         << std::endl
450         << seq.subseq(left, right-left+1) << std::endl;
451       return s.str();
452     }
453   };
454   copySelectedTracks(copy_buffer, AsFasta::formatter);
455 }
456 */
457
458 //! copy sequence from selected tracks as plain sequences
459 void GlSeqBrowser::copySelectedTracksAsString(std::string& copy_buffer)
460 {
461   struct AsString {
462     static string formatter(const Sequence& seq, int left, int right)
463     {
464       stringstream s;
465       s << seq.subseq(left, right-left+1) << std::endl;
466       return s.str();
467     }
468   };
469
470   copySelectedTracks(copy_buffer, AsString::formatter);
471 }
472
473 void GlSeqBrowser::centerOnPath(const vector<int>& paths)
474 {
475   if (paths.size() != track_container.size()) {
476     throw mussa_error("Path length didn't match the number of sequences");
477   }
478
479   for(size_t track_i = 0; track_i != track_container.size(); ++track_i)
480   {
481     // -15 = shift more to the left
482     track_container[track_i].setX((viewport_center-15) - paths[track_i]);
483   }
484 }
485
486 void GlSeqBrowser::update_viewport(float center, double new_zoom)
487 {
488   // limit how close we can get
489   if (new_zoom < 0.01) {
490     new_zoom = 0.01;
491   }
492   double new_width = (new_zoom * (float)viewport_size.x);
493   cur_ortho.left = center-new_width/2.0;
494   cur_ortho.right = center+new_width/2.0;
495 }
496
497 void GlSeqBrowser::update_layout()
498 {
499   typedef std::vector<GlSequence>::iterator glseq_itor_type;
500   float available_height = (float)cur_ortho.top - 2 * (float)border_width;
501   float max_base_pairs = 0;
502   size_t track_count = track_container.size();
503
504   if (track_count > 1) {
505     // we have several sequences
506     float track_spacing = available_height / (track_count-1);
507     float y = available_height + (float)border_width;
508     for(glseq_itor_type seq_i = track_container.begin();
509         seq_i != track_container.end();
510         ++seq_i, y-=track_spacing)
511     {
512       seq_i->setX(0);
513       seq_i->setY(y);
514       if (seq_i->length() > max_base_pairs)
515         max_base_pairs = seq_i->length();
516     }
517   } else if (track_count == 1) {
518     // center the single track
519     glseq_itor_type seq_i = track_container.begin();
520     seq_i->setX(0);
521     seq_i->setY(viewport_size.x /2);
522     max_base_pairs = seq_i->length();
523   } else {
524     // nothing to do as we're empty
525     return;
526   }
527   cur_ortho.right = max_base_pairs + border_width;
528   cur_ortho.left = -border_width;
529   cur_ortho.top = viewport_size.x;
530   cur_ortho.bottom = 0;
531   viewport_center = (cur_ortho.width()/2) + cur_ortho.left;
532   zoomOut();
533 }
534
535 void GlSeqBrowser::draw() const
536 {
537   glMatrixMode(GL_MODELVIEW);
538   glInitNames();
539   glPushName(MussaSegment);
540   draw_segments();
541   glLoadName(MussaTrack);
542   draw_tracks();
543   glPopName();
544   // a selection shouldn't have a glName associated with it
545   draw_selection();
546 }
547
548 void GlSeqBrowser::draw_selection() const
549 {
550   // draw selection box
551   glEnable(GL_BLEND);
552   glDepthMask(GL_FALSE);
553   if (selectedMode) {
554     glColor4f(0.6, 0.6, 0.6, 0.9);
555     glRectf(selectedRegion.left, selectedRegion.top, 
556             selectedRegion.right, selectedRegion.bottom);
557   }
558   glDepthMask(GL_TRUE);
559   glDisable(GL_BLEND);
560 }
561
562 void GlSeqBrowser::draw_tracks() const
563 {
564   for(size_t track_i = 0; track_i != track_container.size(); ++track_i)
565   {
566     glPushName(track_i);
567     track_container[track_i].draw(cur_ortho.left, cur_ortho.right);
568     glPopName();
569   }
570 }
571
572 void GlSeqBrowser::draw_segments() const
573 {
574   glLineWidth(1);
575   // each vector contains path_segment_maps of all the connections
576   // between this track and the next
577   path_segment_map_vector::const_iterator psmv_i;
578   for(psmv_i = path_segments.begin();
579       psmv_i != path_segments.end();
580       ++psmv_i)
581   {
582     path_segment_map_vector::difference_type path_index;
583     path_index = psmv_i - path_segments.begin();
584     // these maps contain the pair index (used so we dont keep drawing the
585     // same segment) and the actual segment structure.
586     pair_segment_map::const_iterator psm_i;
587     for(psm_i = psmv_i->begin();
588         psm_i != psmv_i->end();
589         ++psm_i)
590     {
591       // grab the index into our segment map
592       const segment_key& key = psm_i->first;
593       // the second element of our map pair is a segment
594       const Segment &s = psm_i->second;
595       // need to do something so we can detect our selection
596       vector<int> selected;
597       set_intersection(selected_paths.begin(), selected_paths.end(),
598                        s.path_ids.begin(), s.path_ids.end(),
599                        back_inserter(selected));
600
601       if (not s.reversed) {
602         if (selected_paths.size() == 0 or selected.size() > 0) {
603           glColor3f(1.0, 0.0, 0.0);
604         } else {
605           glColor3f(1.0, 0.8, 0.8);
606         }
607       } else { 
608         if (selected_paths.size() == 0 or selected.size() > 0) {
609           glColor3f(0.0, 0.0, 1.0);
610         } else {
611           glColor3f(0.8, 0.8, 1.0);
612         }
613         /*
614         if (selected_paths.size() == 0 or selected.size() > 0) {
615           glColor3f(0.0, 0.0, 1.0);
616         } else {
617           glColor3f(0.8, 0.8, 1.0);
618         }
619         */
620       }
621       // save the multipart name for our segment
622       glPushName(path_index); glPushName(key.first); glPushName(key.second);
623       glBegin(GL_LINES);
624       float seq_start_x = track_container[path_index].x();
625       float seq_end_x = track_container[path_index+1].x();
626       glVertex3f(s.start.x + seq_start_x, s.start.y, -1);
627       glVertex3f(s.end.x   + seq_end_x  , s.end.y, -1);
628       glEnd();
629       // clear the names
630       glPopName(); glPopName(); glPopName();
631     }
632   }
633 }