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