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 void
  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   int bogo;
  66
  67
  68   win_append = false;
  69   thres_append = false;
  70   seq_files.clear();
  71   fasta_indices.clear();
  72   annot_files.clear();
  73   sub_seq_starts.clear();
  74   sub_seq_ends.clear();
  75
  76   para_file.open(para_file_path.c_str(), ios::in);
  77
  78   getline(para_file,file_data_line);
  79   split_index = file_data_line.find(" ");
  80   param = file_data_line.substr(0,split_index);
  81   value = file_data_line.substr(split_index+1);
  82
  83
  84   while (!para_file.eof())
  85   {
  86     did_seq = false;
  87     if (param == "ANA_NAME")
  88       ana_name = value;
  89     else if (param == "APPEND_WIN")
  90       win_append = true;
  91     else if (param == "APPEND_THRES")
  92       thres_append = true;
  93     else if (param == "SEQUENCE_NUM")
  94       seq_num = atoi(value.c_str());
  95     else if (param == "WINDOW")
  96       window = atoi(value.c_str());
  97     else if (param == "THRESHOLD")
  98       threshold = atoi(value.c_str());
  99     else if (param == "SEQUENCE")
 100     {
 101       seq_files.push_back(value);
 102       fasta_index = 1;
 103       annot_file = "";
 104       sub_seq_start = 0;
 105       sub_seq_end = 0;
 106       seq_params = true;
 107
 108       while ((!para_file.eof()) && seq_params)
 109       {
 110         getline(para_file,file_data_line);
 111         split_index = file_data_line.find(" ");
 112         param = file_data_line.substr(0,split_index);
 113         value = file_data_line.substr(split_index+1);
 114
 115         if (param == "FASTA_INDEX")
 116           fasta_index = atoi(value.c_str());
 117         else if (param == "ANNOTATION")
 118           annot_file = value;
 119         else if (param == "SEQ_START")
 120           sub_seq_start = atoi(value.c_str());
 121         else if (param == "SEQ_END")
 122         {
 123           cout << "hey!  " << atoi(value.c_str()) << endl;
 124           sub_seq_end = atoi(value.c_str());
 125         }
 126         //ignore empty lines or that start with '#'
 127         else if ((param == "") || (param == "#")) {}
 128         else seq_params = false;
 129       }
 130
 131
 132       fasta_indices.push_back(fasta_index);
 133       annot_files.push_back(annot_file);
 134       sub_seq_starts.push_back(sub_seq_start);
 135       sub_seq_ends.push_back(sub_seq_end);
 136       did_seq = true;
 137
 138     }
 139     //ignore empty lines or that start with '#'
 140     else if ((param == "") || (param == "#")) {}
 141     else
 142     {
 143       cout << "Illegal/misplaced mussa parameter in file\n";
 144       cout << param << "\n";
 145     }
 146
 147     if (!did_seq)
 148     {
 149       getline(para_file,file_data_line);
 150       split_index = file_data_line.find(" ");
 151       param = file_data_line.substr(0,split_index);
 152       value = file_data_line.substr(split_index+1);
 153       did_seq = false;
 154     }
 155   }
 156
 157   para_file.close();
 158
 159   cout << "ana_name = " << ana_name << win_append << win_append << "\n";
 160   cout << "window = " << window << " threshold = " << threshold << "\n";
 161 }
 162
 163 //        if (!((param == "") || (param == "#")))
 164 //          cout << value << " = " << param << endl;
 165 void
 166 Mussa::analyze(string para_file_path, int w, int t)
 167 {
 168   time_t t1, t2, begin, end;
 169   double setuptime, seqloadtime, seqcomptime, nwaytime, savetime, totaltime;
 170   int seq_num;
 171
 172   if (w > 0)
 173     window  = w;
 174   if (t > 0)
 175     threshold = t;
 176
 177   begin = time(NULL);
 178
 179   t1 = time(NULL);
 180   load_parameters(para_file_path);
 181   t2 = time(NULL);
 182   setuptime = difftime(t2, t1);
 183
 184
 185   cout << "fee\n";
 186   t1 = time(NULL);
 187   get_Seqs();
 188   t2 = time(NULL);
 189   seqloadtime = difftime(t2, t1);
 190
 191
 192   cout << "fie\n";
 193   t1 = time(NULL);
 194   seqcomp();
 195   t2 = time(NULL);
 196   seqcomptime = difftime(t2, t1);
 197
 198
 199   cout << "foe\n";
 200   t1 = time(NULL);
 201   nway();
 202   t2 = time(NULL);
 203   nwaytime = difftime(t2, t1);
 204
 205
 206   cout << "fum\n";
 207   t1 = time(NULL);
 208   save();
 209   t2 = time(NULL);
 210   savetime = difftime(t2, t1);
 211
 212   end = time(NULL);
 213   totaltime = difftime(end, begin);
 214
 215   cout << "setup\tseqload\tseqcomp\tnway\tsave\ttotal\n";
 216   cout << setuptime << "\t";
 217   cout << seqloadtime << "\t";
 218   cout << seqcomptime << "\t";
 219   cout << nwaytime << "\t";
 220   cout << savetime << "\t";
 221   cout << totaltime << "\n";
 222 }
 223
 224
 225 void
 226 Mussa::get_Seqs()
 227 {
 228   list<string>::iterator seq_files_i, annot_files_i;
 229   list<int>::iterator fasta_indices_i, seq_starts_i, seq_ends_i;
 230   Sequence aSeq;
 231
 232   seq_files_i = seq_files.begin();
 233   fasta_indices_i = fasta_indices.begin();
 234   annot_files_i = annot_files.begin();
 235   seq_starts_i = sub_seq_starts.begin();
 236   seq_ends_i = sub_seq_ends.begin();
 237
 238
 239   while (seq_files_i != seq_files.end())
 240           /* it should be guarenteed that each of the following exist
 241  &&
 242           (fasta_indices_i != fasta_indices.end()) &&
 243           (annot_files_i != annot_files.end())  &&
 244           (seq_starts_i != sub_seq_starts.end())  &&
 245           (seq_ends_i != sub_seq_ends.end())         )
 246           */
 247   {
 248     aSeq.load_fasta(*seq_files_i, *fasta_indices_i,*seq_starts_i,*seq_ends_i);
 249     if (*annot_files_i != "")
 250       aSeq.load_annot(*annot_files_i);
 251     the_Seqs.push_back(aSeq);
 252     cout << aSeq.hdr() << endl;
 253     //cout << aSeq.seq() << endl;
 254     aSeq.clear();
 255     ++seq_files_i;
 256     ++fasta_indices_i;
 257     ++annot_files_i;
 258     ++seq_starts_i;
 259     ++seq_ends_i;
 260   }
 261 }
 262
 263 /*
 264     aSeq.seq();
 265     cout << "\n";
 266     aSeq.hdr();
 267     cout << "\n";
 268 */
 269
 270 void
 271 Mussa::seqcomp()
 272 {
 273   int i, i2;     // loop vars over sequences to analyze
 274   vector<int> seq_lens;
 275   vector<FLPs> empty_FLP_vector;
 276   FLPs dummy_comp;
 277   string save_file_string;
 278
 279   empty_FLP_vector.clear();
 280   for(i = 0; i < seq_num; i++)
 281   {
 282     all_comps.push_back(empty_FLP_vector);
 283     for(i2 = 0; i2 < seq_num; i2++)
 284       all_comps[i].push_back(dummy_comp);
 285   }
 286
 287   for(i = 0; i < seq_num; i++)
 288     seq_lens.push_back(the_Seqs[i].len());
 289
 290   for(i = 0; i < seq_num; i++)
 291     for(i2 = i+1; i2 < seq_num; i2++)
 292     {
 293       all_comps[i][i2].setup("m", window, threshold, seq_lens[i],seq_lens[i2]);
 294       all_comps[i][i2].seqcomp(the_Seqs[i].seq(), the_Seqs[i2].seq(), false);
 295       all_comps[i][i2].seqcomp(the_Seqs[i].seq(),the_Seqs[i2].rev_comp(),true);
 296
 297       all_comps[i][i2].file_save(save_file_string);
 298     }
 299 }
 300
 301
 302
 303 void
 304 Mussa::nway()
 305 {
 306   the_paths.setup(seq_num, window, threshold);
 307   the_paths.find_paths_r(all_comps);
 308   the_paths.simple_refine();
 309 }
 310
 311 void
 312 Mussa::save()
 313 {
 314   string save_path_base, save_path;
 315   fstream save_file;
 316   int i;
 317
 318   // gotta do bit with adding win & thres if to be appended - need itos
 319
 320   // not sure why, but gotta close file each time since can't pass file streams
 321
 322   save_path_base = ana_name;
 323
 324   if (win_append)
 325     save_path_base += "_w" + int_to_str(window);
 326
 327   if (thres_append)
 328     save_path_base += "_t" + int_to_str(threshold);
 329
 330   // save sequence and annots to a special mussa file
 331   save_path = save_path_base + ".museq";
 332   save_file.open(save_path.c_str(), ios::out);
 333   save_file << "<Mussa_Sequence>" << endl;
 334   //save_file.close();
 335
 336   for(i = 0; i < seq_num; i++)
 337     the_Seqs[i].save(save_file);
 338
 339   //save_file.open(save_path.c_str(), ios::app);
 340   save_file << "</Mussa_Sequence>" << endl;
 341   save_file.close();
 342
 343   // save nway paths to its mussa save file
 344   save_path = save_path_base + ".muway";
 345   the_paths.save(save_path);
 346 }
 347
 348 void
 349 Mussa::load(string ana_file)
 350 {
 351   int i;
 352   string load_file_path;
 353   Sequence tmp_seq;
 354
 355
 356   ana_name = ana_file;
 357   load_file_path = ana_name + ".muway";
 358   seq_num = the_paths.load(load_file_path);
 359
 360   load_file_path = ana_name + ".museq";
 361   for (i = 1; i <= seq_num; i++)
 362   {
 363     tmp_seq.clear();
 364     tmp_seq.load_museq(load_file_path, i);
 365     the_Seqs.push_back(tmp_seq);
 366   }
 367 }
 368
 369 /*
 370       cout << "fee\n";
 371       cout << "fie\n";
 372       cout << "foe\n";
 373       cout << "fum\n";
 374 */
 375
 376
 377
 378 void
 379 Mussa::save_old()
 380 {
 381   fstream save_file;
 382   int i;
 383
 384   save_file.open(ana_name.c_str(), ios::out);
 385
 386   for(i = 0; i < seq_num; i++)
 387     save_file << the_Seqs[i].seq() << endl;
 388
 389   save_file << window << endl;
 390   save_file.close();
 391   //note more complex eventually since ana_name may need to have
 392   //window size, threshold and other stuff to modify it...
 393   the_paths.save_old(ana_name);
 394 }
 395
 396
 397 void
 398 Mussa::load_old(char * load_file_path, int s_num)
 399 {
 400   fstream save_file;
 401   string file_data_line;
 402   int i, space_split_i, comma_split_i;
 403   vector<int> loaded_path;
 404   string node_pair, node;
 405   Sequence a_seq;
 406
 407   seq_num = s_num;
 408   the_paths.setup(seq_num, 0, 0);
 409   save_file.open(load_file_path, ios::in);
 410
 411   // currently loads old mussa format
 412
 413   // get sequences
 414   for(i = 0; i < seq_num; i++)
 415   {
 416     getline(save_file, file_data_line);
 417     a_seq.set_seq(file_data_line);
 418     the_Seqs.push_back(a_seq);
 419   }
 420
 421   // get window size
 422   getline(save_file, file_data_line);
 423   window = atoi(file_data_line.c_str());
 424   // get paths
 425
 426   while (!save_file.eof())
 427   {
 428     loaded_path.clear();
 429     getline(save_file, file_data_line);
 430     if (file_data_line != "")
 431     for(i = 0; i < seq_num; i++)
 432     {
 433       space_split_i = file_data_line.find(" ");
 434       node_pair = file_data_line.substr(0,space_split_i);
 435       //cout << "np= " << node_pair;
 436       comma_split_i = node_pair.find(",");
 437       node = node_pair.substr(comma_split_i+1);
 438       //cout << "n= " << node << " ";
 439       loaded_path.push_back(atoi (node.c_str()));
 440       file_data_line = file_data_line.substr(space_split_i+1);
 441     }
 442     //cout << endl;
 443     the_paths.add_path(loaded_path);
 444   }
 445   save_file.close();
 446
 447   //the_paths.save("tmp.save");
 448 }