alg/nway_paths.cxx

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