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