- char outF[STRLEN];
- FILE *fo;
-
- sprintf(modelF, "%s.model", statName);
- model.write(modelF);
-
- vector<int> tlens;
- vector<double> fpkm, tpm, isopct;
- vector<double> glens, gene_eels, gene_counts, gene_tpm, gene_fpkm;
-
- calcExpressionValues(theta, eel, tpm, fpkm);
-
- //calculate IsoPct, etc.
- isopct.assign(M + 1, 0.0);
- tlens.assign(M + 1, 0);
-
- glens.assign(m, 0.0); gene_eels.assign(m, 0.0);
- gene_counts.assign(m, 0.0); gene_tpm.assign(m, 0.0); gene_fpkm.assign(m, 0.0);
-
- for (int i = 0; i < m; i++) {
- int b = gi.spAt(i), e = gi.spAt(i + 1);
- for (int j = b; j < e; j++) {
- const Transcript& transcript = transcripts.getTranscriptAt(j);
- tlens[j] = transcript.getLength();
-
- gene_counts[i] += counts[j];
- gene_tpm[i] += tpm[j];
- gene_fpkm[i] += fpkm[j];
- }
-
- if (gene_tpm[i] < EPSILON) {
- double frac = 1.0 / (e - b);
- for (int j = b; j < e; j++) {
- glens[i] += tlens[j] * frac;
- gene_eels[i] += eel[j] * frac;
- }
- }
- else {
- for (int j = b; j < e; j++) {
- isopct[j] = tpm[j] / gene_tpm[i];
- glens[i] += tlens[j] * isopct[j];
- gene_eels[i] += eel[j] * isopct[j];
- }
- }
- }
-
- //isoform level results
- sprintf(outF, "%s.iso_res", imdName);
- fo = fopen(outF, "w");
- for (int i = 1; i <= M; i++) {
- const Transcript& transcript = transcripts.getTranscriptAt(i);
- fprintf(fo, "%s%c", transcript.getTranscriptID().c_str(), (i < M ? '\t' : '\n'));
- }
- for (int i = 1; i <= M; i++) {
- const Transcript& transcript = transcripts.getTranscriptAt(i);
- fprintf(fo, "%s%c", transcript.getGeneID().c_str(), (i < M ? '\t' : '\n'));
- }
- for (int i = 1; i <= M; i++)
- fprintf(fo, "%d%c", tlens[i], (i < M ? '\t' : '\n'));
- for (int i = 1; i <= M; i++)
- fprintf(fo, "%.2f%c", eel[i], (i < M ? '\t' : '\n'));
- for (int i = 1; i <= M; i++)
- fprintf(fo, "%.2f%c", counts[i], (i < M ? '\t' : '\n'));
- for (int i = 1; i <= M; i++)
- fprintf(fo, "%.2f%c", tpm[i], (i < M ? '\t' : '\n'));
- for (int i = 1; i <= M; i++)
- fprintf(fo, "%.2f%c", fpkm[i], (i < M ? '\t' : '\n'));
- for (int i = 1; i <= M; i++)
- fprintf(fo, "%.2f%c", isopct[i] * 1e2, (i < M ? '\t' : '\n'));
- fclose(fo);
-
- //gene level results
- sprintf(outF, "%s.gene_res", imdName);
- fo = fopen(outF, "w");
- for (int i = 0; i < m; i++) {
- const Transcript& transcript = transcripts.getTranscriptAt(gi.spAt(i));
- fprintf(fo, "%s%c", transcript.getGeneID().c_str(), (i < m - 1 ? '\t' : '\n'));
- }
- for (int i = 0; i < m; i++) {
- int b = gi.spAt(i), e = gi.spAt(i + 1);
- for (int j = b; j < e; j++) {
- fprintf(fo, "%s%c", transcripts.getTranscriptAt(j).getTranscriptID().c_str(), (j < e - 1 ? ',' : (i < m - 1 ? '\t' :'\n')));
- }
- }
- for (int i = 0; i < m; i++)
- fprintf(fo, "%.2f%c", glens[i], (i < m - 1 ? '\t' : '\n'));
- for (int i = 0; i < m; i++)
- fprintf(fo, "%.2f%c", gene_eels[i], (i < m - 1 ? '\t' : '\n'));
- for (int i = 0; i < m; i++)
- fprintf(fo, "%.2f%c", gene_counts[i], (i < m - 1 ? '\t' : '\n'));
- for (int i = 0; i < m; i++)
- fprintf(fo, "%.2f%c", gene_tpm[i], (i < m - 1 ? '\t' : '\n'));
- for (int i = 0; i < m; i++)
- fprintf(fo, "%.2f%c", gene_fpkm[i], (i < m - 1 ? '\t' : '\n'));
- fclose(fo);
-
- if (verbose) { printf("Expression Results are written!\n"); }