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