13 #include "Orientation.h"
17 #include "NoiseProfile.h"
19 #include "ModelParams.h"
22 #include "SingleRead.h"
23 #include "SingleHit.h"
24 #include "ReadReader.h"
30 SingleModel(Refs* refs = NULL) {
32 M = (refs != NULL ? refs->getM() : 0);
33 memset(N, 0, sizeof(N));
35 needCalcConPrb = true;
37 ori = new Orientation();
40 rspd = new RSPD(estRSPD);
42 npro = new NoiseProfile();
44 mean = -1.0; sd = 0.0;
50 //If it is not a master node, only init & update can be used!
51 SingleModel(ModelParams& params, bool isMaster = true) {
53 memcpy(N, params.N, sizeof(params.N));
55 estRSPD = params.estRSPD;
56 mean = params.mean; sd = params.sd;
57 seedLen = params.seedLen;
58 needCalcConPrb = true;
60 ori = NULL; gld = NULL; mld = NULL; rspd = NULL; pro = NULL; npro = NULL;
64 gld = new LenDist(params.minL, params.maxL);
65 if (mean >= EPSILON) {
66 mld = new LenDist(params.mate_minL, params.mate_maxL);
68 if (!estRSPD) { rspd = new RSPD(estRSPD); }
71 ori = new Orientation(params.probF);
72 if (estRSPD) { rspd = new RSPD(estRSPD, params.B); }
73 pro = new Profile(params.maxL);
74 npro = new NoiseProfile();
79 if (ori != NULL) delete ori;
80 if (gld != NULL) delete gld;
81 if (mld != NULL) delete mld;
82 if (rspd != NULL) delete rspd;
83 if (pro != NULL) delete pro;
84 if (npro != NULL) delete npro;
85 if (mw != NULL) delete[] mw;
89 void estimateFromReads(const char*);
91 //if prob is too small, just make it 0
92 double getConPrb(const SingleRead& read, const SingleHit& hit) {
93 if (read.isLowQuality()) return 0.0;
96 int sid = hit.getSid();
97 RefSeq &ref = refs->getRef(sid);
98 int fullLen = ref.getFullLen();
99 int totLen = ref.getTotLen();
100 int dir = hit.getDir();
101 int pos = hit.getPos();
102 int readLen = read.getReadLength();
103 int fpos = (dir == 0 ? pos : totLen - pos - readLen); // the aligned position reported in SAM file, should be a coordinate in forward strand
105 assert(fpos >= 0 && fpos + readLen <= totLen && readLen <= totLen);
106 int seedPos = (dir == 0 ? pos : totLen - pos - seedLen); // the aligned position of the seed in forward strand coordinates
107 if (seedPos >= fullLen || ref.getMask(seedPos)) return 0.0;
113 int minL = std::max(readLen, gld->getMinL());
114 int maxL = std::min(totLen - pos, gld->getMaxL());
115 int pfpos; // possible fpos for fragment
117 for (int fragLen = minL; fragLen <= maxL; fragLen++) {
118 pfpos = (dir == 0 ? pos : totLen - pos - fragLen);
119 effL = std::min(fullLen, totLen - fragLen + 1);
120 value += gld->getAdjustedProb(fragLen, totLen) * rspd->getAdjustedProb(pfpos, effL, fullLen) * mld->getAdjustedProb(readLen, fragLen);
124 effL = std::min(fullLen, totLen - readLen + 1);
125 value = gld->getAdjustedProb(readLen, totLen) * rspd->getAdjustedProb(fpos, effL, fullLen);
128 prob = ori->getProb(dir) * value * pro->getProb(read.getReadSeq(), ref, pos, dir);
130 if (prob < EPSILON) { prob = 0.0; }
133 prob = (mw[sid] < EPSILON ? 0.0 : prob / mw[sid]);
138 double getNoiseConPrb(const SingleRead& read) {
139 if (read.isLowQuality()) return 0.0;
140 double prob = mld != NULL ? mld->getProb(read.getReadLength()) : gld->getProb(read.getReadLength());
141 prob *= npro->getProb(read.getReadSeq());
142 if (prob < EPSILON) { prob = 0.0; }
144 prob = (mw[0] < EPSILON ? 0.0 : prob / mw[0]);
149 double getLogP() { return npro->getLogP(); }
153 void update(const SingleRead& read, const SingleHit& hit, double frac) {
154 if (read.isLowQuality() || frac < EPSILON) return;
156 RefSeq& ref = refs->getRef(hit.getSid());
157 int dir = hit.getDir();
158 int pos = hit.getPos();
161 int fullLen = ref.getFullLen();
163 // Only use one strand to estimate RSPD
164 if (ori->getProb(0) >= ORIVALVE && dir == 0) {
165 rspd->update(pos, fullLen, frac);
168 if (ori->getProb(0) < ORIVALVE && dir == 1) {
169 int totLen = ref.getTotLen();
170 int readLen = read.getReadLength();
175 int minL = std::max(readLen, gld->getMinL());
176 int maxL = std::min(totLen - pos, gld->getMaxL());
178 assert(maxL >= minL);
179 std::vector<double> frag_vec(maxL - minL + 1, 0.0);
181 for (int fragLen = minL; fragLen <= maxL; fragLen++) {
182 pfpos = totLen - pos - fragLen;
183 effL = std::min(fullLen, totLen - fragLen + 1);
184 frag_vec[fragLen - minL] = gld->getAdjustedProb(fragLen, totLen) * rspd->getAdjustedProb(pfpos, effL, fullLen) * mld->getAdjustedProb(readLen, fragLen);
185 sum += frag_vec[fragLen - minL];
187 assert(sum >= EPSILON);
188 for (int fragLen = minL; fragLen <= maxL; fragLen++) {
189 pfpos = totLen - pos - fragLen;
190 rspd->update(pfpos, fullLen, frac * (frag_vec[fragLen - minL] / sum));
194 rspd->update(totLen - pos - readLen, fullLen, frac);
198 pro->update(read.getReadSeq(), ref, pos, dir, frac);
201 void updateNoise(const SingleRead& read, double frac) {
202 if (read.isLowQuality() || frac < EPSILON) return;
204 npro->update(read.getReadSeq(), frac);
209 void collect(const SingleModel&);
211 bool getNeedCalcConPrb() { return needCalcConPrb; }
212 void setNeedCalcConPrb(bool value) { needCalcConPrb = value; }
216 //double* getP1() { return p1; }
217 //double* getP2() { return p2; }
219 void read(const char*);
220 void write(const char*);
222 const LenDist& getGLD() { return *gld; }
224 void startSimulation(simul*, double*);
225 bool simulate(int, SingleRead&, int&);
226 void finishSimulation();
233 int getModelType() const { return model_type; }
236 static const int model_type = 0;
237 static const int read_type = 0;
244 //double *p1, *p2; P_i' & P_i''
246 bool estRSPD; // true if estimate RSPD
247 bool needCalcConPrb; // true need, false does not need
255 simul *sampler; // for simulation
256 double *theta_cdf; // for simulation
258 double *mw; // for masking
263 void SingleModel::estimateFromReads(const char* readFN) {
265 char readFs[2][STRLEN];
268 mld != NULL ? mld->init() : gld->init();
269 for (int i = 0; i < 3; i++)
271 genReadFileNames(readFN, i, read_type, s, readFs);
272 ReadReader<SingleRead> reader(s, readFs);
275 while (reader.next(read)) {
276 mld != NULL ? mld->update(read.getReadLength(), 1.0) : gld->update(read.getReadLength(), 1.0);
277 if (i == 0) { npro->updateC(read.getReadSeq()); }
280 if (verbose && cnt % 1000000 == 0) { printf("%d READS PROCESSED\n", cnt); }
283 if (verbose) { printf("estimateFromReads, N%d finished.\n", i); }
286 mld != NULL ? mld->finish() : gld->finish();
288 if (mean >= EPSILON) {
289 assert(mld->getMaxL() <= gld->getMaxL());
290 gld->setAsNormal(mean, sd, std::max(mld->getMinL(), gld->getMinL()), gld->getMaxL());
292 npro->calcInitParams();
294 mw = new double[M + 1];
298 void SingleModel::init() {
299 if (estRSPD) rspd->init();
304 void SingleModel::finish() {
305 if (estRSPD) rspd->finish();
308 needCalcConPrb = true;
309 if (estRSPD) calcMW();
312 void SingleModel::collect(const SingleModel& o) {
313 if (estRSPD) rspd->collect(*(o.rspd));
314 pro->collect(*(o.pro));
315 npro->collect(*(o.npro));
318 //Only master node can call
319 void SingleModel::read(const char* inpF) {
321 FILE *fi = fopen(inpF, "r");
322 if (fi == NULL) { fprintf(stderr, "Cannot open %s! It may not exist.\n", inpF); exit(-1); }
324 fscanf(fi, "%d", &val);
325 assert(val == model_type);
329 fscanf(fi, "%d", &val);
331 if (mld == NULL) mld = new LenDist();
338 if (fscanf(fi, "%d", &val) == 1) {
341 mw = new double[M + 1];
342 for (int i = 0; i <= M; i++) fscanf(fi, "%lf", &mw[i]);
349 //Only master node can call. Only be called at EM.cpp
350 void SingleModel::write(const char* outF) {
351 FILE *fo = fopen(outF, "w");
353 fprintf(fo, "%d\n", model_type);
356 ori->write(fo); fprintf(fo, "\n");
357 gld->write(fo); fprintf(fo, "\n");
362 else { fprintf(fo, "0\n"); }
364 rspd->write(fo); fprintf(fo, "\n");
365 pro->write(fo); fprintf(fo, "\n");
369 fprintf(fo, "\n%d\n", M);
370 for (int i = 0; i < M; i++) {
371 fprintf(fo, "%.15g ", mw[i]);
373 fprintf(fo, "%.15g\n", mw[M]);
379 void SingleModel::startSimulation(simul* sampler, double* theta) {
380 this->sampler = sampler;
382 theta_cdf = new double[M + 1];
383 for (int i = 0; i <= M; i++) {
384 theta_cdf[i] = theta[i];
385 if (i > 0) theta_cdf[i] += theta_cdf[i - 1];
388 rspd->startSimulation(M, refs);
389 pro->startSimulation();
390 npro->startSimulation();
393 bool SingleModel::simulate(int rid, SingleRead& read, int& sid) {
394 int dir, pos, readLen, fragLen;
397 std::ostringstream strout;
399 sid = sampler->sample(theta_cdf, M + 1);
403 readLen = (mld != NULL ? mld->simulate(sampler, -1) : gld->simulate(sampler, -1));
404 readseq = npro->simulate(sampler, readLen);
407 RefSeq &ref = refs->getRef(sid);
408 dir = ori->simulate(sampler);
409 fragLen = gld->simulate(sampler, ref.getTotLen());
410 if (fragLen < 0) return false;
411 int effL = std::min(ref.getFullLen(), ref.getTotLen() - fragLen + 1);
412 pos = rspd->simulate(sampler, sid, effL);
413 if (pos < 0) return false;
414 if (dir > 0) pos = ref.getTotLen() - pos - fragLen;
417 readLen = mld->simulate(sampler, fragLen);
418 if (readLen < 0) return false;
419 readseq = pro->simulate(sampler, readLen, pos, dir, ref);
422 readseq = pro->simulate(sampler, fragLen, pos, dir, ref);
426 std::ostringstream stdout;
427 stdout<<rid<<"_"<<dir<<"_"<<sid<<"_"<<pos;
430 read = SingleRead(name, readseq);
435 void SingleModel::finishSimulation() {
438 rspd->finishSimulation();
439 pro->finishSimulation();
440 npro->finishSimulation();
443 void SingleModel::calcMW() {
446 assert(seedLen >= OLEN && (mld == NULL ? gld->getMinL() : mld->getMinL()) >= seedLen);
448 memset(mw, 0, sizeof(double) * (M + 1));
452 probF = ori->getProb(0);
453 probR = ori->getProb(1);
455 for (int i = 1; i <= M; i++) {
456 RefSeq& ref = refs->getRef(i);
457 int totLen = ref.getTotLen();
458 int fullLen = ref.getFullLen();
462 int end = std::min(fullLen, totLen - seedLen + 1);
465 for (int seedPos = 0; seedPos < end; seedPos++)
466 if (ref.getMask(seedPos)) {
468 minL = gld->getMinL();
469 maxL = std::min(gld->getMaxL(), totLen - seedPos);
471 for (int fragLen = minL; fragLen <= maxL; fragLen++) {
472 effL = std::min(fullLen, totLen - fragLen + 1);
473 factor = (mld == NULL ? 1.0 : mld->getAdjustedCumulativeProb(std::min(mld->getMaxL(), fragLen), fragLen));
474 value += probF * gld->getAdjustedProb(fragLen, totLen) * rspd->getAdjustedProb(pfpos, effL, fullLen) * factor;
477 minL = gld->getMinL();
478 maxL = std::min(gld->getMaxL(), seedPos + seedLen);
479 for (int fragLen = minL; fragLen <= maxL; fragLen++) {
480 pfpos = seedPos - (fragLen - seedLen);
481 effL = std::min(fullLen, totLen - fragLen + 1);
482 factor = (mld == NULL ? 1.0 : mld->getAdjustedCumulativeProb(std::min(mld->getMaxL(), fragLen), fragLen));
483 value += probR * gld->getAdjustedProb(fragLen, totLen) * rspd->getAdjustedProb(pfpos, effL, fullLen) * factor;
487 //for reverse strand masking
488 for (int seedPos = end; seedPos <= totLen - seedLen; seedPos++) {
489 minL = std::max(gld->getMinL(), seedPos + seedLen - fullLen + 1);
490 maxL = std::min(gld->getMaxL(), seedPos + seedLen);
491 for (int fragLen = minL; fragLen <= maxL; fragLen++) {
492 pfpos = seedPos - (fragLen - seedLen);
493 effL = std::min(fullLen, totLen - fragLen + 1);
494 factor = (mld == NULL ? 1.0 : mld->getAdjustedCumulativeProb(std::min(mld->getMaxL(), fragLen), fragLen));
495 value += probR * gld->getAdjustedProb(fragLen, totLen) * rspd->getAdjustedProb(pfpos, effL, fullLen) * factor;
502 // fprintf(stderr, "Warning: %dth reference sequence is masked for almost all positions!\n", i);
508 #endif /* SINGLEMODEL_H_ */