mussa_class.cc

   1 //                        ----------------------------------------
   2 //                          ---------- mussa_class.cc -----------
   3 //                        ----------------------------------------
   4
   5 #include "mussa_class.hh"
   6
   7 // doesn't do neg ints...
   8 string
   9 int_to_str(int an_int)
  10 {
  11   string converted_int;
  12   int remainder;
  13
  14   converted_int = "";
  15
  16   if (an_int == 0)
  17     converted_int = "0";
  18
  19   while (an_int != 0)
  20   {
  21     remainder = an_int % 10;
  22
  23     if (remainder == 0)
  24       converted_int = "0" + converted_int;
  25     else if (remainder == 1)
  26       converted_int = "1" + converted_int;
  27     else if (remainder == 2)
  28       converted_int = "2" + converted_int;
  29     else if (remainder == 3)
  30       converted_int = "3" + converted_int;
  31     else if (remainder == 4)
  32       converted_int = "4" + converted_int;
  33     else if (remainder == 5)
  34       converted_int = "5" + converted_int;
  35     else if (remainder == 6)
  36       converted_int = "6" + converted_int;
  37     else if (remainder == 7)
  38       converted_int = "7" + converted_int;
  39     else if (remainder == 8)
  40       converted_int = "8" + converted_int;
  41     else if (remainder == 9)
  42       converted_int = "9" + converted_int;
  43
  44     an_int = an_int / 10;
  45   }
  46
  47   return converted_int;
  48 }
  49
  50
  51 Mussa::Mussa()
  52 {
  53 }
  54
  55
  56 string
  57 Mussa::load_parameters(string para_file_path)
  58 {
  59   ifstream para_file;
  60   string file_data_line;
  61   string param, value, annot_file;
  62   int split_index, fasta_index;
  63   int sub_seq_start, sub_seq_end;
  64   bool seq_params, did_seq;
  65   string err_msg;
  66   int bogo;
  67
  68
  69   // initialize values
  70   win_append = false;
  71   thres_append = false;
  72   seq_files.clear();
  73   fasta_indices.clear();
  74   annot_files.clear();
  75   sub_seq_starts.clear();
  76   sub_seq_ends.clear();
  77
  78   para_file.open(para_file_path.c_str(), ios::in);
  79
  80   // if file was opened, read the parameter values
  81   if (para_file)
  82   {
  83     // setup loop by getting file's first line
  84     getline(para_file,file_data_line);
  85     split_index = file_data_line.find(" ");
  86     param = file_data_line.substr(0,split_index);
  87     value = file_data_line.substr(split_index+1);
  88
  89     while (!para_file.eof())
  90     {
  91       did_seq = false;
  92       if (param == "ANA_NAME")
  93         ana_name = value;
  94       else if (param == "APPEND_WIN")
  95         win_append = true;
  96       else if (param == "APPEND_THRES")
  97         thres_append = true;
  98       else if (param == "SEQUENCE_NUM")
  99         seq_num = atoi(value.c_str());
 100       else if (param == "WINDOW")
 101         window = atoi(value.c_str());
 102       else if (param == "THRESHOLD")
 103         threshold = atoi(value.c_str());
 104       else if (param == "SEQUENCE")
 105       {
 106         seq_files.push_back(value);
 107         fasta_index = 1;
 108         annot_file = "";
 109         sub_seq_start = 0;
 110         sub_seq_end = 0;
 111         seq_params = true;
 112
 113         while ((!para_file.eof()) && seq_params)
 114         {
 115           getline(para_file,file_data_line);
 116           split_index = file_data_line.find(" ");
 117           param = file_data_line.substr(0,split_index);
 118           value = file_data_line.substr(split_index+1);
 119
 120           if (param == "FASTA_INDEX")
 121             fasta_index = atoi(value.c_str());
 122           else if (param == "ANNOTATION")
 123             annot_file = value;
 124           else if (param == "SEQ_START")
 125             sub_seq_start = atoi(value.c_str());
 126           else if (param == "SEQ_END")
 127           {
 128             cout << "hey!  " << atoi(value.c_str()) << endl;
 129             sub_seq_end = atoi(value.c_str());
 130           }
 131           //ignore empty lines or that start with '#'
 132           else if ((param == "") || (param == "#")) {}
 133           else seq_params = false;
 134         }
 135
 136         fasta_indices.push_back(fasta_index);
 137         annot_files.push_back(annot_file);
 138         sub_seq_starts.push_back(sub_seq_start);
 139         sub_seq_ends.push_back(sub_seq_end);
 140         did_seq = true;
 141       }
 142       //ignore empty lines or that start with '#'
 143       else if ((param == "") || (param == "#")) {}
 144       else
 145       {
 146         cout << "Illegal/misplaced mussa parameter in file\n";
 147         cout << param << "\n";
 148       }
 149
 150       if (!did_seq)
 151       {
 152         getline(para_file,file_data_line);
 153         split_index = file_data_line.find(" ");
 154         param = file_data_line.substr(0,split_index);
 155         value = file_data_line.substr(split_index+1);
 156         did_seq = false;
 157       }
 158     }
 159
 160     para_file.close();
 161
 162     cout << "ana_name = " << ana_name << win_append << win_append << "\n";
 163     cout << "window = " << window << " threshold = " << threshold << "\n";
 164
 165     return "";
 166   }
 167
 168   // no file was loaded, signal error
 169   else
 170   {
 171     err_msg = "Config File: " + para_file_path + " not found";
 172     return err_msg;
 173   }
 174 }
 175
 176 //        if (!((param == "") || (param == "#")))
 177 //          cout << value << " = " << param << endl;
 178
 179
 180 string
 181 Mussa::analyze(string para_file_path, int w, int t)
 182 {
 183   string err_msg;
 184   time_t t1, t2, begin, end;
 185   double setuptime, seqloadtime, seqcomptime, nwaytime, savetime, totaltime;
 186   int seq_num;
 187
 188   begin = time(NULL);
 189
 190   t1 = time(NULL);
 191   err_msg = load_parameters(para_file_path);
 192   t2 = time(NULL);
 193   setuptime = difftime(t2, t1);
 194
 195   if (w > 0)
 196     window  = w;
 197   if (t > 0)
 198     threshold = t;
 199
 200   if (err_msg == "")
 201   {
 202     cout << "fee\n";
 203     t1 = time(NULL);
 204     err_msg = get_Seqs();
 205     t2 = time(NULL);
 206     seqloadtime = difftime(t2, t1);
 207
 208
 209     if (err_msg == "")
 210     {
 211       cout << "fie\n";
 212       t1 = time(NULL);
 213       seqcomp();
 214       t2 = time(NULL);
 215       seqcomptime = difftime(t2, t1);
 216
 217
 218       cout << "foe\n";
 219       t1 = time(NULL);
 220       nway();
 221       t2 = time(NULL);
 222       nwaytime = difftime(t2, t1);
 223
 224
 225       cout << "fum\n";
 226       t1 = time(NULL);
 227       save();
 228       t2 = time(NULL);
 229       savetime = difftime(t2, t1);
 230
 231       end = time(NULL);
 232       totaltime = difftime(end, begin);
 233
 234       cout << "setup\tseqload\tseqcomp\tnway\tsave\ttotal\n";
 235       cout << setuptime << "\t";
 236       cout << seqloadtime << "\t";
 237       cout << seqcomptime << "\t";
 238       cout << nwaytime << "\t";
 239       cout << savetime << "\t";
 240       cout << totaltime << "\n";
 241     }
 242     else
 243       return err_msg;
 244   }
 245   else
 246     return err_msg;
 247 }
 248
 249
 250 string
 251 Mussa::get_Seqs()
 252 {
 253   list<string>::iterator seq_files_i, annot_files_i;
 254   list<int>::iterator fasta_indices_i, seq_starts_i, seq_ends_i;
 255   Sequence aSeq;
 256   string err_msg;
 257
 258   seq_files_i = seq_files.begin();
 259   fasta_indices_i = fasta_indices.begin();
 260   annot_files_i = annot_files.begin();
 261   seq_starts_i = sub_seq_starts.begin();
 262   seq_ends_i = sub_seq_ends.begin();
 263
 264   err_msg = "";
 265
 266   while ( (seq_files_i != seq_files.end()) && (err_msg == "") )
 267           /* it should be guarenteed that each of the following exist
 268              should I bother checking, and how to deal with if not true...
 269  &&
 270           (fasta_indices_i != fasta_indices.end()) &&
 271           (annot_files_i != annot_files.end())  &&
 272           (seq_starts_i != sub_seq_starts.end())  &&
 273           (seq_ends_i != sub_seq_ends.end())         )
 274           */
 275   {
 276     err_msg = aSeq.load_fasta(*seq_files_i, *fasta_indices_i,*seq_starts_i,
 277                               *seq_ends_i);
 278     if (err_msg == "")
 279     {
 280       if (*annot_files_i != "")
 281         err_msg = aSeq.load_annot(*annot_files_i, *seq_starts_i, *seq_ends_i);
 282
 283       if (err_msg == "")
 284       {
 285         the_Seqs.push_back(aSeq);
 286         cout << aSeq.hdr() << endl;
 287         //cout << aSeq.seq() << endl;
 288         aSeq.clear();
 289         ++seq_files_i;      // advance all the iterators
 290         ++fasta_indices_i;
 291         ++annot_files_i;
 292         ++seq_starts_i;
 293         ++seq_ends_i;
 294       }
 295       else
 296         break;
 297     }
 298     else
 299       break;
 300   }
 301
 302   return err_msg;
 303 }
 304
 305
 306 void
 307 Mussa::seqcomp()
 308 {
 309   int i, i2;     // loop vars over sequences to analyze
 310   vector<int> seq_lens;
 311   vector<FLPs> empty_FLP_vector;
 312   FLPs dummy_comp;
 313   string save_file_string;
 314
 315   empty_FLP_vector.clear();
 316   for(i = 0; i < seq_num; i++)
 317   {
 318     all_comps.push_back(empty_FLP_vector);
 319     for(i2 = 0; i2 < seq_num; i2++)
 320       all_comps[i].push_back(dummy_comp);
 321   }
 322   for(i = 0; i < seq_num; i++)
 323     seq_lens.push_back(the_Seqs[i].len());
 324
 325   for(i = 0; i < seq_num; i++)
 326     for(i2 = i+1; i2 < seq_num; i2++)
 327     {
 328       cout << "seqcomping: " << i << " v. " << i2 << endl;
 329       all_comps[i][i2].setup("m", window, threshold, seq_lens[i],seq_lens[i2]);
 330       all_comps[i][i2].seqcomp(the_Seqs[i].seq(), the_Seqs[i2].seq(), false);
 331       all_comps[i][i2].seqcomp(the_Seqs[i].seq(),the_Seqs[i2].rev_comp(),true);
 332       //all_comps[i][i2].file_save(save_file_string);
 333     }
 334 }
 335
 336
 337
 338 void
 339 Mussa::nway()
 340 {
 341   the_paths.setup(seq_num, window, threshold);
 342   //the_paths.find_paths_r(all_comps);
 343   the_paths.trans_path_search(all_comps);
 344   //the_paths.radiate_path_search(all_comps);
 345   the_paths.simple_refine();
 346 }
 347
 348 void
 349 Mussa::save()
 350 {
 351   string save_path_base, save_path;
 352   fstream save_file;
 353   int i;
 354
 355   // gotta do bit with adding win & thres if to be appended - need itos
 356
 357   // not sure why, but gotta close file each time since can't pass file streams
 358
 359   save_path_base = ana_name;
 360
 361   if (win_append)
 362     save_path_base += "_w" + int_to_str(window);
 363
 364   if (thres_append)
 365     save_path_base += "_t" + int_to_str(threshold);
 366
 367   // save sequence and annots to a special mussa file
 368   save_path = save_path_base + ".museq";
 369   save_file.open(save_path.c_str(), ios::out);
 370   save_file << "<Mussa_Sequence>" << endl;
 371   //save_file.close();
 372
 373   for(i = 0; i < seq_num; i++)
 374     the_Seqs[i].save(save_file);
 375
 376   //save_file.open(save_path.c_str(), ios::app);
 377   save_file << "</Mussa_Sequence>" << endl;
 378   save_file.close();
 379
 380   // save nway paths to its mussa save file
 381   save_path = save_path_base + ".muway";
 382   the_paths.save(save_path);
 383 }
 384
 385 string
 386 Mussa::load(string ana_file)
 387 {
 388   int i;
 389   string load_file_path;
 390   Sequence tmp_seq;
 391   string err_msg;
 392
 393
 394   ana_name = ana_file;
 395   load_file_path = ana_name + ".muway";
 396   err_msg = the_paths.load(load_file_path);
 397
 398   if (err_msg == "")
 399   {
 400     seq_num = the_paths.seq_num();
 401
 402     cout << "No BAM\n";
 403
 404     load_file_path = ana_name + ".museq";
 405
 406     // this is a bit of a hack due to C++ not acting like it should with files
 407     for (i = 1; i <= seq_num; i++)
 408     {
 409       tmp_seq.clear();
 410       tmp_seq.load_museq(load_file_path, i);
 411       the_Seqs.push_back(tmp_seq);
 412     }
 413
 414     return "";
 415   }
 416   else
 417     return err_msg;
 418 }
 419
 420 /*
 421       cout << "fee\n";
 422       cout << "fie\n";
 423       cout << "foe\n";
 424       cout << "fum\n";
 425 */
 426
 427
 428
 429 void
 430 Mussa::save_old()
 431 {
 432   fstream save_file;
 433   int i;
 434
 435   save_file.open(ana_name.c_str(), ios::out);
 436
 437   for(i = 0; i < seq_num; i++)
 438     save_file << the_Seqs[i].seq() << endl;
 439
 440   save_file << window << endl;
 441   save_file.close();
 442   //note more complex eventually since ana_name may need to have
 443   //window size, threshold and other stuff to modify it...
 444   the_paths.save_old(ana_name);
 445 }
 446
 447
 448 void
 449 Mussa::load_old(char * load_file_path, int s_num)
 450 {
 451   fstream save_file;
 452   string file_data_line;
 453   int i, space_split_i, comma_split_i;
 454   vector<int> loaded_path;
 455   string node_pair, node;
 456   Sequence a_seq;
 457
 458   seq_num = s_num;
 459   the_paths.setup(seq_num, 0, 0);
 460   save_file.open(load_file_path, ios::in);
 461
 462   // currently loads old mussa format
 463
 464   // get sequences
 465   for(i = 0; i < seq_num; i++)
 466   {
 467     getline(save_file, file_data_line);
 468     a_seq.set_seq(file_data_line);
 469     the_Seqs.push_back(a_seq);
 470   }
 471
 472   // get window size
 473   getline(save_file, file_data_line);
 474   window = atoi(file_data_line.c_str());
 475   // get paths
 476
 477   while (!save_file.eof())
 478   {
 479     loaded_path.clear();
 480     getline(save_file, file_data_line);
 481     if (file_data_line != "")
 482     for(i = 0; i < seq_num; i++)
 483     {
 484       space_split_i = file_data_line.find(" ");
 485       node_pair = file_data_line.substr(0,space_split_i);
 486       //cout << "np= " << node_pair;
 487       comma_split_i = node_pair.find(",");
 488       node = node_pair.substr(comma_split_i+1);
 489       //cout << "n= " << node << " ";
 490       loaded_path.push_back(atoi (node.c_str()));
 491       file_data_line = file_data_line.substr(space_split_i+1);
 492     }
 493     //cout << endl;
 494     the_paths.add_path(loaded_path);
 495   }
 496   save_file.close();
 497
 498   //the_paths.save("tmp.save");
 499 }