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