PairedEndQModel.h

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