alg/flp.cpp

   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 //                            ---------- flp.cc  -----------
  13 //                        ----------------------------------------
  14
  15 #include <boost/filesystem/operations.hpp>
  16 #include <boost/filesystem/fstream.hpp>
  17 namespace fs = boost::filesystem;
  18
  19 #include "alg/flp.hpp"
  20
  21 #include <iostream>
  22 #include <string>
  23 #include <stdexcept>
  24 #include <cassert>
  25
  26 #include "mussa_exceptions.hpp"
  27 using namespace std;
  28
  29 bool operator==(const FLPs::match& a, const FLPs::match& b)
  30 {
  31   return (a.index == b.index && a.score == b.score);
  32 }
  33
  34 ostream &operator<<(ostream& out, const FLPs::match& m)
  35 {
  36   out << "(" << m.score << ", " << m.index << ")";
  37   return out;
  38 }
  39
  40 FLPs::FLPs() :
  41   window_size(0),
  42   hard_threshold(0),
  43   all_matches(0)
  44 {
  45 }
  46
  47 void
  48 FLPs::setup(int win_size, int hard_thres)
  49 {
  50   window_size = win_size;
  51   hard_threshold = hard_thres;
  52 }
  53
  54 void FLPs::alloc_matches(string::size_type len1)
  55 {
  56   if (all_matches != 0) {
  57     // if we're being called again its likely because
  58     // seqcomp is being called to do the reverse compliment
  59     // so we shouldn't reallocate ourselves.
  60     return;
  61   }
  62   all_matches = new std::vector<list<FLPs::match> >;
  63
  64   list<FLPs::match> empty_match_list;
  65   std::vector<list<FLPs::match> >::size_type window_count;
  66   if (len1 > 0) {
  67     // we actually know how much to preallocate
  68     window_count = len1 - window_size+1;
  69   } else {
  70     // we have no idea, so don't allocate anything
  71     window_count = 0;
  72   }
  73
  74   for (std::vector<list<FLPs::match> >::size_type i=0; i < window_count; ++i) {
  75     all_matches->push_back(empty_match_list);
  76   }
  77   assert (all_matches->size() == window_count);
  78 }
  79
  80
  81 int
  82 FLPs::size() const
  83 {
  84   if (all_matches != 0) {
  85     return all_matches->size();
  86   } else {
  87     return 0;
  88   }
  89 }
  90
  91
  92 /*
  93 bool
  94 FLPs::match_less(match *match1, match *match2)
  95 {
  96   if (match1->score < match2->score)
  97     return true;
  98   else if ( (match1->score == match2->score) &&
  99             (match1->index < match2->index) )
 100     return true;
 101   else
 102     return false;
 103 }
 104
 105 void
 106 FLPs::sort()
 107 {
 108   int i;
 109
 110   for(i = 0; i < seq1_win_num; i++)
 111     if (!all_matches[i].empty())
 112       all_matches[i].sort(&FLPs::match_less);
 113 }
 114 */
 115
 116 list<FLPs::match>
 117 FLPs::matches(int index) const
 118 {
 119   if (all_matches == 0)
 120     throw runtime_error("please call FLPs.seqcomp first");
 121
 122   list<FLPs::match>::iterator list_i, list_end;
 123
 124   index = abs(index);
 125   list_i = (*all_matches)[index].begin();
 126   list_end = (*all_matches)[index].end();
 127   list<FLPs::match> these_matches(list_i, list_end);
 128   return these_matches;
 129 }
 130
 131 list<int>
 132 FLPs::match_locations(int index) const
 133 {
 134   if (all_matches == 0)
 135     throw runtime_error("please call FLPs.seqcomp first");
 136
 137   list<int> these_matches;
 138   list<FLPs::match>::iterator list_i, list_end;
 139
 140   index = abs(index);
 141   list_i = (*all_matches)[index].begin();
 142   list_end = (*all_matches)[index].end();
 143   while (list_i != list_end)
 144   {
 145     these_matches.push_back(list_i->index);
 146     ++list_i;
 147   }
 148   return these_matches;
 149 }
 150
 151 list<int>
 152 FLPs::thres_matches(int index, int thres) const
 153 {
 154   if (all_matches == 0)
 155     throw runtime_error("please call FLPs.seqcomp first");
 156
 157   list<int> thres_matches;
 158   list<FLPs::match>::iterator list_i, list_end;
 159
 160   index = abs(index);
 161   list_i = (*all_matches)[index].begin();
 162   list_end = (*all_matches)[index].end();
 163   thres_matches.clear();
 164
 165   while (list_i != list_end)
 166   {
 167     if (list_i->score >= thres)
 168       thres_matches.push_back(list_i->index);
 169     ++list_i;
 170   }
 171   return thres_matches;
 172 }
 173
 174 void
 175 FLPs::save(fs::path save_file_path)
 176 {
 177   if (all_matches == 0)
 178     throw runtime_error("please call FLPs.seqcomp first");
 179
 180   fs::fstream save_file;
 181   save_file.open(save_file_path, ios::out);
 182
 183   save_file << "<Seqcomp win=" << window_size
 184             << " thres=" << hard_threshold << endl;
 185
 186   for(vector<list<FLPs::match> >::const_iterator win_citor = all_matches->begin();
 187       win_citor != all_matches->end();
 188       ++win_citor)
 189   {
 190     if (!win_citor->empty())
 191     {
 192       for( list<FLPs::match>::const_iterator match_itor = win_citor->begin();
 193            match_itor != win_citor->end();
 194            ++match_itor)
 195       {
 196         save_file << match_itor->index << "," << match_itor->score << " ";
 197       }
 198     }
 199     save_file << endl;
 200   }
 201
 202   save_file << "</Seqcomp>\n";
 203
 204   save_file.close();
 205 }
 206
 207 void
 208 FLPs::load(fs::path file_path)
 209 {
 210   fs::fstream data_file;
 211   string file_data, file_data_line, pair_data, index_data, score_data;
 212   match a_match;
 213   string::size_type split_index, comma_index;
 214   bool tag_open = false;
 215   list<FLPs::match> a_match_list;
 216   // initialize our all_matches pointer
 217   alloc_matches();
 218
 219   if (fs::exists(file_path)) {
 220     data_file.open(file_path, ios::in);
 221
 222     getline(data_file,file_data_line);
 223     // parse seqcomp open tag and parameters
 224     // eg <Seqcomp type=mussa win=30 thres=21>
 225     // if parse successful...
 226     tag_open = true;
 227
 228     while ((data_file.good()) && tag_open)
 229     {
 230       // intialize list to empty
 231       a_match_list.clear();
 232
 233       getline(data_file,file_data_line);
 234       if (file_data_line == "</Seqcomp>")
 235       {
 236         tag_open = false;
 237       }
 238       // parse line of matches
 239       else if (file_data_line == "")
 240       {
 241         //cout << "empty line\n";
 242         all_matches->push_back(a_match_list);
 243       }
 244       else
 245       {
 246         split_index = file_data_line.find(" ");
 247
 248         while (split_index != string::npos)
 249         {
 250           pair_data = file_data_line.substr(0,split_index);
 251           file_data_line = file_data_line.substr(split_index+1);
 252           //cout << "pair_data = " << pair_data << "...";
 253           // parse out the 2 pieces of data, index and score of pair match
 254           comma_index = pair_data.find(",");
 255           index_data = pair_data.substr(0, comma_index);
 256           a_match.index = atoi(index_data.c_str() );
 257           score_data = pair_data.substr(comma_index+1);
 258           a_match.score = atoi(score_data.c_str() );
 259           //cout << a_match.index << "," << a_match.score << " ";
 260
 261           a_match_list.push_back(a_match);
 262
 263           split_index = file_data_line.find(" ");
 264         }
 265         all_matches->push_back(a_match_list);
 266         //cout << all_matches->size() << "\n";
 267       }
 268     }
 269     //cout << "windows in flp = " << all_matches->size() << endl;
 270     data_file.close();
 271   } else {
 272     throw mussa_load_error(file_path.string() + "was not found");
 273   }
 274
 275 }
 276