PairedEndModel.h

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