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