mussa_nway_entropy.cc

   1 /*
   2 #include <list>
   3 #include <vector>
   4 #include <string>
   5 #include <stdlib.h>
   6 #include <algorithm>
   7
   8 #include <iostream>
   9
  10 using namespace std;
  11 */
  12 #include "mussa_nway.hh"
  13 #include <math.h>
  14
  15 // vars that will be availble to entropy function when in its nway class
  16 //vector<char *> c_sequences;
  17 //int window = 5, seq_num;
  18
  19
  20 void
  21 Nway_Paths::setup_ent(double new_entropy_thres, vector<string> some_Seqs)
  22 {
  23   int i;
  24
  25   ent_thres = new_entropy_thres;
  26
  27   c_sequences.clear();
  28   for (i = 0; i < species_num; i++)
  29     c_sequences.push_back((char *)some_Seqs[i].c_str());
  30 }
  31
  32
  33 double
  34 Nway_Paths::path_entropy(vector<int> path)
  35 {
  36   int bp_occurences[5];
  37   double bp_prob, frac_ent, avg_entropy;
  38   vector<double> entropies;
  39   int i, seq_i, i2;
  40   int window = win_size;
  41   int seq_num = species_num;
  42
  43
  44   for(i = 0; i < window; i++)     // loop thru all bp positions
  45   {
  46     for(i2 = 0; i2 < 5; i2++)     // clear old occurences
  47       bp_occurences[i2] = 0;
  48
  49     for(seq_i = 0; seq_i < seq_num; seq_i++) //loop thru all sequences at pos i
  50       if (path[seq_i] > -1)
  51         switch (c_sequences[seq_i][i+path[seq_i]])
  52         {
  53           case 'A':
  54             ++bp_occurences[0];
  55             break;
  56           case 'T':
  57             ++bp_occurences[1];
  58             break;
  59           case 'G':
  60             ++bp_occurences[2];
  61             break;
  62           case 'C':
  63             ++bp_occurences[3];
  64             break;
  65           case 'N':
  66             ++bp_occurences[4];
  67             break;
  68           default:
  69             cout << "error, nonAGCTN on seq: " << seq_i << " at index: ";
  70             cout << i+path[seq_i] << endl;
  71         }
  72       else
  73         switch (c_sequences[seq_i][-1 * (i + path[seq_i] - window + 1) ])
  74         {
  75           case 'A':
  76             ++bp_occurences[1];
  77             break;
  78           case 'T':
  79             ++bp_occurences[0];
  80             break;
  81           case 'G':
  82             ++bp_occurences[3];
  83             break;
  84           case 'C':
  85             ++bp_occurences[2];
  86             break;
  87           case 'N':
  88             ++bp_occurences[4];
  89             break;
  90           default:
  91             cout << "error, nonAGCTN on seq: " << seq_i << " at index: ";
  92             cout << (i + path[seq_i] - window + 1) << endl;
  93         }
  94
  95
  96     entropies.push_back(0.000);
  97     for(i2 = 0; i2 < 5; i2++)    //calculate entropies
  98     {
  99       bp_prob = (double) bp_occurences[i2] / (double) seq_num;
 100       //cout << bp_prob << "::";
 101       if (bp_prob == 0.0000)
 102         frac_ent = 0.0000;
 103       else
 104         frac_ent = bp_prob * ( log10(bp_prob) / log10(2.0) );
 105       //cout << frac_ent << "  ";
 106       entropies[i] += frac_ent;
 107     }
 108     //cout << endl;
 109   }
 110
 111   avg_entropy = 0;
 112   for(i = 0; i < window; i++)
 113   {
 114     avg_entropy += entropies[i];
 115     //cout  << entropies[i] << endl;
 116   }
 117
 118   avg_entropy = -1.00 * (avg_entropy / (double) window);
 119   //cout << "average entropy: " << avg_entropy << endl;
 120   return avg_entropy;
 121 }
 122
 123
 124
 125
 126
 127 void
 128 Nway_Paths::entropy_path_search(vector<vector<FLPs> > all_comparisons)
 129 {
 130   vector<int> path;
 131   double avg_entropy;
 132   int win_i, sp_i, sp_depth, window_num, i, cur_node;
 133   bool some_matches, still_paths, not_advanced, trans_good;
 134   list<int> new_nodes, trans_check_nodes;
 135   vector<list<int> > all_matches;
 136   vector<list<int>::iterator> sp_nodes, sp_nodes_end;
 137   list<int>::iterator debug_i;
 138
 139
 140   pathz.clear();
 141   window_num = all_comparisons[0][1].win_num();
 142   cout << window_num << endl;    // just a sanity check
 143
 144   sp_nodes.reserve(species_num);
 145   sp_nodes_end.reserve(species_num);
 146
 147   // loop thru all windows in first species
 148   for (win_i = 0; win_i < window_num; win_i++)
 149   {
 150     // first we see if the first seq has matches to all other seqs at this win
 151     some_matches = true;
 152     sp_i = 1;
 153     all_matches.clear();
 154     while ( (sp_i < species_num) && (some_matches) )
 155     {
 156       new_nodes.clear();
 157       new_nodes = all_comparisons[0][sp_i].matches(win_i);
 158       if (new_nodes.empty())
 159         some_matches = false;
 160
 161       all_matches.push_back(new_nodes);
 162       sp_i++;
 163     }
 164     //cout << "fee\n";
 165
 166     // if 1st seq does match to all others, make all possible paths
 167     // out of all these matches
 168     if (some_matches)
 169     {
 170       sp_nodes.clear();
 171       sp_nodes_end.clear();
 172       // set each species list of matches to beginning
 173       for (sp_i = 0; sp_i < species_num-1; sp_i++)
 174       {
 175         sp_nodes.push_back(all_matches[sp_i].begin());
 176         sp_nodes_end.push_back(all_matches[sp_i].end());
 177       }
 178
 179       still_paths = true;
 180       //cout << "fie\n";
 181       while (still_paths)
 182       {
 183         // add path that each species iterator is pointing to
 184         path.clear();
 185         path.push_back(win_i);
 186
 187         //cout << win_i;
 188         for (sp_i = 0; sp_i < species_num-1; sp_i++)
 189         {
 190           //cout  << ", " << *(sp_nodes[sp_i]);
 191           path.push_back(*(sp_nodes[sp_i]));
 192         }
 193         //cout << endl;
 194
 195         // check entropy <---------------------------------------------------
 196         avg_entropy = path_entropy(path);
 197         if (avg_entropy <= ent_thres)
 198           pathz.push_back(path);
 199
 200         // now advance the right iterator
 201         not_advanced = true;
 202         sp_depth = species_num - 2;       // this roves up & down species list
 203         while ( (not_advanced) && (sp_depth != -1) )
 204         {
 205           //cout << "foe\n";
 206           //cout << sp_depth << ".." << *(sp_nodes[sp_depth]) << endl;
 207           (sp_nodes[sp_depth])++;   // advance iter at current depth
 208           //cout << sp_depth << ".." << *(sp_nodes[sp_depth]) << endl;
 209
 210           // if we reached the end of the list, reset iter & go up one depth
 211           if (sp_nodes[sp_depth] == sp_nodes_end[sp_depth])
 212           {
 213             //cout << "fum\n";
 214             sp_nodes[sp_depth] = all_matches[sp_depth].begin();
 215             sp_depth--;
 216             //cout << "depth = " << sp_depth << endl;
 217           }
 218           else
 219             not_advanced = false;
 220         }
 221
 222         if (sp_depth == -1)   // jumped up to first species, all paths searched
 223           still_paths = false;
 224         else                 // otherwise just reset to max depth and continue
 225           sp_depth = species_num - 2;
 226       }
 227     }
 228   }
 229 }
 230
 231
 232 //initial coding testing
 233 /*
 234 int main(int argc, char **argv)
 235 {
 236   vector<string> sequences;
 237   vector<int> i_starts;
 238   double avg_entropy;
 239   int seq_i;
 240
 241   sequences.clear();
 242   sequences.push_back("AAAAA");
 243   sequences.push_back("AAAAT");
 244   sequences.push_back("AATTG");
 245   sequences.push_back("ATTGC");
 246
 247   seq_num = 4;
 248
 249   c_sequences.clear();
 250   i_starts.clear();
 251   for(seq_i = 0; seq_i < seq_num; seq_i++)   //loop thru all sequences at pos i
 252   {
 253     c_sequences.push_back((char *)sequences[seq_i].c_str());
 254     i_starts.push_back(0);
 255   }
 256
 257   avg_entropy = Window_Entropy(i_starts);
 258
 259   cout << "average entropy: " << avg_entropy << endl;
 260 }
 261 */