- eel.assign(M + 1, 0.0);
- for (int i = 1; i <= M; i++) {
- int totLen = refs.getRef(i).getTotLen();
- int fullLen = refs.getRef(i).getFullLen();
- int pos1 = max(min(totLen - fullLen + 1, ub) - lb, 0);
- int pos2 = max(min(totLen, ub) - lb, 0);
-
- if (pos2 == 0) { eel[i] = 0.0; continue; }
-
- eel[i] = fullLen * cdf[pos1] + ((cdf[pos2] - cdf[pos1]) * (totLen + 1) - (clen[pos2] - clen[pos1]));
- assert(eel[i] >= 0);
- if (eel[i] < MINEEL) { eel[i] = 0.0; }
- }
-
- delete[] pdf;
- delete[] cdf;
- delete[] clen;
-}
-
-template<class ModelType>
-void writeEstimatedParameters(char* modelF, char* imdName) {
- ModelType model;
- double denom;
- char outF[STRLEN];
- FILE *fo;
-
- model.read(modelF);
-
- calcExpectedEffectiveLengths<ModelType>(model);
-
- denom = pme_theta[0];
- for (int i = 1; i <= M; i++)
- if (eel[i] < EPSILON) pme_theta[i] = 0.0;
- else denom += pme_theta[i];
-
- general_assert(denom >= EPSILON, "No Expected Effective Length is no less than " + ftos(MINEEL, 6) + "?!");
-
- for (int i = 0; i <= M; i++) pme_theta[i] /= denom;
-
- denom = 0.0;
- double *mw = model.getMW();
- for (int i = 0; i <= M; i++) {
- pme_theta[i] = (mw[i] < EPSILON ? 0.0 : pme_theta[i] / mw[i]);
- denom += pme_theta[i];
- }
- assert(denom >= EPSILON);
- for (int i = 0; i <= M; i++) pme_theta[i] /= denom;
-
- //calculate tau values
- double *tau = new double[M + 1];
- memset(tau, 0, sizeof(double) * (M + 1));
-
- denom = 0.0;
- for (int i = 1; i <= M; i++)
- if (eel[i] > EPSILON) {
- tau[i] = pme_theta[i] / eel[i];
- denom += tau[i];
- }
-
- general_assert(denom >= EPSILON, "No alignable reads?!");
-
- for (int i = 1; i <= M; i++) {
- tau[i] /= denom;
- }
-
- //isoform level results
- sprintf(outF, "%s.iso_res", imdName);
- fo = fopen(outF, "a");
- general_assert(fo != NULL, "Cannot open " + cstrtos(outF) + "!");
-
- for (int i = 1; i <= M; i++)
- fprintf(fo, "%.2f%c", pme_c[i], (i < M ? '\t' : '\n'));
- for (int i = 1; i <= M; i++)
- fprintf(fo, "%.15g%c", tau[i], (i < M ? '\t' : '\n'));
-
- fclose(fo);
-
- //gene level results
- sprintf(outF, "%s.gene_res", imdName);
- fo = fopen(outF, "a");
- general_assert(fo != NULL, "Cannot open " + cstrtos(outF) + "!");
-
- for (int i = 0; i < m; i++) {
- double sumC = 0.0; // sum of pme counts
- int b = gi.spAt(i), e = gi.spAt(i + 1);
- for (int j = b; j < e; j++) {
- sumC += pme_c[j];
- }
- fprintf(fo, "%.15g%c", sumC, (i < m - 1 ? '\t' : '\n'));
- }