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