mussa_nway.cc

   1 //                        ----------------------------------------
   2 //                         ---------- mussa_nway.cc  -----------
   3 //                        ----------------------------------------
   4
   5 #include "mussa_nway.hh"
   6 /*
   7       cout << "fee\n";
   8       cout << "fie\n";
   9       cout << "foe\n";
  10       cout << "fum\n";
  11 */
  12
  13
  14 Nway_Paths::Nway_Paths()
  15 {
  16 }
  17
  18 void
  19 Nway_Paths::setup(int sp_num, int w, int t)
  20 {
  21   species_num = sp_num;
  22   threshold = t;
  23   win_size = w;
  24   pathz.clear();
  25 }
  26
  27
  28 void
  29 Nway_Paths::path_search(vector<vector<FLPs> > all_comparisons, vector<int> path, int depth)
  30 {
  31   list<int> new_nodes, trans_check_nodes;
  32   list<int>::iterator new_nodes_i, new_nodes_end;
  33   int i;
  34   bool trans_check_good;
  35
  36   new_nodes = all_comparisons[depth - 1][depth].matches(path[depth-1]);
  37   new_nodes_i = new_nodes.begin();
  38   new_nodes_end = new_nodes.end();
  39   while(new_nodes_i != new_nodes_end)
  40   {
  41     //cout << "    * species " << depth << " node: " << *new_nodes_i << endl;
  42     // check transitivity with previous nodes in path
  43     for(i = 0; i < depth - 1; i++)
  44     {
  45       trans_check_good = true;
  46       trans_check_nodes = all_comparisons[i][depth].matches(path[i]);
  47       if ( (trans_check_nodes.end() == find(trans_check_nodes.begin(),
  48                                             trans_check_nodes.end(),
  49                                             *new_nodes_i) ) &&
  50            (trans_check_nodes.end() == find(trans_check_nodes.begin(),
  51                                            trans_check_nodes.end(),
  52                                            *new_nodes_i * -1) ) )
  53         trans_check_good = false;
  54     }
  55
  56     if (trans_check_good)
  57     {
  58       // this makes sure path nodes are recorded with RC status relative to
  59       // the base species
  60       if ( path[depth-1] >= 0)
  61         path.push_back(*new_nodes_i);
  62       else
  63         path.push_back(*new_nodes_i * -1);
  64
  65       if (depth < species_num - 1)
  66         path_search(all_comparisons, path, depth + 1);
  67       else
  68         pathz.push_back(path);
  69       path.pop_back();
  70     }
  71     ++new_nodes_i;
  72   }
  73 }
  74 //          cout << "    ****I have the power...\n";
  75
  76 /* use this if I ever get the friggin seqcomp match lists to sort...
  77       if (binary_search(trans_check_nodes.begin(), trans_check_nodes.end(),
  78                         *new_nodes_i))
  79 */
  80
  81 void
  82 Nway_Paths::find_paths_r(vector<vector<FLPs> > all_comparisons)
  83 {
  84   vector<int> path;
  85   int win_i, window_num;
  86   list<int> new_nodes;
  87   list<int>::iterator new_nodes_i, new_nodes_end;
  88
  89   pathz.clear();
  90   window_num = all_comparisons[0][1].win_num();
  91   cout << window_num << endl;
  92   // loop thru all windows in first species
  93   for (win_i = 0; win_i < window_num; win_i++)
  94   {
  95     path.clear();
  96     path.push_back(win_i);
  97     new_nodes = all_comparisons[0][1].matches(path[0]);
  98     new_nodes_i = new_nodes.begin();
  99     new_nodes_end = new_nodes.end();
 100     //if (new_nodes_i != new_nodes_end)
 101     //cout << "* species 0 node: " << win_i << endl;
 102     //  cout << "foookin hell\n";
 103     path.push_back(0);
 104     while(new_nodes_i != new_nodes_end)
 105     {
 106       //cout << "  * species 1 node: " << *new_nodes_i << endl;
 107       path[1] = *new_nodes_i;
 108       path_search(all_comparisons, path, 2);
 109       ++new_nodes_i;
 110     }
 111   }
 112 }
 113
 114 // dumbly goes thru and combines path windows exactly adjacent (ie + 1 index)
 115 // doesn't deal with interleaved adjacency
 116 void
 117 Nway_Paths::simple_refine()
 118 {
 119   // ext_path remembers the first window set in an extending path
 120   vector<int> ext_path, cur_path, next_path, new_path;
 121   list<vector<int> >::iterator pathz_i;
 122   int i2, win_ext_len = 0;
 123   bool extending;
 124
 125
 126   refined_pathz.clear();
 127
 128   pathz_i = pathz.begin();
 129   ext_path = *pathz_i;
 130   while(pathz_i != pathz.end())
 131   {
 132     // keep track of current path and advance to next path
 133     cur_path = *pathz_i;
 134     ++pathz_i;
 135     next_path = *pathz_i;
 136
 137     // if node for each seq is equal to the next node+1 then for all
 138     // sequences then we are extending
 139     extending = true;
 140     for(i2 = 0; i2 < species_num; i2++)
 141     {
 142       //cout << cur_path[i2] << " vs " << next_path[i2] << endl;
 143       if (cur_path[i2] + 1 != next_path[i2])
 144         extending = false;
 145     }
 146     if (extending)
 147     {
 148       //cout << "Raaaawwwrrr!  I am extending\n";
 149       win_ext_len++;
 150     }
 151     else
 152     {
 153       // add the extend window length as first element and add as refined
 154       new_path = ext_path;
 155       new_path.insert(new_path.begin(),win_size + win_ext_len);
 156       refined_pathz.push_back(new_path);
 157       // reset stuff
 158       win_ext_len = 0;
 159       ext_path = next_path;
 160     }
 161   }
 162 }
 163
 164
 165 void
 166 Nway_Paths::add_path(vector<int> loaded_path)
 167 {
 168   pathz.push_back(loaded_path);
 169 }
 170
 171
 172 void
 173 Nway_Paths::save(string save_file_path)
 174 {
 175   fstream save_file;
 176   list<vector<int> >::iterator path_i, paths_end;
 177   vector<int> a_path;
 178   int i;
 179
 180   save_file.open(save_file_path.c_str(), ios::out);
 181
 182   save_file << "<Mussa type=flp seq_num=" << species_num;
 183   save_file << " win=" << win_size;
 184   // add a function para new_thres defaults to -1 to later deal with
 185   // reanalysis with higher thres - if statement whether to record base thres
 186   // or new thres (ie if -1, then base)
 187   save_file << " thres=" << threshold << " >\n";
 188
 189   path_i = refined_pathz.begin();
 190   paths_end = refined_pathz.end();
 191   //path_i = pathz.begin();
 192   //paths_end = pathz.end();
 193   while (path_i != paths_end)
 194   {
 195     a_path = *path_i;
 196     //cout << a_path.size() << endl;
 197     //first entry is the window length of the windows in the path
 198     save_file << a_path[0] << ":";
 199     for(i = 1; i <= species_num; ++i)
 200     {
 201       save_file << a_path[i];
 202       if (i != species_num)
 203         save_file << ",";
 204     }
 205     save_file << endl;
 206     ++path_i;
 207   }
 208
 209   save_file << "</Mussa>\n";
 210   save_file.close();
 211 }
 212
 213 /*
 214   if (path_i == paths_end)
 215     cout << "Arrrrrrgggghhhhhh\n";
 216 */
 217
 218 int
 219 Nway_Paths::load(string load_file_path)
 220 {
 221   fstream load_file;
 222   string file_data_line, header_data, data, path_node, path_length;
 223   int i, space_split_i, equal_split_i, comma_split_i, colon_split_i;
 224   vector<int> loaded_path;
 225
 226
 227   load_file.open(load_file_path.c_str(), ios::in);
 228
 229   // get header data
 230   // grab mussa tag - discard for now...maybe check in future...
 231   getline(load_file,file_data_line);
 232   space_split_i = file_data_line.find(" ");
 233   file_data_line = file_data_line.substr(space_split_i+1);
 234   // grab type tag - need in future to distinguish between flp and vlp paths
 235   space_split_i = file_data_line.find(" ");
 236   file_data_line = file_data_line.substr(space_split_i+1);
 237   // get species/seq number
 238   space_split_i = file_data_line.find(" ");
 239   header_data = file_data_line.substr(0,space_split_i);
 240   equal_split_i = header_data.find("=");
 241   data = file_data_line.substr(equal_split_i+1);
 242   species_num = atoi (data.c_str());
 243   file_data_line = file_data_line.substr(space_split_i+1);
 244   // get window size
 245   space_split_i = file_data_line.find(" ");
 246   header_data = file_data_line.substr(0,space_split_i);
 247   equal_split_i = header_data.find("=");
 248   data = file_data_line.substr(equal_split_i+1);
 249   win_size = atoi (data.c_str());
 250   file_data_line = file_data_line.substr(space_split_i+1);
 251   // get window size
 252   space_split_i = file_data_line.find(" ");
 253   header_data = file_data_line.substr(0,space_split_i);
 254   equal_split_i = header_data.find("=");
 255   data = file_data_line.substr(equal_split_i+1);
 256   threshold = atoi (data.c_str());
 257   file_data_line = file_data_line.substr(space_split_i+1);
 258
 259   cout << "seq_num=" << species_num << " win=" << win_size;
 260   cout << " thres=" << threshold << endl;
 261
 262   // clear out the current data
 263   refined_pathz.clear();
 264
 265   int temp;
 266
 267   getline(load_file,file_data_line);
 268   while ( (!load_file.eof()) && (file_data_line != "</Mussa>") )
 269   {
 270     if (file_data_line != "")
 271     {
 272       loaded_path.clear();
 273       colon_split_i = file_data_line.find(":");
 274       path_length = file_data_line.substr(0,colon_split_i);
 275       loaded_path.push_back(atoi (path_length.c_str()));
 276       file_data_line = file_data_line.substr(colon_split_i+1);
 277       for(i = 0; i < species_num; i++)
 278       {
 279         comma_split_i = file_data_line.find(",");
 280         path_node = file_data_line.substr(0, comma_split_i);
 281         temp = atoi (path_node.c_str());
 282         loaded_path.push_back(temp);
 283         file_data_line = file_data_line.substr(comma_split_i+1);
 284       }
 285       refined_pathz.push_back(loaded_path);
 286     }
 287     getline(load_file,file_data_line);
 288   }
 289
 290
 291   load_file.close();
 292
 293   return species_num;
 294 }
 295
 296
 297
 298 void
 299 Nway_Paths::save_old(string save_file_path)
 300 {
 301   fstream save_file;
 302   list<vector<int> >::iterator path_i, paths_end;
 303   vector<int> a_path;
 304   int i;
 305
 306   save_file.open(save_file_path.c_str(), ios::app);
 307
 308   path_i = pathz.begin();
 309   paths_end = pathz.end();
 310   while(path_i != paths_end)
 311   {
 312     a_path = *path_i;
 313     //cout << a_path.size() << endl;
 314     for(i = 0; i < species_num; ++i)
 315       save_file << i << "," << a_path[i] << " ";
 316     save_file << endl;
 317     ++path_i;
 318   }
 319   save_file.close();
 320 }
 321
 322
 323 /*
 324 void
 325 Nway_Paths::find_paths(vector<vector<FLPs> > all_comparisons)
 326 {
 327   int win_i, sp_i;
 328   vector<list<list<int> > > path_src_tree;
 329   list<int> new_nodes;
 330   list<int>::iterator node_i, node_end;
 331   <list<list<int> > >::iterator branch_i, branch_end;
 332
 333   pathz.clear();
 334   path_src_tree.reserve(species_num - 1);
 335
 336
 337   // loop thru all windows in first species
 338   for (win_i = 0; win_i < window_num; win_i++)
 339   {
 340     // clear the path search tree
 341     for(i = 0; i < species_num; i++)
 342       path_src_tree[i].clear();
 343
 344     // top level kept empty even tho implicity has one entry of the first
 345     // species at this window - why bother, adds a little speed
 346
 347     // get connection list for first species, creating a list of nodes
 348     // of second species connected to the first species at this window
 349     new_nodes = all_comparisons[0][1];
 350     path_src_tree[1].push_back(new_nodes);
 351
 352     // loop thru rest of species for this window to see if any paths of matches
 353     // go across all species
 354     // if path search tree becomes empty, break out of loop, no reason to search further
 355     sp_i = 1;
 356     while ((sp_i < species_num) && (path tree not empty))
 357     {
 358       branch_i = path_src_tree[1].begin();
 359       branch_end = path_src_tree[1].end();
 360       while (branch_i != branch_end)
 361       {
 362         node_i = branch_i->begin();
 363         node_end = branch_i->end();
 364       }
 365
 366
 367       // loop over all current nodes
 368          // get connection list for each node
 369          // loop over each previous node in list
 370             // get those nodes connection list
 371             // intersect previous node connections with current
 372
 373       ++sp_i;
 374     }
 375
 376     // insert any of the paths found into the master list of paths
 377
 378     // add no paths if tmp_pathz is empty...
 379   }
 380 }
 381
 382 void Nway_Paths::refine()
 383 {
 384 }
 385
 386 void Nway_Paths::print()
 387 {
 388   list<list<int> >::iterator pathz_i;
 389
 390   cout << "printing list of lists\n";
 391   for (pathz_i = pathz.begin(); pathz_i != pathz.end(); ++pathz_i)
 392   {
 393     for (path_i = pathz_i->begin(); path_i != pathz_i->end(); ++path_i)
 394       cout << *path_i << " ";
 395     cout << endl;
 396   }
 397 }
 398 */
 399
 400
 401
 402