mussa_class.cc

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