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