QProfile.h

   1 #ifndef QPROFILE_H_
   2 #define QPROFILE_H_
   3
   4 #include<cstdio>
   5 #include<cstring>
   6 #include<cassert>
   7
   8 #include "utils.h"
   9 #include "RefSeq.h"
  10 #include "simul.h"
  11
  12
  13 class QProfile {
  14 public:
  15         QProfile();
  16         QProfile& operator=(const QProfile&);
  17
  18         void init();
  19         void update(const std::string&, const std::string&, const RefSeq&, int, int, double);
  20         void finish();
  21
  22         double getProb(const std::string&, const std::string&, const RefSeq&, int, int);
  23
  24         void collect(const QProfile&);
  25
  26         void read(FILE*);
  27         void write(FILE*);
  28
  29         void startSimulation();
  30         std::string simulate(simul*, int, int, int, const std::string&, const RefSeq&);
  31         void finishSimulation();
  32
  33 private:
  34         static const int NCODES = 5; // number of possible codes
  35         static const int SIZE = 100;
  36
  37         double p[SIZE][NCODES][NCODES]; // p[q][r][c] = p(c|r,q)
  38
  39         //make sure that quality score in [0, 93]
  40         int c2q(char c) { assert(c >= 33 && c <= 126); return c - 33; }
  41
  42         double (*pc)[NCODES][NCODES]; // for simulation
  43 };
  44
  45 QProfile::QProfile() {
  46         memset(p, 0, sizeof(p));
  47
  48         //make initialized parameters
  49         //ASSUME order of A, C, G, T, N
  50         int N = NCODES - 1;
  51         double probN = 1e-5;
  52         double probC, probO; // current, other
  53
  54         for (int i = 0; i < SIZE; i++) {
  55                 for (int j = 0; j < NCODES - 1; j++) {
  56                         p[i][j][N] = probN;
  57
  58                         probO = exp(-i / 10.0 * log(10.0));
  59                         probC = 1.0 - probO;
  60                         probO /= (NCODES - 2);
  61
  62                         probC *= (1.0 - probN);
  63                         probO *= (1.0 - probN);
  64
  65                         assert(probC >= 0.0 && probO >= 0.0);
  66
  67                         for (int k = 0; k < NCODES - 1; k++) {
  68                                 if (j == k) p[i][j][k] = probC;
  69                                 else p[i][j][k] = probO;
  70                         }
  71                 }
  72                 p[i][N][N] = probN;
  73                 for (int k = 0; k < NCODES - 1; k++)
  74                         p[i][N][k] = (1.0 - probN) / (NCODES - 1);
  75         }
  76 }
  77
  78 QProfile& QProfile::operator=(const QProfile& rv) {
  79         if (this == &rv) return *this;
  80         memcpy(p, rv.p, sizeof(rv.p));
  81         return *this;
  82 }
  83
  84 void QProfile::init() {
  85         memset(p, 0, sizeof(p));
  86 }
  87
  88 void QProfile::update(const std::string& readseq, const std::string& qual, const RefSeq& refseq, int pos, int dir, double frac) {
  89         int len = readseq.size();
  90         for (int i = 0; i < len; i++) {
  91           p[c2q(qual[i])][refseq.get_id(i + pos, dir)][get_base_id(readseq[i])] += frac;
  92         }
  93 }
  94
  95 void QProfile::finish() {
  96         double sum;
  97
  98         for (int i = 0; i < SIZE; i++) {
  99                 for (int j = 0; j < NCODES; j++) {
 100                         sum = 0.0;
 101                         for (int k = 0; k < NCODES; k++) sum += p[i][j][k];
 102                         if (sum < EPSILON) {
 103                                 for (int k = 0; k < NCODES; k++) p[i][j][k] = 0.0;
 104                                 continue;
 105                         }
 106                         for (int k = 0; k < NCODES; k++) p[i][j][k] /= sum;
 107                 }
 108         }
 109 }
 110
 111 double QProfile::getProb(const std::string& readseq, const std::string& qual, const RefSeq& refseq, int pos, int dir) {
 112         double prob = 1.0;
 113         int len = readseq.size();
 114
 115         for (int i = 0; i < len; i++) {
 116                 prob *= p[c2q(qual[i])][refseq.get_id(i + pos, dir)][get_base_id(readseq[i])];
 117         }
 118
 119         return prob;
 120 }
 121
 122 void QProfile::collect(const QProfile& o) {
 123         for (int i = 0; i < SIZE; i++)
 124                 for (int j = 0; j < NCODES; j++)
 125                         for (int k = 0; k < NCODES; k++)
 126                                 p[i][j][k] += o.p[i][j][k];
 127 }
 128
 129 void QProfile::read(FILE *fi) {
 130         int tmp_size, tmp_ncodes;
 131         assert(fscanf(fi, "%d %d", &tmp_size, &tmp_ncodes) == 2);
 132         assert(tmp_size == SIZE && tmp_ncodes == NCODES);
 133         for (int i = 0; i < SIZE; i++)
 134                 for (int j = 0; j < NCODES; j++)
 135                         for (int k = 0; k < NCODES; k++)
 136                           assert(fscanf(fi, "%lf", &p[i][j][k]) == 1);
 137 }
 138
 139 void QProfile::write(FILE *fo) {
 140         fprintf(fo, "%d %d\n", SIZE, NCODES);
 141         for (int i = 0; i < SIZE; i++) {
 142                 for (int j = 0; j < NCODES; j++) {
 143                         for (int k = 0; k < NCODES - 1; k++)
 144                                 fprintf(fo, "%.10g ", p[i][j][k]);
 145                         fprintf(fo, "%.10g\n", p[i][j][NCODES - 1]);
 146                 }
 147                 if (i < SIZE - 1) { fprintf(fo, "\n"); }
 148         }
 149 }
 150
 151 void QProfile::startSimulation() {
 152         pc = new double[SIZE][NCODES][NCODES];
 153         for (int i = 0; i < SIZE; i++) {
 154           for (int j = 0; j < NCODES; j++)
 155             for (int k = 0; k < NCODES; k++) {
 156               pc[i][j][k] = p[i][j][k];
 157               if (k > 0) pc[i][j][k] += pc[i][j][k - 1];
 158             }
 159
 160           //avoid sampling from 0.0!!!
 161           double cp_sum, cp_d, cp_n;
 162           cp_sum = cp_d = cp_n = 0.0;
 163
 164           for (int j = 0; j < NCODES - 1; j++) {
 165             cp_sum += pc[i][j][NCODES - 1];
 166             cp_d += p[i][j][j];
 167             cp_n += p[i][j][NCODES - 1];
 168           }
 169
 170           if (cp_sum == 0.0) continue;
 171
 172           double p_d, p_o, p_n;
 173           p_d = cp_d / cp_sum;
 174           p_n = cp_n / cp_sum;
 175           p_o = (1.0 - p_d - p_n) / (NCODES - 2);
 176           for (int j = 0; j < NCODES - 1; j++) {
 177             if (pc[i][j][NCODES - 1] > 0.0) continue;
 178             for (int k = 0; k < NCODES; k++) {
 179               if (k == j) pc[i][j][k] = p_d;
 180               else if (k == NCODES - 1) pc[i][j][k] = p_n;
 181               else pc[i][j][k] = p_o;
 182               if (k > 0) pc[i][j][k] += pc[i][j][k - 1];
 183             }
 184           }
 185           if (pc[i][NCODES - 1][NCODES - 1] == 0.0) {
 186             p_o = (1.0 - p_n) / (NCODES - 1);
 187             for (int k = 0; k < NCODES; k++) {
 188               pc[i][NCODES - 1][k] = (k < NCODES - 1 ? p_o : p_n);
 189               if (k > 0) pc[i][NCODES - 1][k] += pc[i][NCODES - 1][k - 1];
 190             }
 191           }
 192         }
 193 }
 194
 195 std::string QProfile::simulate(simul* sampler, int len, int pos, int dir, const std::string& qual, const RefSeq& refseq) {
 196         std::string readseq = "";
 197
 198         for (int i = 0; i < len; i++) {
 199                 readseq.push_back(getCharacter(sampler->sample(pc[c2q(qual[i])][refseq.get_id(i + pos, dir)], NCODES)));
 200         }
 201         return readseq;
 202 }
 203
 204 void QProfile::finishSimulation() {
 205         delete[] pc;
 206 }
 207
 208 #endif /* QPROFILE_H_ */