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