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