Profile.h

   1 #ifndef PROFILE_H_
   2 #define PROFILE_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 Profile {
  14 public:
  15         Profile(int = 1000);
  16         ~Profile() {
  17                 delete[] p;
  18         }
  19
  20         Profile& operator=(const Profile&);
  21
  22         void init();
  23         void update(const std::string&, const RefSeq&, int, int, double);
  24         void finish();
  25
  26         double getProb(const std::string&, const RefSeq&, int, int);
  27
  28         void collect(const Profile&);
  29
  30         void read(FILE*);
  31         void write(FILE*);
  32
  33         void startSimulation();
  34         std::string simulate(simul*, int, int, int, const RefSeq&);
  35         void finishSimulation();
  36
  37 private:
  38         static const int NCODES = 5;
  39
  40         int proLen; // profile length
  41         int size; // # of items in p;
  42         double (*p)[NCODES][NCODES]; //profile matrices
  43
  44         double (*pc)[NCODES][NCODES]; // for simulation
  45 };
  46
  47 Profile::Profile(int maxL) {
  48         proLen = maxL;
  49         size = proLen * NCODES * NCODES;
  50         p = new double[proLen][NCODES][NCODES];
  51         memset(p, 0, sizeof(double) * size);
  52
  53         //set initial parameters
  54         int N = NCODES - 1;
  55         double probN = 1e-5, portionC = 0.99; //portionC, among ACGT, the portion of probability mass the correct base takes
  56         double probC, probO;
  57
  58         for (int i = 0; i < proLen; i++) {
  59                 for (int j = 0; j < NCODES - 1; j++) {
  60                         p[i][j][N] = probN;
  61                         probC = portionC * (1.0 - probN);
  62                         probO = (1.0 - portionC) / (NCODES - 2) * (1.0 - probN);
  63
  64                         for (int k = 0; k < NCODES - 1; k++) {
  65                                 p[i][j][k] = (j == k ? probC : probO);
  66                         }
  67                 }
  68                 p[i][N][N] = probN;
  69                 for (int k = 0; k < NCODES - 1; k++)
  70                                 p[i][N][k] = (1.0 - probN) / (NCODES - 1);
  71         }
  72 }
  73
  74 Profile& Profile::operator=(const Profile& rv) {
  75         if (this == &rv) return *this;
  76         if (proLen != rv.proLen) {
  77                 delete[] p;
  78                 proLen = rv.proLen;
  79                 size = rv.size;
  80                 p = new double[rv.proLen][NCODES][NCODES];
  81         }
  82         memcpy(p, rv.p, sizeof(double) * rv.size);
  83
  84         return *this;
  85 }
  86
  87 void Profile::init() {
  88         memset(p, 0, sizeof(double) * size);
  89 }
  90
  91 void Profile::update(const std::string& readseq, const RefSeq& refseq, int pos, int dir, double frac) {
  92         int len = readseq.size();
  93         for (int i = 0; i < len; i++) {
  94                 p[i][refseq.get_id(i + pos, dir)][get_base_id(readseq[i])] += frac;
  95         }
  96 }
  97
  98 void Profile::finish() {
  99         double sum;
 100
 101         for (int i = 0; i < proLen; i++) {
 102                 for (int j = 0; j < NCODES; j++) {
 103                         sum = 0.0;
 104                         for (int k = 0; k < NCODES; k++) sum += p[i][j][k];
 105                         if (sum < EPSILON) {
 106                                 for (int k = 0; k < NCODES; k++) p[i][j][k] = 0.0;
 107                                 continue;
 108                         }
 109                         for (int k = 0; k < NCODES; k++) p[i][j][k] /= sum;
 110                 }
 111         }
 112 }
 113
 114 double Profile::getProb(const std::string& readseq, const RefSeq& refseq, int pos, int dir) {
 115         double prob = 1.0;
 116         int len = readseq.size();
 117
 118         for (int i = 0; i < len; i++) {
 119                 prob *= p[i][refseq.get_id(i + pos, dir)][get_base_id(readseq[i])];
 120
 121         }
 122
 123         return prob;
 124 }
 125
 126 void Profile::collect(const Profile& o) {
 127         for (int i = 0; i < proLen; i++)
 128                 for (int j = 0; j < NCODES; j++)
 129                         for (int k = 0; k < NCODES; k++)
 130                                 p[i][j][k] += o.p[i][j][k];
 131 }
 132
 133 void Profile::read(FILE *fi) {
 134         int tmp_prolen, tmp_ncodes;
 135         assert(fscanf(fi, "%d %d", &tmp_prolen, &tmp_ncodes) == 2);
 136         assert(tmp_ncodes == NCODES);
 137         if (tmp_prolen != proLen) {
 138                 delete[] p;
 139                 proLen = tmp_prolen;
 140                 size = proLen * NCODES * NCODES;
 141                 p = new double[proLen][NCODES][NCODES];
 142                 memset(p, 0, sizeof(double) * size);
 143         }
 144
 145         for (int i = 0; i < proLen; i++)
 146                 for (int j = 0; j < NCODES; j++)
 147                         for (int k = 0; k < NCODES; k++)
 148                           assert(fscanf(fi, "%lf", &p[i][j][k]) == 1);
 149 }
 150
 151 void Profile::write(FILE* fo) {
 152         fprintf(fo, "%d %d\n", proLen, NCODES);
 153         for (int i = 0; i < proLen; i++) {
 154                 for (int j = 0; j < NCODES; j++) {
 155                         for (int k = 0; k < NCODES - 1; k++)
 156                                 fprintf(fo, "%.10g ", p[i][j][k]);
 157                         fprintf(fo, "%.10g\n", p[i][j][NCODES - 1]);
 158                 }
 159                 if (i < proLen - 1) { fprintf(fo, "\n"); }
 160         }
 161 }
 162
 163 void Profile::startSimulation() {
 164         pc = new double[proLen][NCODES][NCODES];
 165         for (int i = 0; i < proLen; i++) {
 166                 for (int j = 0; j < NCODES; j++)
 167                         for (int k = 0; k < NCODES; k++) {
 168                                 pc[i][j][k] = p[i][j][k];
 169                                 if (k > 0) pc[i][j][k] += pc[i][j][k - 1];
 170                         }
 171                 //avoid sampling from 0.0!!!
 172                 double cp_sum, cp_d, cp_n;
 173                 cp_sum = cp_d = cp_n = 0.0;
 174
 175                 for (int j = 0; j < NCODES - 1; j++) {
 176                   cp_sum += pc[i][j][NCODES - 1];
 177                   cp_d += p[i][j][j];
 178                   cp_n += p[i][j][NCODES - 1];
 179                 }
 180
 181                 if (cp_sum == 0.0) continue;
 182
 183                 double p_d, p_o, p_n;
 184                 p_d = cp_d / cp_sum;
 185                 p_n = cp_n / cp_sum;
 186                 p_o = (1.0 - p_d - p_n) / (NCODES - 2);
 187                 for (int j = 0; j < NCODES - 1; j++) {
 188                   if (pc[i][j][NCODES - 1] > 0.0) continue;
 189                   for (int k = 0; k < NCODES; k++) {
 190                     if (k == j) pc[i][j][k] = p_d;
 191                     else if (k == NCODES - 1) pc[i][j][k] = p_n;
 192                     else pc[i][j][k] = p_o;
 193                     if (k > 0) pc[i][j][k] += pc[i][j][k - 1];
 194                   }
 195                 }
 196                 if (pc[i][NCODES - 1][NCODES - 1] == 0.0) {
 197                   p_o = (1.0 - p_n) / (NCODES - 1);
 198                   for (int k = 0; k < NCODES; k++) {
 199                     pc[i][NCODES - 1][k] = (k < NCODES - 1 ? p_o : p_n);
 200                     if (k > 0) pc[i][NCODES - 1][k] += pc[i][NCODES - 1][k - 1];
 201                   }
 202                 }
 203         }
 204
 205 }
 206
 207 std::string Profile::simulate(simul* sampler, int len, int pos, int dir, const RefSeq& refseq) {
 208         std::string readseq = "";
 209
 210         for (int i = 0; i < len; i++) {
 211                 readseq.push_back(getCharacter(sampler->sample(pc[i][refseq.get_id(i + pos, dir)], NCODES)));
 212         }
 213         return readseq;
 214 }
 215
 216 void Profile::finishSimulation() {
 217         delete[] pc;
 218 }
 219
 220 #endif /* PROFILE_H_ */