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