alg/mussa_class.cxx

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