1 #ifndef PAIREDENDQMODEL_H_
2 #define PAIREDENDQMODEL_H_
13 #include "Orientation.h"
18 #include "NoiseQProfile.h"
20 #include "ModelParams.h"
23 #include "SingleReadQ.h"
24 #include "PairedEndReadQ.h"
25 #include "PairedEndHit.h"
26 #include "ReadReader.h"
30 class PairedEndQModel {
32 PairedEndQModel(Refs* refs = NULL) {
34 M = (refs != NULL ? refs->getM() : 0);
35 memset(N, 0, sizeof(N));
37 needCalcConPrb = true;
39 ori = new Orientation();
41 rspd = new RSPD(estRSPD);
43 qpro = new QProfile();
44 nqpro = new NoiseQProfile();
51 //If it is not a master node, only init & update can be used!
52 PairedEndQModel(ModelParams& params, bool isMaster = true) {
54 memcpy(N, params.N, sizeof(params.N));
56 estRSPD = params.estRSPD;
57 seedLen = params.seedLen;
58 needCalcConPrb = true;
60 ori = NULL; gld = NULL; rspd = NULL; qd = NULL; qpro = NULL; nqpro = NULL; mld = NULL;
64 if (!estRSPD) rspd = new RSPD(estRSPD);
66 mld = new LenDist(params.mate_minL, params.mate_maxL);
69 ori = new Orientation(params.probF);
70 gld = new LenDist(params.minL, params.maxL);
71 if (estRSPD) rspd = new RSPD(estRSPD, params.B);
72 qpro = new QProfile();
73 nqpro = new NoiseQProfile();
78 if (ori != NULL) delete ori;
79 if (gld != NULL) delete gld;
80 if (rspd != NULL) delete rspd;
81 if (qd != NULL) delete qd;
82 if (qpro != NULL) delete qpro;
83 if (nqpro != NULL) delete nqpro;
84 if (mld != NULL) delete mld;
85 if (mw != NULL) delete mw;
88 void estimateFromReads(const char*);
90 //if prob is too small, just make it 0
91 double getConPrb(const PairedEndReadQ& read, const PairedEndHit& hit) {
92 if (read.isLowQuality()) return 0.0;
95 int sid = hit.getSid();
96 RefSeq &ref = refs->getRef(sid);
97 int dir = hit.getDir();
98 int pos = hit.getPos();
99 int fullLen = ref.getFullLen();
100 int totLen = ref.getTotLen();
101 int insertLen = hit.getInsertL();
103 int fpos = (dir == 0 ? pos : totLen - pos - insertLen); // the aligned position reported in SAM file, should be a coordinate in forward strand
104 int effL = std::min(fullLen, totLen - insertLen + 1);
106 assert(fpos >= 0 && fpos + insertLen <= totLen && insertLen <= totLen);
107 if (fpos >= fullLen || ref.getMask(fpos)) return 0.0; // For paired-end model, fpos is the seedPos
109 prob = ori->getProb(dir) * gld->getAdjustedProb(insertLen, totLen) *
110 rspd->getAdjustedProb(fpos, effL, fullLen);
112 const SingleReadQ& mate1 = read.getMate1();
113 prob *= mld->getAdjustedProb(mate1.getReadLength(), insertLen) *
114 qpro->getProb(mate1.getReadSeq(), mate1.getQScore(), ref, pos, dir);
116 const SingleReadQ& mate2 = read.getMate2();
117 int m2pos = totLen - pos - insertLen;
119 prob *= mld->getAdjustedProb(mate2.getReadLength(), hit.getInsertL()) *
120 qpro->getProb(mate2.getReadSeq(), mate2.getQScore(), ref, m2pos, m2dir);
122 if (prob < EPSILON) { prob = 0.0; }
124 prob = (mw[sid] < EPSILON ? 0.0 : prob / mw[sid]);
129 double getNoiseConPrb(const PairedEndReadQ& read) {
130 if (read.isLowQuality()) return 0.0;
133 const SingleReadQ& mate1 = read.getMate1();
134 const SingleReadQ& mate2 = read.getMate2();
136 prob = mld->getProb(mate1.getReadLength()) * nqpro->getProb(mate1.getReadSeq(), mate1.getQScore());
137 prob *= mld->getProb(mate2.getReadLength()) * nqpro->getProb(mate2.getReadSeq(), mate2.getQScore());
139 if (prob < EPSILON) { prob = 0.0; }
141 prob = (mw[0] < EPSILON ? 0.0: prob / mw[0]);
146 double getLogP() { return nqpro->getLogP(); }
150 void update(const PairedEndReadQ& read, const PairedEndHit& hit, double frac) {
151 if (read.isLowQuality() || frac < EPSILON) return;
153 RefSeq& ref = refs->getRef(hit.getSid());
154 const SingleReadQ& mate1 = read.getMate1();
155 const SingleReadQ& mate2 = read.getMate2();
157 gld->update(hit.getInsertL(), frac);
159 int fpos = (hit.getDir() == 0 ? hit.getPos() : ref.getTotLen() - hit.getPos() - hit.getInsertL());
160 rspd->update(fpos, ref.getFullLen(), frac);
162 qpro->update(mate1.getReadSeq(), mate1.getQScore(), ref, hit.getPos(), hit.getDir(), frac);
164 int m2pos = ref.getTotLen() - hit.getPos() - hit.getInsertL();
165 int m2dir = !hit.getDir();
166 qpro->update(mate2.getReadSeq(), mate2.getQScore(), ref, m2pos, m2dir, frac);
169 void updateNoise(const PairedEndReadQ& read, double frac) {
170 if (read.isLowQuality() || frac < EPSILON) return;
172 const SingleReadQ& mate1 = read.getMate1();
173 const SingleReadQ& mate2 = read.getMate2();
175 nqpro->update(mate1.getReadSeq(), mate1.getQScore(), frac);
176 nqpro->update(mate2.getReadSeq(), mate2.getQScore(), frac);
181 void collect(const PairedEndQModel&);
183 bool getNeedCalcConPrb() { return needCalcConPrb; }
184 void setNeedCalcConPrb(bool value) { needCalcConPrb = value; }
186 void read(const char*);
187 void write(const char*);
189 const LenDist& getGLD() { return *gld; }
191 void startSimulation(simul*, double*);
192 bool simulate(int, PairedEndReadQ&, int&);
193 void finishSimulation();
195 //Use it after function 'read' or 'estimateFromReads'
201 int getModelType() const { return model_type; }
204 static const int model_type = 3;
205 static const int read_type = 3;
213 bool needCalcConPrb; //true need, false does not need
216 LenDist *gld, *mld; //mld1 mate_length_dist
220 NoiseQProfile *nqpro;
222 simul *sampler; // for simulation
223 double *theta_cdf; // for simulation
225 double *mw; // for masking
230 void PairedEndQModel::estimateFromReads(const char* readFN) {
232 char readFs[2][STRLEN];
236 for (int i = 0; i < 3; i++)
238 genReadFileNames(readFN, i, read_type, s, readFs);
239 ReadReader<PairedEndReadQ> reader(s, readFs, refs->hasPolyA(), seedLen); // allow calculation of calc_lq() function
242 while (reader.next(read)) {
243 SingleReadQ mate1 = read.getMate1();
244 SingleReadQ mate2 = read.getMate2();
246 if (!read.isLowQuality()) {
247 mld->update(mate1.getReadLength(), 1.0);
248 mld->update(mate2.getReadLength(), 1.0);
250 qd->update(mate1.getQScore());
251 qd->update(mate2.getQScore());
254 nqpro->updateC(mate1.getReadSeq(), mate1.getQScore());
255 nqpro->updateC(mate2.getReadSeq(), mate2.getQScore());
258 else if (verbose && (mate1.getReadLength() < seedLen || mate2.getReadLength() < seedLen)) {
259 printf("Warning: Read %s is ignored due to at least one of the mates' length < seed length %d!\n", read.getName().c_str(), seedLen);
263 if (verbose && cnt % 1000000 == 0) { printf("%d READS PROCESSED\n", cnt); }
266 if (verbose) { printf("estimateFromReads, N%d finished.\n", i); }
271 nqpro->calcInitParams();
273 mw = new double[M + 1];
277 void PairedEndQModel::init() {
279 if (estRSPD) rspd->init();
284 void PairedEndQModel::finish() {
286 if (estRSPD) rspd->finish();
289 needCalcConPrb = true;
293 void PairedEndQModel::collect(const PairedEndQModel& o) {
294 gld->collect(*(o.gld));
295 if (estRSPD) rspd->collect(*(o.rspd));
296 qpro->collect(*(o.qpro));
297 nqpro->collect(*(o.nqpro));
300 //Only master node can call
301 void PairedEndQModel::read(const char* inpF) {
303 FILE *fi = fopen(inpF, "r");
304 if (fi == NULL) { fprintf(stderr, "Cannot open %s! It may not exist.\n", inpF); exit(-1); }
306 assert(fscanf(fi, "%d", &val) == 1);
307 assert(val == model_type);
317 if (fscanf(fi, "%d", &val) == 1) {
320 mw = new double[M + 1];
321 for (int i = 0; i <= M; i++) assert(fscanf(fi, "%lf", &mw[i]) == 1);
329 //Only master node can call. Only be called at EM.cpp
330 void PairedEndQModel::write(const char* outF) {
331 FILE *fo = fopen(outF, "w");
333 fprintf(fo, "%d\n", model_type);
336 ori->write(fo); fprintf(fo, "\n");
337 gld->write(fo); fprintf(fo, "\n");
338 mld->write(fo); fprintf(fo, "\n");
339 rspd->write(fo); fprintf(fo, "\n");
340 qd->write(fo); fprintf(fo, "\n");
341 qpro->write(fo); fprintf(fo, "\n");
345 fprintf(fo, "\n%d\n", M);
346 for (int i = 0; i < M; i++) {
347 fprintf(fo, "%.15g ", mw[i]);
349 fprintf(fo, "%.15g\n", mw[M]);
355 void PairedEndQModel::startSimulation(simul* sampler, double* theta) {
356 this->sampler = sampler;
358 theta_cdf = new double[M + 1];
359 for (int i = 0; i <= M; i++) {
360 theta_cdf[i] = theta[i];
361 if (i > 0) theta_cdf[i] += theta_cdf[i - 1];
364 rspd->startSimulation(M, refs);
365 qd->startSimulation();
366 qpro->startSimulation();
367 nqpro->startSimulation();
370 bool PairedEndQModel::simulate(int rid, PairedEndReadQ& read, int& sid) {
372 int insertL, mateL1, mateL2;
374 std::string qual1, qual2, readseq1, readseq2;
375 std::ostringstream strout;
377 sid = sampler->sample(theta_cdf, M + 1);
380 dir = pos = insertL = 0;
381 mateL1 = mld->simulate(sampler, -1);
382 qual1 = qd->simulate(sampler, mateL1);
383 readseq1 = nqpro->simulate(sampler, mateL1, qual1);
385 mateL2 = mld->simulate(sampler, -1);
386 qual2 = qd->simulate(sampler, mateL2);
387 readseq2 = nqpro->simulate(sampler, mateL2, qual2);
390 RefSeq &ref = refs->getRef(sid);
391 dir = ori->simulate(sampler);
392 insertL = gld->simulate(sampler, ref.getTotLen());
393 if (insertL < 0) return false;
394 int effL = std::min(ref.getFullLen(), ref.getTotLen() - insertL + 1);
395 pos = rspd->simulate(sampler, sid, effL);
396 if (pos < 0) return false;
397 if (dir > 0) pos = ref.getTotLen() - pos - insertL;
399 mateL1 = mld->simulate(sampler, insertL);
400 qual1 = qd->simulate(sampler, mateL1);
401 readseq1 = qpro->simulate(sampler, mateL1, pos, dir, qual1, ref);
403 int m2pos = ref.getTotLen() - pos - insertL;
406 mateL2 = mld->simulate(sampler, insertL);
407 qual2 = qd->simulate(sampler, mateL2);
408 readseq2 = qpro->simulate(sampler, mateL2, m2pos, m2dir, qual2, ref);
411 strout<<rid<<"_"<<dir<<"_"<<sid<<"_"<<pos<<"_"<<insertL;
414 read = PairedEndReadQ(SingleReadQ(name + "/1", readseq1, qual1), SingleReadQ(name + "/2", readseq2, qual2));
419 void PairedEndQModel::finishSimulation() {
422 rspd->finishSimulation();
423 qd->finishSimulation();
424 qpro->finishSimulation();
425 nqpro->finishSimulation();
429 void PairedEndQModel::calcMW() {
430 assert(mld->getMinL() >= seedLen);
432 memset(mw, 0, sizeof(double) * (M + 1));
435 for (int i = 1; i <= M; i++) {
436 RefSeq& ref = refs->getRef(i);
437 int totLen = ref.getTotLen();
438 int fullLen = ref.getFullLen();
439 int end = std::min(fullLen, totLen - gld->getMinL() + 1);
444 //seedPos is fpos here
445 for (int seedPos = 0; seedPos < end; seedPos++)
446 if (ref.getMask(seedPos)) {
447 minL = gld->getMinL();
448 maxL = std::min(gld->getMaxL(), totLen - seedPos);
450 for (int fragLen = minL; fragLen <= maxL; fragLen++) {
451 effL = std::min(fullLen, totLen - fragLen + 1);
452 value += gld->getAdjustedProb(fragLen, totLen) * rspd->getAdjustedProb(pfpos, effL, fullLen);
459 // fprintf(stderr, "Warning: %dth reference sequence is masked for almost all positions!\n", i);
465 #endif /* PAIREDENDQMODEL_H_ */