1 #ifndef PAIREDENDQMODEL_H_
2 #define PAIREDENDQMODEL_H_
14 #include "my_assert.h"
15 #include "Orientation.h"
20 #include "NoiseQProfile.h"
22 #include "ModelParams.h"
25 #include "SingleReadQ.h"
26 #include "PairedEndReadQ.h"
27 #include "PairedEndHit.h"
28 #include "ReadReader.h"
32 class PairedEndQModel {
34 PairedEndQModel(Refs* refs = NULL) {
36 M = (refs != NULL ? refs->getM() : 0);
37 memset(N, 0, sizeof(N));
39 needCalcConPrb = true;
41 ori = new Orientation();
43 rspd = new RSPD(estRSPD);
45 qpro = new QProfile();
46 nqpro = new NoiseQProfile();
53 //If it is not a master node, only init & update can be used!
54 PairedEndQModel(ModelParams& params, bool isMaster = true) {
56 memcpy(N, params.N, sizeof(params.N));
58 estRSPD = params.estRSPD;
59 seedLen = params.seedLen;
60 needCalcConPrb = true;
62 ori = NULL; gld = NULL; rspd = NULL; qd = NULL; qpro = NULL; nqpro = NULL; mld = NULL;
66 if (!estRSPD) rspd = new RSPD(estRSPD);
68 mld = new LenDist(params.mate_minL, params.mate_maxL);
71 ori = new Orientation(params.probF);
72 gld = new LenDist(params.minL, params.maxL);
73 if (estRSPD) rspd = new RSPD(estRSPD, params.B);
74 qpro = new QProfile();
75 nqpro = new NoiseQProfile();
80 if (ori != NULL) delete ori;
81 if (gld != NULL) delete gld;
82 if (rspd != NULL) delete rspd;
83 if (qd != NULL) delete qd;
84 if (qpro != NULL) delete qpro;
85 if (nqpro != NULL) delete nqpro;
86 if (mld != NULL) delete mld;
87 if (mw != NULL) delete mw;
90 void estimateFromReads(const char*);
92 //if prob is too small, just make it 0
93 double getConPrb(const PairedEndReadQ& read, const PairedEndHit& hit) {
94 if (read.isLowQuality()) return 0.0;
97 int sid = hit.getSid();
98 RefSeq &ref = refs->getRef(sid);
99 int dir = hit.getDir();
100 int pos = hit.getPos();
101 int fullLen = ref.getFullLen();
102 int totLen = ref.getTotLen();
103 int insertLen = hit.getInsertL();
105 int fpos = (dir == 0 ? pos : totLen - pos - insertLen); // the aligned position reported in SAM file, should be a coordinate in forward strand
106 int effL = std::min(fullLen, totLen - insertLen + 1);
108 general_assert(fpos >= 0, "The alignment of fragment " + 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 + insertLen <= totLen,"Fragment " + 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(insertLen <= totLen, "Fragment " + read.getName() + " has length " + itos(insertLen) + ", but it is aligned to transcript " \
114 + itos(sid) + ", whose length (" + itos(totLen) + ") is shorter than the fragment's length!");
116 if (fpos >= fullLen || ref.getMask(fpos)) return 0.0; // For paired-end model, fpos is the seedPos
118 prob = ori->getProb(dir) * gld->getAdjustedProb(insertLen, totLen) *
119 rspd->getAdjustedProb(fpos, effL, fullLen);
121 const SingleReadQ& mate1 = read.getMate1();
122 prob *= mld->getAdjustedProb(mate1.getReadLength(), insertLen) *
123 qpro->getProb(mate1.getReadSeq(), mate1.getQScore(), ref, pos, dir);
125 const SingleReadQ& mate2 = read.getMate2();
126 int m2pos = totLen - pos - insertLen;
128 prob *= mld->getAdjustedProb(mate2.getReadLength(), hit.getInsertL()) *
129 qpro->getProb(mate2.getReadSeq(), mate2.getQScore(), ref, m2pos, m2dir);
131 if (prob < EPSILON) { prob = 0.0; }
133 prob = (mw[sid] < EPSILON ? 0.0 : prob / mw[sid]);
138 double getNoiseConPrb(const PairedEndReadQ& read) {
139 if (read.isLowQuality()) return 0.0;
142 const SingleReadQ& mate1 = read.getMate1();
143 const SingleReadQ& mate2 = read.getMate2();
145 prob = mld->getProb(mate1.getReadLength()) * nqpro->getProb(mate1.getReadSeq(), mate1.getQScore());
146 prob *= mld->getProb(mate2.getReadLength()) * nqpro->getProb(mate2.getReadSeq(), mate2.getQScore());
148 if (prob < EPSILON) { prob = 0.0; }
150 prob = (mw[0] < EPSILON ? 0.0: prob / mw[0]);
155 double getLogP() { return nqpro->getLogP(); }
159 void update(const PairedEndReadQ& read, const PairedEndHit& hit, double frac) {
160 if (read.isLowQuality() || frac < EPSILON) return;
162 RefSeq& ref = refs->getRef(hit.getSid());
163 const SingleReadQ& mate1 = read.getMate1();
164 const SingleReadQ& mate2 = read.getMate2();
166 gld->update(hit.getInsertL(), frac);
168 int fpos = (hit.getDir() == 0 ? hit.getPos() : ref.getTotLen() - hit.getPos() - hit.getInsertL());
169 rspd->update(fpos, ref.getFullLen(), frac);
171 qpro->update(mate1.getReadSeq(), mate1.getQScore(), ref, hit.getPos(), hit.getDir(), frac);
173 int m2pos = ref.getTotLen() - hit.getPos() - hit.getInsertL();
174 int m2dir = !hit.getDir();
175 qpro->update(mate2.getReadSeq(), mate2.getQScore(), ref, m2pos, m2dir, frac);
178 void updateNoise(const PairedEndReadQ& read, double frac) {
179 if (read.isLowQuality() || frac < EPSILON) return;
181 const SingleReadQ& mate1 = read.getMate1();
182 const SingleReadQ& mate2 = read.getMate2();
184 nqpro->update(mate1.getReadSeq(), mate1.getQScore(), frac);
185 nqpro->update(mate2.getReadSeq(), mate2.getQScore(), frac);
190 void collect(const PairedEndQModel&);
192 bool getNeedCalcConPrb() { return needCalcConPrb; }
193 void setNeedCalcConPrb(bool value) { needCalcConPrb = value; }
195 void read(const char*);
196 void write(const char*);
198 const LenDist& getGLD() { return *gld; }
200 void startSimulation(simul*, double*);
201 bool simulate(READ_INT_TYPE, PairedEndReadQ&, int&);
202 void finishSimulation();
204 //Use it after function 'read' or 'estimateFromReads'
210 int getModelType() const { return model_type; }
213 static const int model_type = 3;
214 static const int read_type = 3;
222 bool needCalcConPrb; //true need, false does not need
225 LenDist *gld, *mld; //mld1 mate_length_dist
229 NoiseQProfile *nqpro;
231 simul *sampler; // for simulation
232 double *theta_cdf; // for simulation
234 double *mw; // for masking
239 void PairedEndQModel::estimateFromReads(const char* readFN) {
241 char readFs[2][STRLEN];
245 for (int i = 0; i < 3; i++)
247 genReadFileNames(readFN, i, read_type, s, readFs);
248 ReadReader<PairedEndReadQ> reader(s, readFs, refs->hasPolyA(), seedLen); // allow calculation of calc_lq() function
250 READ_INT_TYPE cnt = 0;
251 while (reader.next(read)) {
252 SingleReadQ mate1 = read.getMate1();
253 SingleReadQ mate2 = read.getMate2();
255 if (!read.isLowQuality()) {
256 mld->update(mate1.getReadLength(), 1.0);
257 mld->update(mate2.getReadLength(), 1.0);
259 qd->update(mate1.getQScore());
260 qd->update(mate2.getQScore());
263 nqpro->updateC(mate1.getReadSeq(), mate1.getQScore());
264 nqpro->updateC(mate2.getReadSeq(), mate2.getQScore());
267 else if (verbose && (mate1.getReadLength() < seedLen || mate2.getReadLength() < seedLen)) {
268 std::cout<< "Warning: Read "<< read.getName()<< " is ignored due to at least one of the mates' length < seed length "<< seedLen<< "!"<< std::endl;
272 if (verbose && cnt % 1000000 == 0) { std::cout<< cnt<< " READS PROCESSED"<< std::endl; }
275 if (verbose) { std::cout<<"estimateFromReads, N"<< i<<" finished."<< std::endl; }
280 nqpro->calcInitParams();
282 mw = new double[M + 1];
286 void PairedEndQModel::init() {
288 if (estRSPD) rspd->init();
293 void PairedEndQModel::finish() {
295 if (estRSPD) rspd->finish();
298 needCalcConPrb = true;
302 void PairedEndQModel::collect(const PairedEndQModel& o) {
303 gld->collect(*(o.gld));
304 if (estRSPD) rspd->collect(*(o.rspd));
305 qpro->collect(*(o.qpro));
306 nqpro->collect(*(o.nqpro));
309 //Only master node can call
310 void PairedEndQModel::read(const char* inpF) {
312 FILE *fi = fopen(inpF, "r");
313 if (fi == NULL) { fprintf(stderr, "Cannot open %s! It may not exist.\n", inpF); exit(-1); }
315 assert(fscanf(fi, "%d", &val) == 1);
316 assert(val == model_type);
326 if (fscanf(fi, "%d", &val) == 1) {
329 mw = new double[M + 1];
330 for (int i = 0; i <= M; i++) assert(fscanf(fi, "%lf", &mw[i]) == 1);
338 //Only master node can call. Only be called at EM.cpp
339 void PairedEndQModel::write(const char* outF) {
340 FILE *fo = fopen(outF, "w");
342 fprintf(fo, "%d\n", model_type);
345 ori->write(fo); fprintf(fo, "\n");
346 gld->write(fo); fprintf(fo, "\n");
347 mld->write(fo); fprintf(fo, "\n");
348 rspd->write(fo); fprintf(fo, "\n");
349 qd->write(fo); fprintf(fo, "\n");
350 qpro->write(fo); fprintf(fo, "\n");
354 fprintf(fo, "\n%d\n", M);
355 for (int i = 0; i < M; i++) {
356 fprintf(fo, "%.15g ", mw[i]);
358 fprintf(fo, "%.15g\n", mw[M]);
364 void PairedEndQModel::startSimulation(simul* sampler, double* theta) {
365 this->sampler = sampler;
367 theta_cdf = new double[M + 1];
368 for (int i = 0; i <= M; i++) {
369 theta_cdf[i] = theta[i];
370 if (i > 0) theta_cdf[i] += theta_cdf[i - 1];
373 rspd->startSimulation(M, refs);
374 qd->startSimulation();
375 qpro->startSimulation();
376 nqpro->startSimulation();
379 bool PairedEndQModel::simulate(READ_INT_TYPE rid, PairedEndReadQ& read, int& sid) {
381 int insertL, mateL1, mateL2;
383 std::string qual1, qual2, readseq1, readseq2;
384 std::ostringstream strout;
386 sid = sampler->sample(theta_cdf, M + 1);
389 dir = pos = insertL = 0;
390 mateL1 = mld->simulate(sampler, -1);
391 qual1 = qd->simulate(sampler, mateL1);
392 readseq1 = nqpro->simulate(sampler, mateL1, qual1);
394 mateL2 = mld->simulate(sampler, -1);
395 qual2 = qd->simulate(sampler, mateL2);
396 readseq2 = nqpro->simulate(sampler, mateL2, qual2);
399 RefSeq &ref = refs->getRef(sid);
400 dir = ori->simulate(sampler);
401 insertL = gld->simulate(sampler, ref.getTotLen());
402 if (insertL < 0) return false;
403 int effL = std::min(ref.getFullLen(), ref.getTotLen() - insertL + 1);
404 pos = rspd->simulate(sampler, sid, effL);
405 if (pos < 0) return false;
406 if (dir > 0) pos = ref.getTotLen() - pos - insertL;
408 mateL1 = mld->simulate(sampler, insertL);
409 qual1 = qd->simulate(sampler, mateL1);
410 readseq1 = qpro->simulate(sampler, mateL1, pos, dir, qual1, ref);
412 int m2pos = ref.getTotLen() - pos - insertL;
415 mateL2 = mld->simulate(sampler, insertL);
416 qual2 = qd->simulate(sampler, mateL2);
417 readseq2 = qpro->simulate(sampler, mateL2, m2pos, m2dir, qual2, ref);
420 strout<<rid<<"_"<<dir<<"_"<<sid<<"_"<<pos<<"_"<<insertL;
423 read = PairedEndReadQ(SingleReadQ(name + "/1", readseq1, qual1), SingleReadQ(name + "/2", readseq2, qual2));
428 void PairedEndQModel::finishSimulation() {
431 rspd->finishSimulation();
432 qd->finishSimulation();
433 qpro->finishSimulation();
434 nqpro->finishSimulation();
438 void PairedEndQModel::calcMW() {
439 assert(mld->getMinL() >= seedLen);
441 memset(mw, 0, sizeof(double) * (M + 1));
444 for (int i = 1; i <= M; i++) {
445 RefSeq& ref = refs->getRef(i);
446 int totLen = ref.getTotLen();
447 int fullLen = ref.getFullLen();
448 int end = std::min(fullLen, totLen - gld->getMinL() + 1);
453 //seedPos is fpos here
454 for (int seedPos = 0; seedPos < end; seedPos++)
455 if (ref.getMask(seedPos)) {
456 minL = gld->getMinL();
457 maxL = std::min(gld->getMaxL(), totLen - seedPos);
459 for (int fragLen = minL; fragLen <= maxL; fragLen++) {
460 effL = std::min(fullLen, totLen - fragLen + 1);
461 value += gld->getAdjustedProb(fragLen, totLen) * rspd->getAdjustedProb(pfpos, effL, fullLen);
468 // fprintf(stderr, "Warning: %dth reference sequence is masked for almost all positions!\n", i);
474 #endif /* PAIREDENDQMODEL_H_ */