alg/mussa.cpp

   1 //  This file is part of the Mussa source distribution.
   2 //  http://mussa.caltech.edu/
   3 //  Contact author: Tristan  De Buysscher, tristan@caltech.edu
   4
   5 // This program and all associated source code files are Copyright (C) 2005
   6 // the California Institute of Technology, Pasadena, CA, 91125 USA.  It is
   7 // under the GNU Public License; please see the included LICENSE.txt
   8 // file for more information, or contact Tristan directly.
   9
  10
  11 //                        ----------------------------------------
  12 //                          ---------- mussa_class.cc -----------
  13 //                        ----------------------------------------
  14 #include <boost/filesystem/operations.hpp>
  15 #include <boost/filesystem/fstream.hpp>
  16 namespace fs = boost::filesystem;
  17
  18 #include <iostream>
  19 #include <sstream>
  20
  21 #include "mussa_exceptions.hpp"
  22 #include "alg/mussa.hpp"
  23 #include "alg/flp.hpp"
  24
  25 using namespace std;
  26
  27 Mussa::Mussa()
  28 {
  29   clear();
  30 }
  31
  32 Mussa::Mussa(const Mussa& m)
  33   : analysis_name(m.analysis_name),
  34     window(m.window),
  35     threshold(m.threshold),
  36     soft_thres(m.soft_thres),
  37     ana_mode(m.ana_mode),
  38     win_append(m.win_append),
  39     thres_append(m.thres_append),
  40     motif_sequences(m.motif_sequences),
  41     color_mapper(m.color_mapper)
  42 {
  43 }
  44
  45 // set all parameters to null state
  46 void
  47 Mussa::clear()
  48 {
  49   analysis_name = "";
  50   ana_mode = TransitiveNway;
  51   window = 0;
  52   threshold = 0;
  53   soft_thres = 0;
  54   win_append = false;
  55   thres_append = false;
  56   motif_sequences.clear();
  57   color_mapper.clear();
  58 }
  59
  60 // these 5 simple methods manually set the parameters for doing an analysis
  61 // used so that the gui can take input from user and setup the analysis
  62 // note - still need a set_append(bool, bool) method...
  63 void
  64 Mussa::set_name(string a_name)
  65 {
  66   analysis_name = a_name;
  67 }
  68
  69 string Mussa::get_name()
  70 {
  71   return analysis_name;
  72 }
  73
  74 int
  75 Mussa::size() const
  76 {
  77   if (the_seqs.size() > 0)
  78     return the_seqs.size();
  79   else
  80     return 0;
  81 }
  82
  83 void
  84 Mussa::set_window(int a_window)
  85 {
  86   window = a_window;
  87 }
  88
  89 int Mussa::get_window() const
  90 {
  91   return window;
  92 }
  93
  94 void
  95 Mussa::set_threshold(int a_threshold)
  96 {
  97   threshold = a_threshold;
  98   if (a_threshold > soft_thres)
  99     soft_thres = a_threshold;
 100 }
 101
 102 int Mussa::get_threshold() const
 103 {
 104   return threshold;
 105 }
 106
 107 void
 108 Mussa::set_soft_threshold(int new_threshold)
 109 {
 110   if (new_threshold < threshold) {
 111     soft_thres = threshold;
 112   } else if (new_threshold > window) {
 113     soft_thres = window;
 114   } else {
 115     soft_thres = new_threshold;
 116   }
 117 }
 118
 119 int Mussa::get_soft_threshold() const
 120 {
 121   return soft_thres;
 122 }
 123
 124 void
 125 Mussa::set_analysis_mode(enum analysis_modes new_ana_mode)
 126 {
 127   ana_mode = new_ana_mode;
 128 }
 129
 130 enum Mussa::analysis_modes Mussa::get_analysis_mode() const
 131 {
 132   return ana_mode;
 133 }
 134
 135 string Mussa::get_analysis_mode_name() const
 136 {
 137   switch (ana_mode)
 138   {
 139   case TransitiveNway:
 140     return string("Transitive");
 141     break;
 142   case RadialNway:
 143     return string("Radial");
 144     break;
 145   case EntropyNway:
 146     return string("Entropy");
 147     break;
 148   case RecursiveNway:
 149     return string("[deprecated] Recursive");
 150     break;
 151   default:
 152     throw runtime_error("invalid analysis mode type");
 153     break;
 154   }
 155 }
 156
 157 const NwayPaths& Mussa::paths() const
 158 {
 159   return the_paths;
 160 }
 161
 162 //template <class IteratorT>
 163 //void Mussa::createLocalAlignment(IteratorT begin, IteratorT end)
 164 void Mussa::createLocalAlignment(std::list<ConservedPath>::iterator begin,
 165                                  std::list<ConservedPath>::iterator end,
 166                                  std::list<ConservedPath::path_type>& result,
 167                                  std::list<std::vector<bool> >& reversed)
 168 {
 169   const vector<Sequence>& raw_seq = the_seqs;
 170   ConservedPath::path_type aligned_path;
 171   size_t i2, i3;
 172   int x_start, x_end;
 173   int window_length, win_i;
 174   int rc_1 = 0;
 175   int rc_2 = 0;
 176   vector<bool> rc_list;
 177   bool full_match;
 178   vector<bool> matched;
 179   int align_counter;
 180
 181   align_counter = 0;
 182   for(std::list<ConservedPath>::iterator pathz_i=begin; pathz_i != end; ++pathz_i)
 183   {
 184     ConservedPath& a_path = *pathz_i;
 185     window_length = a_path.window_size;
 186     // determine which parts of the path are RC relative to first species
 187     rc_list = a_path.reverseComplimented();
 188
 189     // loop over each bp in the conserved region for all sequences
 190     for(win_i = 0; win_i < window_length; win_i++)
 191     {
 192       aligned_path.clear();
 193       // determine which exact base pairs match between the sequences
 194       full_match = true;
 195       for(i2 = 0; i2 < a_path.size()-1; i2++)
 196       {
 197         // assume not rc as most likely, adjust below
 198         rc_1 = 0;
 199         rc_2 = 0;
 200         // no matter the case, any RC node needs adjustments
 201         if (a_path[i2] < 0)
 202           rc_1 = window_length-1;
 203         if (a_path[i2+1] < 0)
 204           rc_2 = window_length-1;
 205
 206         x_start = (abs(a_path[i2]-rc_1+win_i));
 207         x_end =   (abs(a_path[i2+1]-rc_2+win_i));
 208
 209         // RC case handling
 210         // ugh, and xor...only want rc coloring if just one of the nodes is rc
 211         // if both nodes are rc, then they are 'normal' relative to each other
 212         if((rc_list[i2] || rc_list[i2+1] )&&!(rc_list[i2] && rc_list[i2+1]))
 213         { //the hideous rc matching logic - not complex, but annoying
 214           if(!(( (raw_seq[i2][x_start]=='A')&&(raw_seq[i2+1][x_end]=='T')) ||
 215                 ((raw_seq[i2][x_start]=='T')&&(raw_seq[i2+1][x_end]=='A')) ||
 216                 ((raw_seq[i2][x_start]=='G')&&(raw_seq[i2+1][x_end]=='C')) ||
 217                 ((raw_seq[i2][x_start]=='C')&&(raw_seq[i2+1][x_end]=='G'))) )
 218           {
 219             full_match = false;
 220           } else {
 221             aligned_path.push_back(x_start);
 222           }
 223         }
 224         else
 225         { // forward match
 226           if (!( (raw_seq[i2][x_start] == raw_seq[i2+1][x_end]) &&
 227                 (raw_seq[i2][x_start] != 'N') &&
 228                 (raw_seq[i2+1][x_end] != 'N') ) ) {
 229             full_match = false;
 230           } else {
 231             aligned_path.push_back(x_start);
 232           }
 233         }
 234       }
 235       // grab the last part of our path, assuming we matched
 236       if (full_match)
 237         aligned_path.push_back(x_end);
 238
 239       if (aligned_path.size() == a_path.size()) {
 240         result.push_back(aligned_path);
 241         reversed.push_back(rc_list);
 242       }
 243     }
 244     align_counter++;
 245   }
 246 }
 247
 248
 249 // takes a string and sets it as the next seq
 250 void
 251 Mussa::add_a_seq(string a_seq)
 252 {
 253   Sequence aSeq;
 254
 255   aSeq.set_seq(a_seq);
 256   the_seqs.push_back(aSeq);
 257 }
 258
 259 const vector<Sequence>&
 260 Mussa::sequences() const
 261 {
 262   return the_seqs;
 263 }
 264
 265 void Mussa::load_sequence(fs::path seq_file, fs::path annot_file,
 266                           int fasta_index, int sub_seq_start, int sub_seq_end)
 267 {
 268   Sequence aseq;
 269   aseq.load_fasta(seq_file, fasta_index, sub_seq_start, sub_seq_end);
 270   if ( not annot_file.empty() ) {
 271     aseq.load_annot(annot_file, sub_seq_start, sub_seq_end);
 272   }
 273   the_seqs.push_back(aseq);
 274 }
 275
 276 void
 277 Mussa::load_mupa_file(fs::path para_file_path)
 278 {
 279   fs::ifstream para_file;
 280   string file_data_line;
 281   string param, value;
 282   fs::path annot_file;
 283   int split_index, fasta_index;
 284   int sub_seq_start, sub_seq_end;
 285   bool seq_params, did_seq;
 286   string err_msg;
 287   bool parsing_path;
 288   string::size_type new_index, dir_index;
 289
 290   // initialize values
 291   clear();
 292
 293   // if file was opened, read the parameter values
 294   if (fs::exists(para_file_path))
 295   {
 296     para_file.open(para_file_path, ios::in);
 297
 298     // what directory is the mupa file in?
 299     fs::path file_path_base = para_file_path.branch_path();
 300
 301     // setup loop by getting file's first line
 302     getline(para_file,file_data_line);
 303     split_index = file_data_line.find(" ");
 304     param = file_data_line.substr(0,split_index);
 305     value = file_data_line.substr(split_index+1);
 306
 307     while (!para_file.eof())
 308     {
 309       did_seq = false;
 310       if (param == "ANA_NAME")
 311         analysis_name = value;
 312       else if (param == "APPEND_WIN")
 313         win_append = true;
 314       else if (param == "APPEND_THRES")
 315         thres_append = true;
 316       else if (param == "SEQUENCE_NUM")
 317         ; // ignore sequence_num now
 318       else if (param == "WINDOW")
 319         window = atoi(value.c_str());
 320       else if (param == "THRESHOLD")
 321         threshold = atoi(value.c_str());
 322       else if (param == "SEQUENCE")
 323       {
 324         fs::path seq_file = file_path_base / value;
 325         //cout << "seq_file_name " << seq_files.back() << endl;
 326         fasta_index = 1;
 327         annot_file = "";
 328         sub_seq_start = 0;
 329         sub_seq_end = 0;
 330         seq_params = true;
 331
 332         while ((!para_file.eof()) && seq_params)
 333         {
 334           getline(para_file,file_data_line);
 335           split_index = file_data_line.find(" ");
 336           param = file_data_line.substr(0,split_index);
 337           value = file_data_line.substr(split_index+1);
 338
 339           if (param == "FASTA_INDEX")
 340             fasta_index = atoi(value.c_str());
 341           else if (param == "ANNOTATION")
 342             annot_file = file_path_base / value;
 343           else if (param == "SEQ_START")
 344             sub_seq_start = atoi(value.c_str());
 345           else if (param == "SEQ_END")
 346           {
 347             sub_seq_end = atoi(value.c_str());
 348           }
 349           //ignore empty lines or that start with '#'
 350           else if ((param == "") || (param == "#")) {}
 351           else seq_params = false;
 352         }
 353         load_sequence(seq_file, annot_file, fasta_index, sub_seq_start,
 354                       sub_seq_end);
 355         did_seq = true;
 356       }
 357       //ignore empty lines or that start with '#'
 358       else if ((param == "") || (param == "#")) {}
 359       else
 360       {
 361         clog << "Illegal/misplaced mussa parameter in file\n";
 362         clog << param << "\n";
 363       }
 364
 365       if (!did_seq)
 366       {
 367         getline(para_file,file_data_line);
 368         split_index = file_data_line.find(" ");
 369         param = file_data_line.substr(0,split_index);
 370         value = file_data_line.substr(split_index+1);
 371         did_seq = false;
 372       }
 373     }
 374
 375     para_file.close();
 376
 377     soft_thres = threshold;
 378     //cout << "nway mupa: analysis_name = " << analysis_name
 379     //     << " window = " << window
 380     //     << " threshold = " << threshold << endl;
 381   }
 382   // no file was loaded, signal error
 383   else
 384   {
 385     throw mussa_load_error("Config File: " + para_file_path.string() + " not found");
 386   }
 387 }
 388
 389
 390 void
 391 Mussa::analyze(int w, int t, enum Mussa::analysis_modes the_ana_mode, double new_ent_thres)
 392 {
 393   time_t t1, t2, begin, end;
 394   double seqloadtime, seqcomptime, nwaytime, savetime, totaltime;
 395
 396   begin = time(NULL);
 397
 398   ana_mode = the_ana_mode;
 399   ent_thres = new_ent_thres;
 400   if (w > 0)
 401     window  = w;
 402   if (t > 0)
 403   {
 404     threshold = t;
 405     soft_thres = t;
 406   }
 407
 408   t1 = time(NULL);
 409
 410   if (the_seqs.size() < 2) {
 411     throw mussa_analysis_error("you need to have at least 2 sequences to "
 412                                "do an analysis.");
 413   }
 414   //cout << "nway ana: seq_num = " << the_seqs.size() << endl;
 415
 416   t2 = time(NULL);
 417   seqloadtime = difftime(t2, t1);
 418
 419   t1 = time(NULL);
 420   seqcomp();
 421   t2 = time(NULL);
 422   seqcomptime = difftime(t2, t1);
 423
 424
 425   t1 = time(NULL);
 426   the_paths.setup(window, threshold);
 427   nway();
 428   t2 = time(NULL);
 429   nwaytime = difftime(t2, t1);
 430
 431   t1 = time(NULL);
 432   save();
 433   t2 = time(NULL);
 434   savetime = difftime(t2, t1);
 435
 436   end = time(NULL);
 437   totaltime = difftime(end, begin);
 438
 439
 440   //cout << "seqload\tseqcomp\tnway\tsave\ttotal\n";
 441   //cout << seqloadtime << "\t";
 442   //cout << seqcomptime << "\t";
 443   //cout << nwaytime << "\t";
 444   //cout << savetime << "\t";
 445   //cout << totaltime << "\n";
 446 }
 447
 448 void
 449 Mussa::seqcomp()
 450 {
 451   vector<int> seq_lens;
 452   vector<FLPs> empty_FLP_vector;
 453   FLPs dummy_comp;
 454   string save_file_string;
 455
 456   empty_FLP_vector.clear();
 457   for(vector<Sequence>::size_type i = 0; i < the_seqs.size(); i++)
 458   {
 459     all_comps.push_back(empty_FLP_vector);
 460     for(vector<Sequence>::size_type i2 = 0; i2 < the_seqs.size(); i2++)
 461       all_comps[i].push_back(dummy_comp);
 462   }
 463   for(vector<Sequence>::size_type i = 0; i < the_seqs.size(); i++)
 464     seq_lens.push_back(the_seqs[i].size());
 465
 466   for(vector<Sequence>::size_type i = 0; i < the_seqs.size(); i++)
 467     for(vector<Sequence>::size_type i2 = i+1; i2 < the_seqs.size(); i2++)
 468     {
 469       //cout << "seqcomping: " << i << " v. " << i2 << endl;
 470       all_comps[i][i2].setup(window, threshold);
 471       all_comps[i][i2].seqcomp(the_seqs[i].get_seq(), the_seqs[i2].get_seq(), false);
 472       all_comps[i][i2].seqcomp(the_seqs[i].get_seq(),the_seqs[i2].rev_comp(),true);
 473     }
 474 }
 475
 476 void
 477 Mussa::nway()
 478 {
 479   vector<string> some_Seqs;
 480
 481   the_paths.set_soft_threshold(soft_thres);
 482
 483   if (ana_mode == TransitiveNway) {
 484     the_paths.trans_path_search(all_comps);
 485   }
 486   else if (ana_mode == RadialNway) {
 487     the_paths.radiate_path_search(all_comps);
 488   }
 489   else if (ana_mode == EntropyNway)
 490   {
 491     //unlike other methods, entropy needs to look at the sequence at this stage
 492     some_Seqs.clear();
 493     for(vector<Sequence>::size_type i = 0; i < the_seqs.size(); i++)
 494     {
 495       some_Seqs.push_back(the_seqs[i].get_seq());
 496     }
 497
 498     the_paths.setup_ent(ent_thres, some_Seqs); // ent analysis extra setup
 499     the_paths.entropy_path_search(all_comps);
 500   }
 501
 502   // old recursive transitive analysis function
 503   else if (ana_mode == RecursiveNway)
 504     the_paths.find_paths_r(all_comps);
 505
 506   the_paths.simple_refine();
 507 }
 508
 509 void
 510 Mussa::save()
 511 {
 512   string save_name;
 513   fs::path flp_filepath;
 514   fs::fstream save_file;
 515   ostringstream append_info;
 516   int dir_create_status;
 517
 518   if (not analysis_name.empty()) {
 519     // not sure why, but gotta close file each time since can't pass
 520     // file streams
 521     save_name = analysis_name;
 522
 523     // gotta do bit with adding win & thres if to be appended
 524     if (win_append)
 525     {
 526       append_info.str("");
 527       append_info <<  "_w" << window;
 528       save_name += append_info.str();
 529     }
 530
 531     if (thres_append)
 532     {
 533       append_info.str("");
 534       append_info <<  "_t" << threshold;
 535       save_name += append_info.str();
 536     }
 537     fs::path save_path( save_name);
 538
 539     if (not fs::exists(save_path)) {
 540       fs::create_directory(save_path);
 541     }
 542     // save sequence and annots to a special mussa file
 543     save_file.open(save_path / (save_name+".museq"), ios::out);
 544     save_file << "<Mussa_Sequence>" << endl;
 545
 546     for(vector<Sequence>::size_type i = 0; i < the_seqs.size(); i++)
 547     {
 548       the_seqs[i].save(save_file);
 549     }
 550
 551     save_file << "</Mussa_Sequence>" << endl;
 552     save_file.close();
 553
 554     // save nway paths to its mussa save file
 555     the_paths.save(save_path / (save_name + ".muway"));
 556
 557     for(vector<Sequence>::size_type i = 0; i < the_seqs.size(); i++)
 558       for(vector<Sequence>::size_type i2 = i+1; i2 < the_seqs.size(); i2++)
 559       {
 560         append_info.str("");
 561         append_info <<  "_sp_" << i << "v" << i2;
 562         all_comps[i][i2].save(save_path/(save_name+append_info.str()+".flp"));
 563       }
 564   }
 565 }
 566
 567 void
 568 Mussa::save_muway(fs::path save_path)
 569 {
 570   the_paths.save(save_path);
 571 }
 572
 573 void
 574 Mussa::load(fs::path ana_file)
 575 {
 576   int i, i2;
 577   fs::path file_path_base;
 578   fs::path a_file_path;
 579   fs::path ana_path(ana_file);
 580   bool parsing_path;
 581   Sequence tmp_seq;
 582   string err_msg;
 583   ostringstream append_info;
 584   vector<FLPs> empty_FLP_vector;
 585   FLPs dummy_comp;
 586
 587   //cout << "ana_file name " << ana_file.string() << endl;
 588   analysis_name = ana_path.leaf();
 589   //cout << " ana_name " << analysis_name << endl;
 590   file_path_base =  ana_path.branch_path() / analysis_name;
 591   a_file_path = file_path_base / (analysis_name + ".muway");
 592   //cout << " loading museq: " << a_file_path.string() << endl;
 593   the_paths.load(a_file_path);
 594   // perhaps this could be more elegent, but at least this'll let
 595   // us know what our threshold and window sizes were when we load a muway
 596   window = the_paths.get_window();
 597   threshold = the_paths.get_threshold();
 598   soft_thres = threshold;
 599
 600   int seq_num = the_paths.sequence_count();
 601
 602   a_file_path = file_path_base / (analysis_name + ".museq");
 603
 604   // this is a bit of a hack due to C++ not acting like it should with files
 605   for (i = 1; i <= seq_num; i++)
 606   {
 607     tmp_seq.clear();
 608     //cout << "mussa_class: loading museq frag... " << a_file_path.string() << endl;
 609     tmp_seq.load_museq(a_file_path, i);
 610     the_seqs.push_back(tmp_seq);
 611   }
 612
 613   empty_FLP_vector.clear();
 614   for(i = 0; i < seq_num; i++)
 615   {
 616     all_comps.push_back(empty_FLP_vector);
 617     for(i2 = 0; i2 < seq_num; i2++)
 618       all_comps[i].push_back(dummy_comp);
 619   }
 620
 621   for(i = 0; i < seq_num; i++)
 622   {
 623     for(i2 = i+1; i2 < seq_num; i2++)
 624     {
 625       append_info.str("");
 626       append_info << analysis_name <<  "_sp_" << i << "v" << i2 << ".flp";
 627       //cout << append_info.str() << endl;
 628       a_file_path = file_path_base / append_info.str();
 629       //cout << "path " << a_file_path.string() << endl;
 630       all_comps[i][i2].load(a_file_path);
 631       //cout << "real size = " << all_comps[i][i2].size() << endl;
 632     }
 633   }
 634 }
 635
 636
 637 void
 638 Mussa::save_old()
 639 {
 640   fs::fstream save_file;
 641
 642   save_file.open(analysis_name, ios::out);
 643
 644   for(vector<Sequence>::size_type i = 0; i < the_seqs.size(); i++)
 645     save_file << the_seqs[i].get_seq() << endl;
 646
 647   save_file << window << endl;
 648   save_file.close();
 649   //note more complex eventually since analysis_name may need to have
 650   //window size, threshold and other stuff to modify it...
 651   the_paths.save_old(analysis_name);
 652 }
 653
 654
 655 void
 656 Mussa::load_old(char * load_file_path, int s_num)
 657 {
 658   fstream save_file;
 659   string file_data_line;
 660   int i, space_split_i, comma_split_i;
 661   vector<int> loaded_path;
 662   string node_pair, node;
 663   Sequence a_seq;
 664
 665   int seq_num = s_num;
 666   the_paths.setup(0, 0);
 667   save_file.open(load_file_path, ios::in);
 668
 669   // currently loads old mussa format
 670
 671   // get sequences
 672   for(i = 0; i < seq_num; i++)
 673   {
 674     getline(save_file, file_data_line);
 675     a_seq.set_seq(file_data_line);
 676     the_seqs.push_back(a_seq);
 677   }
 678
 679   // get window size
 680   getline(save_file, file_data_line);
 681   window = atoi(file_data_line.c_str());
 682   // get paths
 683
 684   while (!save_file.eof())
 685   {
 686     loaded_path.clear();
 687     getline(save_file, file_data_line);
 688     if (file_data_line != "")
 689     for(i = 0; i < seq_num; i++)
 690     {
 691       space_split_i = file_data_line.find(" ");
 692       node_pair = file_data_line.substr(0,space_split_i);
 693       //cout << "np= " << node_pair;
 694       comma_split_i = node_pair.find(",");
 695       node = node_pair.substr(comma_split_i+1);
 696       //cout << "n= " << node << " ";
 697       loaded_path.push_back(atoi (node.c_str()));
 698       file_data_line = file_data_line.substr(space_split_i+1);
 699     }
 700     //cout << endl;
 701     // FIXME: do we have any information about what the threshold should be?
 702     the_paths.add_path(0, loaded_path);
 703   }
 704   save_file.close();
 705
 706   //the_paths.save("tmp.save");
 707 }
 708
 709 void Mussa::add_motifs(const vector<string>& motifs,
 710                        const vector<Color>& colors)
 711 {
 712   if (motifs.size() != colors.size()) {
 713     throw mussa_error("motif and color vectors must be the same size");
 714   }
 715
 716   for(size_t i = 0; i != motifs.size(); ++i)
 717   {
 718     motif_sequences.insert(motifs[i]);
 719     color_mapper.appendInstanceColor("motif", motifs[i], colors[i]);
 720   }
 721   update_sequences_motifs();
 722 }
 723
 724 // I mostly split the ifstream out so I can use a stringstream to test it.
 725 void Mussa::load_motifs(std::istream &in)
 726 {
 727   string seq;
 728   float red;
 729   float green;
 730   float blue;
 731
 732   while(in.good())
 733   {
 734     in >> seq >> red >> green >> blue;
 735     // if we couldn't read this line 'cause we're like at the end of the file
 736     // try to exit the loop
 737     if (!in.good())
 738       break;
 739     try {
 740       seq = Sequence::motif_normalize(seq);
 741     } catch(motif_normalize_error e) {
 742       clog << "unable to parse " << seq << " skipping" << endl;
 743       clog << e.what() << endl;
 744       continue;
 745     }
 746     if (red < 0.0 or red > 1.0) {
 747       clog << "invalid red value " << red << ". must be in range [0..1]"
 748            << endl;
 749       continue;
 750     }
 751     if (green < 0.0 or green > 1.0) {
 752       clog << "invalid green value " << green << ". must be in range [0..1]"
 753            << endl;
 754       continue;
 755     }
 756     if (blue < 0.0 or blue > 1.0) {
 757       clog << "invalid blue value " << blue << ". must be in range [0..1]"
 758            << endl;
 759       continue;
 760     }
 761     if (motif_sequences.find(seq) == motif_sequences.end()) {
 762       // sequence wasn't found
 763       motif_sequences.insert(seq);
 764       Color c(red, green, blue);
 765       color_mapper.appendInstanceColor("motif", seq, c);
 766     } else {
 767       clog << "sequence " << seq << " was already defined skipping"
 768            << endl;
 769       continue;
 770     }
 771   }
 772   update_sequences_motifs();
 773 }
 774
 775 void Mussa::load_motifs(fs::path filename)
 776 {
 777   fs::ifstream f;
 778   f.open(filename, ifstream::in);
 779   load_motifs(f);
 780 }
 781
 782 void Mussa::update_sequences_motifs()
 783 {
 784   // once we've loaded all the motifs from the file,
 785   // lets attach them to the sequences
 786   for(vector<Sequence>::iterator seq_i = the_seqs.begin();
 787       seq_i != the_seqs.end();
 788       ++seq_i)
 789   {
 790     // clear out old motifs
 791     seq_i->clear_motifs();
 792     // for all the motifs in our set, attach them to the current sequence
 793     for(set<string>::iterator motif_i = motif_sequences.begin();
 794         motif_i != motif_sequences.end();
 795         ++motif_i)
 796     {
 797       seq_i->add_motif(*motif_i);
 798     }
 799   }
 800 }
 801
 802 const set<string>& Mussa::motifs() const
 803 {
 804   return motif_sequences;
 805 }
 806
 807 AnnotationColors& Mussa::colorMapper()
 808 {
 809   return color_mapper;
 810 }