NoiseQProfile.h

   1 #ifndef NOISEQPROFILE_H_
   2 #define NOISEQPROFILE_H_
   3
   4 #include<cmath>
   5 #include<cstdio>
   6 #include<cstring>
   7 #include<string>
   8
   9 #include "utils.h"
  10 #include "RefSeq.h"
  11 #include "simul.h"
  12
  13 class NoiseQProfile {
  14 public:
  15         NoiseQProfile() {
  16                 logp = 0.0;
  17                 memset(c, 0, sizeof(c));
  18                 memset(p, 0, sizeof(p));
  19         }
  20
  21         NoiseQProfile& operator=(const NoiseQProfile&);
  22
  23         void init();
  24         void updateC(const std::string&, const std::string&);
  25         void update(const std::string&, const std::string&, double frac);
  26         void finish();
  27         void calcInitParams();
  28
  29         double getProb(const std::string&, const std::string&);
  30         double getLogP() { return logp; }
  31
  32         void collect(const NoiseQProfile&);
  33
  34         void read(FILE*);
  35         void write(FILE*);
  36
  37         void startSimulation();
  38         std::string simulate(simul*, int, const std::string&);
  39         void finishSimulation();
  40
  41 private:
  42         static const int NCODES = 5; // number of possible codes
  43         static const int SIZE = 100;
  44
  45         double logp; //log prob;
  46         double c[SIZE][NCODES]; //counts in N0;
  47         double p[SIZE][NCODES]; //p[q][c] = p(c|q)
  48
  49         int c2q(char c) { assert(c >= 33 && c <= 126); return c - 33; }
  50
  51         double (*pc)[NCODES]; // for simulation
  52 };
  53
  54 NoiseQProfile& NoiseQProfile::operator=(const NoiseQProfile& rv) {
  55         if (this == &rv) return *this;
  56         logp = rv.logp;
  57         memcpy(c, rv.c, sizeof(rv.c));
  58         memcpy(p, rv.p, sizeof(rv.p));
  59         return *this;
  60 }
  61
  62 void NoiseQProfile::init() {
  63         memset(p, 0, sizeof(p));
  64 }
  65
  66 void NoiseQProfile::updateC(const std::string& readseq, const std::string& qual) {
  67         int len = readseq.size();
  68         for (int i = 0; i < len; i++) {
  69                 ++c[c2q(qual[i])][get_base_id(readseq[i])];
  70         }
  71 }
  72
  73 void NoiseQProfile::update(const std::string& readseq, const std::string& qual, double frac) {
  74         int len = readseq.size();
  75         for (int i = 0; i < len; i++) {
  76                 p[c2q(qual[i])][get_base_id(readseq[i])] += frac;
  77         }
  78 }
  79
  80 void NoiseQProfile::finish() {
  81         double sum;
  82
  83         //If N0 is 0, p(c|q) = 0 for all c, q
  84         logp = 0.0;
  85         for (int i = 0; i < SIZE; i++) {
  86                 sum = 0.0;
  87                 for (int j = 0; j < NCODES; j++) sum += (p[i][j] + c[i][j]);
  88                 if (sum <= 0.0) continue;
  89                 //if (isZero(sum)) continue;
  90                 for (int j = 0; j < NCODES; j++) {
  91                         p[i][j] = (p[i][j] + c[i][j]) /sum;
  92                         if (c[i][j] > 0.0) { logp += c[i][j] * log(p[i][j]); }
  93                 }
  94         }
  95 }
  96
  97 //make init parameters not zero
  98 void NoiseQProfile::calcInitParams() {
  99         double sum;
 100
 101         logp = 0.0;
 102         for (int i = 0; i < SIZE; i++) {
 103                 sum = 0.0;
 104                 for (int j = 0; j < NCODES; j++) sum += (1.0 + c[i][j]); // 1.0 pseudo count
 105                 for (int j = 0; j < NCODES; j++) {
 106                         p[i][j] = (c[i][j] + 1.0) / sum;
 107                         if (c[i][j] > 0.0) { logp += c[i][j] * log(p[i][j]); }
 108                 }
 109         }
 110 }
 111
 112 double NoiseQProfile::getProb(const std::string& readseq, const std::string& qual) {
 113         double prob = 1.0;
 114         int len = readseq.size();
 115
 116         for (int i = 0; i < len; i++) {
 117                 prob *= p[c2q(qual[i])][get_base_id(readseq[i])];
 118         }
 119
 120         return prob;
 121 }
 122
 123 void NoiseQProfile::collect(const NoiseQProfile& o) {
 124         for (int i = 0; i < SIZE; i++) {
 125                 for (int j = 0; j < NCODES; j++)
 126                         p[i][j] += o.p[i][j];
 127         }
 128 }
 129
 130 //If read from file, assume do not need to estimate from data
 131 void NoiseQProfile::read(FILE *fi) {
 132         int tmp_size, tmp_ncodes;
 133
 134         /* DEBUG */
 135         memset(c, 0, sizeof(c));
 136         //logp = 0.0;
 137         /* END */
 138
 139         fscanf(fi, "%d %d", &tmp_size, &tmp_ncodes);
 140         assert(tmp_size == SIZE && tmp_ncodes == NCODES);
 141         for (int i = 0; i < SIZE; i++) {
 142                 for (int j = 0; j < NCODES; j++) {
 143                         fscanf(fi, "%lf", &p[i][j]);
 144                         /* DEBUG */
 145                         //if (c[i][j] > 0) logp += c[i][j] * log(p[i][j]);
 146                         /* END */
 147                 }
 148         }
 149 }
 150
 151 void NoiseQProfile::write(FILE *fo) {
 152         fprintf(fo, "%d %d\n", SIZE, NCODES);
 153         for (int i = 0; i < SIZE; i++) {
 154                 for (int j = 0; j < NCODES - 1; j++) { fprintf(fo, "%.10g ", p[i][j]); }
 155                 fprintf(fo, "%.10g\n", p[i][NCODES - 1]);
 156         }
 157 }
 158
 159 void NoiseQProfile::startSimulation() {
 160         pc = new double[SIZE][NCODES];
 161
 162         for (int i = 0; i < SIZE; i++)
 163                 for (int j = 0; j < NCODES; j++) {
 164                         pc[i][j] = p[i][j];
 165                         if (j > 0) pc[i][j] += pc[i][j - 1];
 166                 }
 167 }
 168
 169 std::string NoiseQProfile::simulate(simul* sampler, int len, const std::string& qual) {
 170         std::string readseq = "";
 171
 172         for (int i = 0; i < len; i++) {
 173                 readseq.push_back(getCharacter(sampler->sample(pc[c2q(qual[i])], NCODES)));
 174         }
 175         return readseq;
 176 }
 177
 178 void NoiseQProfile::finishSimulation() {
 179         delete[] pc;
 180 }
 181
 182 #endif /* NOISEQPROFILE_H_ */