SingleQModel.h

   1 #ifndef SINGLEQMODEL_H_
   2 #define SINGLEQMODEL_H_
   3
   4 #include<cmath>
   5 #include<cstdio>
   6 #include<cassert>
   7 #include<cstring>
   8 #include<string>
   9 #include<algorithm>
  10 #include<sstream>
  11
  12 #include "utils.h"
  13 #include "Orientation.h"
  14 #include "LenDist.h"
  15 #include "RSPD.h"
  16 #include "QualDist.h"
  17 #include "QProfile.h"
  18 #include "NoiseQProfile.h"
  19
  20 #include "ModelParams.h"
  21 #include "RefSeq.h"
  22 #include "Refs.h"
  23 #include "SingleReadQ.h"
  24 #include "SingleHit.h"
  25 #include "ReadReader.h"
  26
  27 #include "simul.h"
  28
  29 class SingleQModel {
  30 public:
  31         SingleQModel(Refs* refs = NULL) {
  32                 this->refs = refs;
  33                 M = (refs != NULL ? refs->getM() : 0);
  34                 memset(N, 0, sizeof(N));
  35                 estRSPD = false;
  36                 needCalcConPrb = true;
  37
  38                 ori = new Orientation();
  39                 gld = new LenDist();
  40                 mld = NULL;
  41                 rspd = new RSPD(estRSPD);
  42                 qd = new QualDist();
  43                 qpro = new QProfile();
  44                 nqpro = new NoiseQProfile();
  45
  46                 mean = -1.0; sd = 0.0;
  47                 mw = NULL;
  48
  49                 seedLen = 0;
  50         }
  51
  52         //If it is not a master node, only init & update can be used!
  53         SingleQModel(ModelParams& params, bool isMaster = true) {
  54                 M = params.M;
  55                 memcpy(N, params.N, sizeof(params.N));
  56                 refs = params.refs;
  57                 estRSPD = params.estRSPD;
  58                 mean = params.mean; sd = params.sd;
  59                 seedLen = params.seedLen;
  60                 needCalcConPrb = true;
  61
  62                 ori = NULL; gld = NULL; mld = NULL; rspd = NULL; qd = NULL; qpro = NULL; nqpro = NULL;
  63                 mw = NULL;
  64
  65                 if (isMaster) {
  66                         ori = new Orientation(params.probF);
  67                         gld = new LenDist(params.minL, params.maxL);
  68                         if (mean >= EPSILON) {
  69                                 mld = new LenDist(params.mate_minL, params.mate_maxL);
  70                         }
  71                         if (!estRSPD) { rspd = new RSPD(estRSPD); }
  72                         qd = new QualDist();
  73                 }
  74
  75                 if (estRSPD) { rspd = new RSPD(estRSPD, params.B); }
  76                 qpro = new QProfile();
  77                 nqpro = new NoiseQProfile();
  78         }
  79
  80         ~SingleQModel() {
  81                 refs = NULL;
  82                 if (ori != NULL) delete ori;
  83                 if (gld != NULL) delete gld;
  84                 if (mld != NULL) delete mld;
  85                 if (rspd != NULL) delete rspd;
  86                 if (qd != NULL) delete qd;
  87                 if (qpro != NULL) delete qpro;
  88                 if (nqpro != NULL) delete nqpro;
  89                 if (mw != NULL) delete[] mw;
  90                 /*delete[] p1, p2;*/
  91         }
  92
  93         //SingleQModel& operator=(const SingleQModel&);
  94
  95         void estimateFromReads(const char*);
  96
  97         //if prob is too small, just make it 0
  98         double getConPrb(const SingleReadQ& read, const SingleHit& hit) const {
  99                 if (read.isLowQuality()) return 0.0;
 100
 101                 double prob;
 102                 int sid = hit.getSid();
 103                 RefSeq &ref = refs->getRef(sid);
 104                 int fullLen = ref.getFullLen();
 105                 int totLen = ref.getTotLen();
 106                 int dir = hit.getDir();
 107                 int pos = hit.getPos();
 108                 int readLen = read.getReadLength();
 109                 int fpos = (dir == 0 ? pos : totLen - pos - readLen); // the aligned position reported in SAM file, should be a coordinate in forward strand
 110
 111                 assert(fpos >= 0 && fpos + readLen <= totLen && readLen <= totLen);
 112                 int seedPos = (dir == 0 ? pos : totLen - pos - seedLen); // the aligned position of the seed in forward strand coordinates
 113                 if (seedPos >= fullLen || ref.getMask(seedPos)) return 0.0;
 114
 115                 int effL;
 116                 double value;
 117
 118                 if (mld != NULL) {
 119                         int minL = std::max(readLen, gld->getMinL());
 120                         int maxL = std::min(totLen - pos, gld->getMaxL());
 121                         int pfpos; // possible fpos for fragment
 122                         value = 0.0;
 123                         for (int fragLen = minL; fragLen <= maxL; fragLen++) {
 124                                 pfpos = (dir == 0 ? pos : totLen - pos - fragLen);
 125                                 effL = std::min(fullLen, totLen - fragLen + 1);
 126                                 value += gld->getAdjustedProb(fragLen, totLen) * rspd->getAdjustedProb(pfpos, effL, fullLen) * mld->getAdjustedProb(readLen, fragLen);
 127                         }
 128                 }
 129                 else {
 130                         effL = std::min(fullLen, totLen - readLen + 1);
 131                         value = gld->getAdjustedProb(readLen, totLen) * rspd->getAdjustedProb(fpos, effL, fullLen);
 132                 }
 133
 134                 prob = ori->getProb(dir) * value * qpro->getProb(read.getReadSeq(), read.getQScore(), ref, pos, dir);
 135
 136                 if (prob < EPSILON) { prob = 0.0; }
 137
 138                 prob = (mw[sid] < EPSILON ? 0.0 : prob / mw[sid]);
 139
 140                 return prob;
 141         }
 142
 143         double getNoiseConPrb(const SingleReadQ& read) {
 144                 if (read.isLowQuality()) return 0.0;
 145                 double prob = mld != NULL ? mld->getProb(read.getReadLength()) : gld->getProb(read.getReadLength());
 146                 prob *= nqpro->getProb(read.getReadSeq(), read.getQScore());
 147                 if (prob < EPSILON) { prob = 0.0; }
 148
 149                 prob = (mw[0] < EPSILON ? 0.0 : prob / mw[0]);
 150
 151                 return prob;
 152         }
 153
 154         double getLogP() { return nqpro->getLogP(); }
 155
 156         void init();
 157
 158         void update(const SingleReadQ& read, const SingleHit& hit, double frac) {
 159                 if (read.isLowQuality() || frac < EPSILON) return;
 160
 161                 RefSeq& ref = refs->getRef(hit.getSid());
 162                 int dir = hit.getDir();
 163                 int pos = hit.getPos();
 164
 165                 if (estRSPD) {
 166                         int fullLen = ref.getFullLen();
 167
 168                         // Only use one strand to estimate RSPD
 169                         if (ori->getProb(0) >= ORIVALVE && dir == 0) {
 170                                 rspd->update(pos, fullLen, frac);
 171                         }
 172
 173                         if (ori->getProb(0) < ORIVALVE && dir == 1) {
 174                                 int totLen = ref.getTotLen();
 175                                 int readLen = read.getReadLength();
 176
 177                                 int pfpos, effL;
 178
 179                                 if (mld != NULL) {
 180                                         int minL = std::max(readLen, gld->getMinL());
 181                                         int maxL = std::min(totLen - pos, gld->getMaxL());
 182                                         double sum = 0.0;
 183                                         assert(maxL >= minL);
 184                                         std::vector<double> frag_vec(maxL - minL + 1, 0.0);
 185
 186                                         for (int fragLen = minL; fragLen <= maxL; fragLen++) {
 187                                                 pfpos = totLen - pos - fragLen;
 188                                                 effL = std::min(fullLen, totLen - fragLen + 1);
 189                                                 frag_vec[fragLen - minL] = gld->getAdjustedProb(fragLen, totLen) * rspd->getAdjustedProb(pfpos, effL, fullLen) * mld->getAdjustedProb(readLen, fragLen);
 190                                                 sum += frag_vec[fragLen - minL];
 191                                         }
 192                                         assert(sum >= EPSILON);
 193                                         for (int fragLen = minL; fragLen <= maxL; fragLen++) {
 194                                                 pfpos = totLen - pos - fragLen;
 195                                                 rspd->update(pfpos, fullLen, frac * (frag_vec[fragLen - minL] / sum));
 196                                         }
 197                                 }
 198                                 else {
 199                                         rspd->update(totLen - pos - readLen, fullLen, frac);
 200                                 }
 201                         }
 202                 }
 203                 qpro->update(read.getReadSeq(), read.getQScore(), ref, pos, dir, frac);
 204         }
 205
 206         void updateNoise(const SingleReadQ& read, double frac) {
 207                 if (read.isLowQuality() || frac < EPSILON) return;
 208
 209                 nqpro->update(read.getReadSeq(), read.getQScore(), frac);
 210         }
 211
 212         void finish();
 213
 214         void collect(const SingleQModel&);
 215
 216         //void copy(const SingleQModel&);
 217
 218         bool getNeedCalcConPrb() { return needCalcConPrb; }
 219         void setNeedCalcConPrb(bool value) { needCalcConPrb = value; }
 220
 221         //void calcP1();
 222         //void calcP2();
 223         //double* getP1() { return p1; }
 224         //double* getP2() { return p2; }
 225
 226         void read(const char*);
 227         void write(const char*);
 228
 229         const LenDist& getGLD() { return *gld; }
 230
 231         void startSimulation(simul*, double*);
 232         bool simulate(int, SingleReadQ&, int&);
 233         void finishSimulation();
 234
 235         //Use it after function 'read' or 'estimateFromReads'
 236         double* getMW() {
 237           assert(mw != NULL);
 238           return mw;
 239         }
 240
 241         int getModelType() const { return model_type; }
 242
 243 private:
 244         static const int model_type = 1;
 245         static const int read_type = 1;
 246
 247         int M;
 248         int N[3];
 249         Refs *refs;
 250         double mean, sd;
 251         int seedLen;
 252         //double *p1, *p2; P_i' & P_i'';
 253
 254         bool estRSPD; // true if estimate RSPD
 255         bool needCalcConPrb; //true need, false does not need
 256
 257         Orientation *ori;
 258         LenDist *gld, *mld;
 259         RSPD *rspd;
 260         QualDist *qd;
 261         QProfile *qpro;
 262         NoiseQProfile *nqpro;
 263
 264         simul *sampler; // for simulation
 265         double *theta_cdf; // for simulation
 266
 267         double *mw; // for masking
 268
 269         void calcMW();
 270 };
 271
 272 void SingleQModel::estimateFromReads(const char* readFN) {
 273         int s;
 274         char readFs[2][STRLEN];
 275         SingleReadQ read;
 276
 277         mld != NULL ? mld->init() : gld->init();
 278         for (int i = 0; i < 3; i++)
 279                 if (N[i] > 0) {
 280                         genReadFileNames(readFN, i, read_type, s, readFs);
 281                         ReadReader<SingleReadQ> reader(s, readFs);
 282
 283                         int cnt = 0;
 284                         while (reader.next(read)) {
 285                                 mld != NULL ? mld->update(read.getReadLength(), 1.0) : gld->update(read.getReadLength(), 1.0);
 286                                 qd->update(read.getQScore());
 287                                 if (i == 0) { nqpro->updateC(read.getReadSeq(), read.getQScore()); }
 288
 289                                 ++cnt;
 290                                 if (verbose && cnt % 1000000 == 0) { printf("%d READS PROCESSED\n", cnt); }
 291                         }
 292
 293                         if (verbose) { printf("estimateFromReads, N%d finished.\n", i); }
 294                 }
 295
 296         mld != NULL ? mld->finish() : gld->finish();
 297
 298         //mean should be > 0
 299         if (mean >= EPSILON) {
 300           assert(mld->getMaxL() <= gld->getMaxL());
 301           gld->setAsNormal(mean, sd, std::max(mld->getMinL(), gld->getMinL()), gld->getMaxL());
 302         }
 303         qd->finish();
 304         nqpro->calcInitParams();
 305
 306         mw = new double[M + 1];
 307         calcMW();
 308 }
 309
 310 void SingleQModel::init() {
 311         if (estRSPD) rspd->init();
 312         qpro->init();
 313         nqpro->init();
 314 }
 315
 316 void SingleQModel::finish() {
 317         if (estRSPD) rspd->finish();
 318         qpro->finish();
 319         nqpro->finish();
 320         needCalcConPrb = true;
 321         if (estRSPD) calcMW();
 322 }
 323
 324 void SingleQModel::collect(const SingleQModel& o) {
 325         if (estRSPD) rspd->collect(*(o.rspd));
 326         qpro->collect(*(o.qpro));
 327         nqpro->collect(*(o.nqpro));
 328 }
 329
 330 //Only master node can call
 331 void SingleQModel::read(const char* inpF) {
 332         int val;
 333         FILE *fi = fopen(inpF, "r");
 334         if (fi == NULL) { fprintf(stderr, "Cannot open %s! It may not exist.\n", inpF); exit(-1); }
 335
 336         fscanf(fi, "%d", &val);
 337         assert(val == model_type);
 338
 339         ori->read(fi);
 340         gld->read(fi);
 341         fscanf(fi, "%d", &val);
 342         if (val > 0) {
 343                 if (mld == NULL) mld = new LenDist();
 344                 mld->read(fi);
 345         }
 346         rspd->read(fi);
 347         qd->read(fi);
 348         qpro->read(fi);
 349         nqpro->read(fi);
 350
 351         if (fscanf(fi, "%d", &M) == 1) {
 352           mw = new double[M + 1];
 353           for (int i = 0; i <= M; i++) fscanf(fi, "%lf", &mw[i]);
 354         }
 355
 356         fclose(fi);
 357 }
 358
 359 //Only master node can call
 360 void SingleQModel::write(const char* outF) {
 361         FILE *fo = fopen(outF, "w");
 362
 363         fprintf(fo, "%d\n", model_type);
 364         fprintf(fo, "\n");
 365
 366         ori->write(fo);  fprintf(fo, "\n");
 367         gld->write(fo);  fprintf(fo, "\n");
 368         if (mld != NULL) {
 369                 fprintf(fo, "1\n");
 370                 mld->write(fo);
 371         }
 372         else { fprintf(fo, "0\n"); }
 373         fprintf(fo, "\n");
 374         rspd->write(fo); fprintf(fo, "\n");
 375         qd->write(fo);   fprintf(fo, "\n");
 376         qpro->write(fo); fprintf(fo, "\n");
 377         nqpro->write(fo);
 378
 379         if (mw != NULL) {
 380           fprintf(fo, "\n%d\n", M);
 381           for (int i = 0; i < M; i++) {
 382             fprintf(fo, "%.15g ", mw[i]);
 383           }
 384           fprintf(fo, "%.15g\n", mw[M]);
 385         }
 386
 387         fclose(fo);
 388 }
 389
 390 void SingleQModel::startSimulation(simul* sampler, double* theta) {
 391         this->sampler = sampler;
 392
 393         theta_cdf = new double[M + 1];
 394         for (int i = 0; i <= M; i++) {
 395                 theta_cdf[i] = theta[i];
 396                 if (i > 0) theta_cdf[i] += theta_cdf[i - 1];
 397         }
 398
 399         rspd->startSimulation(M, refs);
 400         qd->startSimulation();
 401         qpro->startSimulation();
 402         nqpro->startSimulation();
 403 }
 404
 405 bool SingleQModel::simulate(int rid, SingleReadQ& read, int& sid) {
 406         int dir, pos, readLen, fragLen;
 407         std::string name;
 408         std::string qual, readseq;
 409         std::ostringstream strout;
 410
 411         sid = sampler->sample(theta_cdf, M + 1);
 412
 413         if (sid == 0) {
 414                 dir = pos = 0;
 415                 readLen = (mld != NULL ? mld->simulate(sampler, -1) : gld->simulate(sampler, -1));
 416                 qual = qd->simulate(sampler, readLen);
 417                 readseq = nqpro->simulate(sampler, readLen, qual);
 418         }
 419         else {
 420                 RefSeq &ref = refs->getRef(sid);
 421                 dir = ori->simulate(sampler);
 422                 fragLen = gld->simulate(sampler, ref.getTotLen());
 423                 if (fragLen < 0) return false;
 424
 425                 int effL = std::min(ref.getFullLen(), ref.getTotLen() - fragLen + 1);
 426                 pos = rspd->simulate(sampler, sid, effL);
 427                 if (pos < 0) return false;
 428                 if (dir > 0) pos = ref.getTotLen() - pos - fragLen;
 429
 430                 if (mld != NULL) {
 431                         readLen = mld->simulate(sampler, fragLen);
 432                         if (readLen < 0) return false;
 433                         qual = qd->simulate(sampler, readLen);
 434                         readseq = qpro->simulate(sampler, readLen, pos, dir, qual, ref);
 435                 }
 436                 else {
 437                         qual = qd->simulate(sampler, fragLen);
 438                         readseq = qpro->simulate(sampler, fragLen, pos, dir, qual, ref);
 439                 }
 440         }
 441
 442         std::ostringstream stdout;
 443         stdout<<rid<<"_"<<dir<<"_"<<sid<<"_"<<pos;
 444         name = stdout.str();
 445
 446         read = SingleReadQ(name, readseq, qual);
 447
 448         return true;
 449 }
 450
 451 void SingleQModel::finishSimulation() {
 452         delete[] theta_cdf;
 453
 454         rspd->finishSimulation();
 455         qd->finishSimulation();
 456         qpro->finishSimulation();
 457         nqpro->finishSimulation();
 458 }
 459
 460 void SingleQModel::calcMW() {
 461   double probF, probR;
 462
 463   assert(seedLen >= OLEN && (mld == NULL ? gld->getMinL() : mld->getMinL()) >= seedLen);
 464
 465   memset(mw, 0, sizeof(double) * (M + 1));
 466   mw[0] = 1.0;
 467
 468
 469   probF = ori->getProb(0);
 470   probR = ori->getProb(1);
 471
 472   for (int i = 1; i <= M; i++) {
 473     RefSeq& ref = refs->getRef(i);
 474     int totLen = ref.getTotLen();
 475     int fullLen = ref.getFullLen();
 476     double value = 0.0;
 477     int minL, maxL;
 478     int effL, pfpos;
 479     int end = std::min(fullLen, totLen - seedLen + 1);
 480     double factor;
 481
 482     for (int seedPos = 0; seedPos < end; seedPos++)
 483       if (ref.getMask(seedPos)) {
 484         //forward
 485         minL = gld->getMinL();
 486         maxL = std::min(gld->getMaxL(), totLen - seedPos);
 487         pfpos = seedPos;
 488         for (int fragLen = minL; fragLen <= maxL; fragLen++) {
 489           effL = std::min(fullLen, totLen - fragLen + 1);
 490           factor = (mld == NULL ? 1.0 : mld->getAdjustedCumulativeProb(std::min(mld->getMaxL(), fragLen), fragLen));
 491           value += probF * gld->getAdjustedProb(fragLen, totLen) * rspd->getAdjustedProb(pfpos, effL, fullLen) * factor;
 492         }
 493         //reverse
 494         minL = gld->getMinL();
 495         maxL = std::min(gld->getMaxL(), seedPos + seedLen);
 496         for (int fragLen = minL; fragLen <= maxL; fragLen++) {
 497           pfpos = seedPos - (fragLen - seedLen);
 498           effL = std::min(fullLen, totLen - fragLen + 1);
 499           factor = (mld == NULL ? 1.0 : mld->getAdjustedCumulativeProb(std::min(mld->getMaxL(), fragLen), fragLen));
 500           value += probR * gld->getAdjustedProb(fragLen, totLen) * rspd->getAdjustedProb(pfpos, effL, fullLen) * factor;
 501         }
 502       }
 503
 504     //for reverse strand masking
 505     for (int seedPos = end; seedPos <= totLen - seedLen; seedPos++) {
 506       minL = std::max(gld->getMinL(), seedPos + seedLen - fullLen + 1);
 507       maxL = std::min(gld->getMaxL(), seedPos + seedLen);
 508       for (int fragLen = minL; fragLen <= maxL; fragLen++) {
 509         pfpos = seedPos - (fragLen - seedLen);
 510         effL = std::min(fullLen, totLen - fragLen + 1);
 511         factor = (mld == NULL ? 1.0 : mld->getAdjustedCumulativeProb(std::min(mld->getMaxL(), fragLen), fragLen));
 512         value += probR * gld->getAdjustedProb(fragLen, totLen) * rspd->getAdjustedProb(pfpos, effL, fullLen) * factor;
 513       }
 514     }
 515
 516     mw[i] = 1.0 - value;
 517
 518     if (mw[i] < 1e-8) {
 519       //      fprintf(stderr, "Warning: %dth reference sequence is masked for almost all positions!\n", i);
 520       mw[i] = 0.0;
 521     }
 522   }
 523 }
 524
 525 #endif /* SINGLEQMODEL_H_ */