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