SingleModel.h

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