mussa_class.cc

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