mussa_nway_entropy.cc

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