1 #ifndef WRITERESULTS_H_
2 #define WRITERESULTS_H_
12 #include "my_assert.h"
13 #include "GroupInfo.h"
14 #include "Transcript.h"
15 #include "Transcripts.h"
19 #include "SingleModel.h"
20 #include "SingleQModel.h"
21 #include "PairedEndModel.h"
22 #include "PairedEndQModel.h"
24 template<class ModelType>
25 void calcExpectedEffectiveLengths(int M, Refs& refs, ModelType& model, std::vector<double>& eel) {
27 double *pdf = NULL, *cdf = NULL, *clen = NULL; // clen[i] = sigma_{j=1}^{i}pdf[i]*(lb+i)
29 model.getGLD().copyTo(pdf, cdf, lb, ub, span);
30 clen = new double[span + 1];
32 for (int i = 1; i <= span; i++) {
33 clen[i] = clen[i - 1] + pdf[i] * (lb + i);
36 eel.assign(M + 1, 0.0);
37 for (int i = 1; i <= M; i++) {
38 int totLen = refs.getRef(i).getTotLen();
39 int fullLen = refs.getRef(i).getFullLen();
40 int pos1 = std::max(std::min(totLen - fullLen + 1, ub) - lb, 0);
41 int pos2 = std::max(std::min(totLen, ub) - lb, 0);
43 if (pos2 == 0) { eel[i] = 0.0; continue; }
45 eel[i] = fullLen * cdf[pos1] + ((cdf[pos2] - cdf[pos1]) * (totLen + 1) - (clen[pos2] - clen[pos1]));
47 if (eel[i] < MINEEL) { eel[i] = 0.0; }
55 void polishTheta(int M, std::vector<double>& theta, const std::vector<double>& eel, const double* mw) {
58 /* The reason that for noise gene, mw value is 1 is :
59 * currently, all masked positions are for poly(A) sites, which in theory should be filtered out.
60 * So the theta0 does not containing reads from any masked position
63 for (int i = 0; i <= M; i++) {
65 if (i > 0 && (mw[i] < EPSILON || eel[i] < EPSILON)) {
69 theta[i] = theta[i] / mw[i];
72 // currently is OK, since no transcript should be masked totally, only the poly(A) tail related part will be masked
73 general_assert(sum >= EPSILON, "No effective length is no less than" + ftos(MINEEL, 6) + " !");
74 for (int i = 0; i <= M; i++) theta[i] /= sum;
77 void calcExpressionValues(int M, const std::vector<double>& theta, const std::vector<double>& eel, std::vector<double>& tpm, std::vector<double>& fpkm) {
79 std::vector<double> frac;
81 //calculate fraction of count over all mappabile reads
83 frac.assign(M + 1, 0.0);
84 for (int i = 1; i <= M; i++)
85 if (eel[i] >= EPSILON) {
89 general_assert(denom >= EPSILON, "No alignable reads?!");
90 for (int i = 1; i <= M; i++) frac[i] /= denom;
93 fpkm.assign(M + 1, 0.0);
94 for (int i = 1; i <= M; i++)
95 if (eel[i] >= EPSILON) fpkm[i] = frac[i] * 1e9 / eel[i];
98 tpm.assign(M + 1, 0.0);
100 for (int i = 1; i <= M; i++) denom += fpkm[i];
101 for (int i = 1; i <= M; i++) tpm[i] = fpkm[i] / denom * 1e6;
104 inline bool isAlleleSpecific(const char* refName, GroupInfo* gt = NULL, GroupInfo* ta = NULL) {
106 char gtF[STRLEN], taF[STRLEN];
108 sprintf(gtF, "%s.gt", refName);
109 sprintf(taF, "%s.ta", refName);
110 std::ifstream gtIF(gtF), taIF(taF);
111 alleleS = gtIF.is_open() && taIF.is_open();
112 if (gtIF.is_open()) gtIF.close();
113 if (taIF.is_open()) taIF.close();
116 if (gt != NULL) gt->load(gtF);
117 if (ta != NULL) ta->load(taF);
123 void writeResultsEM(int M, const char* refName, const char* imdName, Transcripts& transcripts, std::vector<double>& theta, std::vector<double>& eel, double* counts) {
131 std::vector<int> tlens;
132 std::vector<double> fpkm, tpm, isopct;
133 std::vector<double> glens, gene_eels, gene_counts, gene_tpm, gene_fpkm;
136 sprintf(groupF, "%s.grp", refName);
140 // For allele-specific expression
143 std::vector<double> trans_lens, trans_eels, trans_counts, trans_tpm, trans_fpkm, ta_pct, gt_pct;
145 bool alleleS = isAlleleSpecific(refName, >, &ta); // if allele-specific
147 calcExpressionValues(M, theta, eel, tpm, fpkm);
149 //calculate IsoPct, etc.
150 isopct.assign(M + 1, 0.0);
151 tlens.assign(M + 1, 0);
153 glens.assign(m, 0.0); gene_eels.assign(m, 0.0);
154 gene_counts.assign(m, 0.0); gene_tpm.assign(m, 0.0); gene_fpkm.assign(m, 0.0);
156 for (int i = 0; i < m; i++) {
157 int b = gi.spAt(i), e = gi.spAt(i + 1);
158 for (int j = b; j < e; j++) {
159 const Transcript& transcript = transcripts.getTranscriptAt(j);
160 tlens[j] = transcript.getLength();
162 gene_counts[i] += counts[j];
163 gene_tpm[i] += tpm[j];
164 gene_fpkm[i] += fpkm[j];
167 if (gene_tpm[i] < EPSILON) {
168 double frac = 1.0 / (e - b);
169 for (int j = b; j < e; j++) {
170 glens[i] += tlens[j] * frac;
171 gene_eels[i] += eel[j] * frac;
175 for (int j = b; j < e; j++) {
176 isopct[j] = tpm[j] / gene_tpm[i];
177 glens[i] += tlens[j] * isopct[j];
178 gene_eels[i] += eel[j] * isopct[j];
185 ta_pct.assign(M + 1, 0.0);
186 trans_lens.assign(m_trans, 0.0); trans_eels.assign(m_trans, 0.0);
187 trans_counts.assign(m_trans, 0.0); trans_tpm.assign(m_trans, 0.0); trans_fpkm.assign(m_trans, 0.0);
189 for (int i = 0; i < m_trans; i++) {
190 int b = ta.spAt(i), e = ta.spAt(i + 1);
191 for (int j = b; j < e; j++) {
192 trans_counts[i] += counts[j];
193 trans_tpm[i] += tpm[j];
194 trans_fpkm[i] += fpkm[j];
197 if (trans_tpm[i] < EPSILON) {
198 double frac = 1.0 / (e - b);
199 for (int j = b; j < e; j++) {
200 trans_lens[i] += tlens[j] * frac;
201 trans_eels[i] += eel[j] * frac;
205 for (int j = b; j < e; j++) {
206 ta_pct[j] = tpm[j] / trans_tpm[i];
207 trans_lens[i] += tlens[j] * ta_pct[j];
208 trans_eels[i] += eel[j] * ta_pct[j];
213 gt_pct.assign(m_trans, 0.0);
214 for (int i = 0; i < m; i++)
215 if (gene_tpm[i] >= EPSILON) {
216 int b = gt.spAt(i), e = gt.spAt(i + 1);
217 for (int j = b; j < e; j++) gt_pct[j] = trans_tpm[j] / gene_tpm[i];
222 //isoform level results
223 sprintf(outF, "%s.iso_res", imdName);
224 fo = fopen(outF, "w");
225 for (int i = 1; i <= M; i++) {
226 const Transcript& transcript = transcripts.getTranscriptAt(i);
227 fprintf(fo, "%s%c", transcript.getTranscriptID().c_str(), (i < M ? '\t' : '\n'));
229 for (int i = 1; i <= M; i++) {
230 const Transcript& transcript = transcripts.getTranscriptAt(i);
231 fprintf(fo, "%s%c", transcript.getGeneID().c_str(), (i < M ? '\t' : '\n'));
233 for (int i = 1; i <= M; i++)
234 fprintf(fo, "%d%c", tlens[i], (i < M ? '\t' : '\n'));
235 for (int i = 1; i <= M; i++)
236 fprintf(fo, "%.2f%c", eel[i], (i < M ? '\t' : '\n'));
237 for (int i = 1; i <= M; i++)
238 fprintf(fo, "%.2f%c", counts[i], (i < M ? '\t' : '\n'));
239 for (int i = 1; i <= M; i++)
240 fprintf(fo, "%.2f%c", tpm[i], (i < M ? '\t' : '\n'));
241 for (int i = 1; i <= M; i++)
242 fprintf(fo, "%.2f%c", fpkm[i], (i < M ? '\t' : '\n'));
243 for (int i = 1; i <= M; i++)
244 fprintf(fo, "%.2f%c", isopct[i] * 1e2, (i < M ? '\t' : '\n'));
248 // allele level results
249 sprintf(outF, "%s.allele_res", imdName);
250 fo = fopen(outF, "w");
251 for (int i = 1; i <= M; i++) {
252 const Transcript& transcript = transcripts.getTranscriptAt(i);
253 fprintf(fo, "%s%c", transcript.getSeqName().c_str(), (i < M ? '\t' : '\n'));
255 for (int i = 1; i <= M; i++) {
256 const Transcript& transcript = transcripts.getTranscriptAt(i);
257 fprintf(fo, "%s%c", transcript.getTranscriptID().c_str(), (i < M ? '\t' : '\n'));
259 for (int i = 1; i <= M; i++) {
260 const Transcript& transcript = transcripts.getTranscriptAt(i);
261 fprintf(fo, "%s%c", transcript.getGeneID().c_str(), (i < M ? '\t' : '\n'));
263 for (int i = 1; i <= M; i++)
264 fprintf(fo, "%d%c", tlens[i], (i < M ? '\t' : '\n'));
265 for (int i = 1; i <= M; i++)
266 fprintf(fo, "%.2f%c", eel[i], (i < M ? '\t' : '\n'));
267 for (int i = 1; i <= M; i++)
268 fprintf(fo, "%.2f%c", counts[i], (i < M ? '\t' : '\n'));
269 for (int i = 1; i <= M; i++)
270 fprintf(fo, "%.2f%c", tpm[i], (i < M ? '\t' : '\n'));
271 for (int i = 1; i <= M; i++)
272 fprintf(fo, "%.2f%c", fpkm[i], (i < M ? '\t' : '\n'));
273 for (int i = 1; i <= M; i++)
274 fprintf(fo, "%.2f%c", ta_pct[i] * 1e2, (i < M ? '\t' : '\n'));
275 for (int i = 1; i <= M; i++)
276 fprintf(fo, "%.2f%c", isopct[i] * 1e2, (i < M ? '\t' : '\n'));
279 // isoform level results
280 sprintf(outF, "%s.iso_res", imdName);
281 fo = fopen(outF, "w");
282 for (int i = 0; i < m_trans; i++) {
283 const Transcript& transcript = transcripts.getTranscriptAt(ta.spAt(i));
284 fprintf(fo, "%s%c", transcript.getTranscriptID().c_str(), (i < m_trans - 1 ? '\t' : '\n'));
286 for (int i = 0; i < m_trans; i++) {
287 const Transcript& transcript = transcripts.getTranscriptAt(ta.spAt(i));
288 fprintf(fo, "%s%c", transcript.getGeneID().c_str(), (i < m_trans - 1 ? '\t' : '\n'));
290 for (int i = 0; i < m_trans; i++)
291 fprintf(fo, "%.2f%c", trans_lens[i], (i < m_trans - 1 ? '\t' : '\n'));
292 for (int i = 0; i < m_trans; i++)
293 fprintf(fo, "%.2f%c", trans_eels[i], (i < m_trans - 1 ? '\t' : '\n'));
294 for (int i = 0; i < m_trans; i++)
295 fprintf(fo, "%.2f%c", trans_counts[i], (i < m_trans - 1 ? '\t' : '\n'));
296 for (int i = 0; i < m_trans; i++)
297 fprintf(fo, "%.2f%c", trans_tpm[i], (i < m_trans - 1 ? '\t' : '\n'));
298 for (int i = 0; i < m_trans; i++)
299 fprintf(fo, "%.2f%c", trans_fpkm[i], (i < m_trans - 1 ? '\t' : '\n'));
300 for (int i = 0; i < m_trans; i++)
301 fprintf(fo, "%.2f%c", gt_pct[i] * 1e2, (i < m_trans - 1 ? '\t' : '\n'));
306 sprintf(outF, "%s.gene_res", imdName);
307 fo = fopen(outF, "w");
308 for (int i = 0; i < m; i++) {
309 const Transcript& transcript = transcripts.getTranscriptAt(gi.spAt(i));
310 fprintf(fo, "%s%c", transcript.getGeneID().c_str(), (i < m - 1 ? '\t' : '\n'));
312 for (int i = 0; i < m; i++) {
313 int b = gi.spAt(i), e = gi.spAt(i + 1);
314 std::string curtid = "", tid;
315 for (int j = b; j < e; j++) {
316 tid = transcripts.getTranscriptAt(j).getTranscriptID();
318 if (curtid != "") fprintf(fo, ",");
319 fprintf(fo, "%s", tid.c_str());
323 fprintf(fo, "%c", (i < m - 1 ? '\t' : '\n'));
325 for (int i = 0; i < m; i++)
326 fprintf(fo, "%.2f%c", glens[i], (i < m - 1 ? '\t' : '\n'));
327 for (int i = 0; i < m; i++)
328 fprintf(fo, "%.2f%c", gene_eels[i], (i < m - 1 ? '\t' : '\n'));
329 for (int i = 0; i < m; i++)
330 fprintf(fo, "%.2f%c", gene_counts[i], (i < m - 1 ? '\t' : '\n'));
331 for (int i = 0; i < m; i++)
332 fprintf(fo, "%.2f%c", gene_tpm[i], (i < m - 1 ? '\t' : '\n'));
333 for (int i = 0; i < m; i++)
334 fprintf(fo, "%.2f%c", gene_fpkm[i], (i < m - 1 ? '\t' : '\n'));
337 if (verbose) { printf("Expression Results are written!\n"); }
340 void writeResultsGibbs(int M, int m, int m_trans, GroupInfo& gi, GroupInfo >, GroupInfo &ta, bool alleleS, char* imdName, std::vector<double>& pme_c, std::vector<double>& pme_fpkm, std::vector<double>& pme_tpm, std::vector<double>& pve_c, std::vector<double>& pve_c_genes, std::vector<double>& pve_c_trans) {
344 std::vector<double> isopct;
345 std::vector<double> gene_counts, gene_tpm, gene_fpkm;
347 // For allele-specific expression
348 std::vector<double> trans_counts, trans_tpm, trans_fpkm, ta_pct, gt_pct;
350 //calculate IsoPct, etc.
351 isopct.assign(M + 1, 0.0);
352 gene_counts.assign(m, 0.0); gene_tpm.assign(m, 0.0); gene_fpkm.assign(m, 0.0);
354 for (int i = 0; i < m; i++) {
355 int b = gi.spAt(i), e = gi.spAt(i + 1);
356 for (int j = b; j < e; j++) {
357 gene_counts[i] += pme_c[j];
358 gene_tpm[i] += pme_tpm[j];
359 gene_fpkm[i] += pme_fpkm[j];
361 if (gene_tpm[i] < EPSILON) continue;
362 for (int j = b; j < e; j++)
363 isopct[j] = pme_tpm[j] / gene_tpm[i];
367 ta_pct.assign(M + 1, 0.0);
368 trans_counts.assign(m_trans, 0.0); trans_tpm.assign(m_trans, 0.0); trans_fpkm.assign(m_trans, 0.0);
370 for (int i = 0; i < m_trans; i++) {
371 int b = ta.spAt(i), e = ta.spAt(i + 1);
372 for (int j = b; j < e; j++) {
373 trans_counts[i] += pme_c[j];
374 trans_tpm[i] += pme_tpm[j];
375 trans_fpkm[i] += pme_fpkm[j];
377 if (trans_tpm[i] < EPSILON) continue;
378 for (int j = b; j < e; j++)
379 ta_pct[j] = pme_tpm[j] / trans_tpm[i];
382 gt_pct.assign(m_trans, 0.0);
383 for (int i = 0; i < m; i++)
384 if (gene_tpm[i] >= EPSILON) {
385 int b = gt.spAt(i), e = gt.spAt(i + 1);
386 for (int j = b; j < e; j++) gt_pct[j] = trans_tpm[j] / gene_tpm[i];
391 //isoform level results
392 sprintf(outF, "%s.iso_res", imdName);
393 fo = fopen(outF, "a");
394 general_assert(fo != NULL, "Cannot open " + cstrtos(outF) + "!");
396 for (int i = 1; i <= M; i++)
397 fprintf(fo, "%.2f%c", pme_c[i], (i < M ? '\t' : '\n'));
398 for (int i = 1; i <= M; i++)
399 fprintf(fo, "%.2f%c", sqrt(pve_c[i]), (i < M ? '\t' : '\n'));
400 for (int i = 1; i <= M; i++)
401 fprintf(fo, "%.2f%c", pme_tpm[i], (i < M ? '\t' : '\n'));
402 for (int i = 1; i <= M; i++)
403 fprintf(fo, "%.2f%c", pme_fpkm[i], (i < M ? '\t' : '\n'));
404 for (int i = 1; i <= M; i++)
405 fprintf(fo, "%.2f%c", isopct[i] * 1e2, (i < M ? '\t' : '\n'));
409 //allele level results
410 sprintf(outF, "%s.allele_res", imdName);
411 fo = fopen(outF, "a");
412 general_assert(fo != NULL, "Cannot open " + cstrtos(outF) + "!");
414 for (int i = 1; i <= M; i++)
415 fprintf(fo, "%.2f%c", pme_c[i], (i < M ? '\t' : '\n'));
416 for (int i = 1; i <= M; i++)
417 fprintf(fo, "%.2f%c", sqrt(pve_c[i]), (i < M ? '\t' : '\n'));
418 for (int i = 1; i <= M; i++)
419 fprintf(fo, "%.2f%c", pme_tpm[i], (i < M ? '\t' : '\n'));
420 for (int i = 1; i <= M; i++)
421 fprintf(fo, "%.2f%c", pme_fpkm[i], (i < M ? '\t' : '\n'));
422 for (int i = 1; i <= M; i++)
423 fprintf(fo, "%.2f%c", ta_pct[i] * 1e2, (i < M ? '\t' : '\n'));
424 for (int i = 1; i <= M; i++)
425 fprintf(fo, "%.2f%c", isopct[i] * 1e2, (i < M ? '\t' : '\n'));
428 //isoform level results
429 sprintf(outF, "%s.iso_res", imdName);
430 fo = fopen(outF, "a");
431 general_assert(fo != NULL, "Cannot open " + cstrtos(outF) + "!");
433 for (int i = 0; i < m_trans; i++)
434 fprintf(fo, "%.2f%c", trans_counts[i], (i < m_trans - 1 ? '\t' : '\n'));
435 for (int i = 0; i < m_trans; i++)
436 fprintf(fo, "%.2f%c", sqrt(pve_c_trans[i]), (i < m_trans - 1 ? '\t' : '\n'));
437 for (int i = 0; i < m_trans; i++)
438 fprintf(fo, "%.2f%c", trans_tpm[i], (i < m_trans - 1 ? '\t' : '\n'));
439 for (int i = 0; i < m_trans; i++)
440 fprintf(fo, "%.2f%c", trans_fpkm[i], (i < m_trans - 1 ? '\t' : '\n'));
441 for (int i = 0; i < m_trans; i++)
442 fprintf(fo, "%.2f%c", gt_pct[i] * 1e2, (i < m_trans - 1 ? '\t' : '\n'));
447 sprintf(outF, "%s.gene_res", imdName);
448 fo = fopen(outF, "a");
449 general_assert(fo != NULL, "Cannot open " + cstrtos(outF) + "!");
451 for (int i = 0; i < m; i++)
452 fprintf(fo, "%.2f%c", gene_counts[i], (i < m - 1 ? '\t' : '\n'));
453 for (int i = 0; i < m; i++)
454 fprintf(fo, "%.2f%c", sqrt(pve_c_genes[i]), (i < m - 1 ? '\t' : '\n'));
455 for (int i = 0; i < m; i++)
456 fprintf(fo, "%.2f%c", gene_tpm[i], (i < m - 1 ? '\t' : '\n'));
457 for (int i = 0; i < m; i++)
458 fprintf(fo, "%.2f%c", gene_fpkm[i], (i < m - 1 ? '\t' : '\n'));
461 if (verbose) { printf("Gibbs based expression values are written!\n"); }
464 void writeResultsSimulation(int M, char* refName, char* outFN, Transcripts& transcripts, std::vector<double>& eel, std::vector<double>& counts) {
473 sprintf(groupF, "%s.grp", refName);
477 std::vector<int> tlens;
478 std::vector<double> tpm, fpkm, isopct;
479 std::vector<double> glens, gene_eels, gene_counts, gene_tpm, gene_fpkm;
481 // For allele-specific expression
484 std::vector<double> trans_lens, trans_eels, trans_counts, trans_tpm, trans_fpkm, ta_pct, gt_pct;
486 bool alleleS = isAlleleSpecific(refName, >, &ta); // if allele-specific
488 for (int i = 1; i <= M; i++)
489 general_assert(eel[i] > EPSILON || counts[i] <= EPSILON, "An isoform whose effecitve length < " + ftos(MINEEL, 6) + " got sampled!");
491 calcExpressionValues(M, counts, eel, tpm, fpkm);
493 //calculate IsoPct, etc.
494 isopct.assign(M + 1, 0.0);
495 tlens.assign(M + 1, 0);
497 glens.assign(m, 0.0); gene_eels.assign(m, 0.0);
498 gene_counts.assign(m, 0.0); gene_tpm.assign(m, 0.0); gene_fpkm.assign(m, 0.0);
500 for (int i = 0; i < m; i++) {
501 int b = gi.spAt(i), e = gi.spAt(i + 1);
502 for (int j = b; j < e; j++) {
503 const Transcript& transcript = transcripts.getTranscriptAt(j);
504 tlens[j] = transcript.getLength();
506 gene_counts[i] += counts[j];
507 gene_tpm[i] += tpm[j];
508 gene_fpkm[i] += fpkm[j];
511 if (gene_tpm[i] < EPSILON) {
512 double frac = 1.0 / (e - b);
513 for (int j = b; j < e; j++) {
514 glens[i] += tlens[j] * frac;
515 gene_eels[i] += eel[j] * frac;
519 for (int j = b; j < e; j++) {
520 isopct[j] = tpm[j] / gene_tpm[i];
521 glens[i] += tlens[j] * isopct[j];
522 gene_eels[i] += eel[j] * isopct[j];
529 ta_pct.assign(M + 1, 0.0);
530 trans_lens.assign(m_trans, 0.0); trans_eels.assign(m_trans, 0.0);
531 trans_counts.assign(m_trans, 0.0); trans_tpm.assign(m_trans, 0.0); trans_fpkm.assign(m_trans, 0.0);
533 for (int i = 0; i < m_trans; i++) {
534 int b = ta.spAt(i), e = ta.spAt(i + 1);
535 for (int j = b; j < e; j++) {
536 trans_counts[i] += counts[j];
537 trans_tpm[i] += tpm[j];
538 trans_fpkm[i] += fpkm[j];
541 if (trans_tpm[i] < EPSILON) {
542 double frac = 1.0 / (e - b);
543 for (int j = b; j < e; j++) {
544 trans_lens[i] += tlens[j] * frac;
545 trans_eels[i] += eel[j] * frac;
549 for (int j = b; j < e; j++) {
550 ta_pct[j] = tpm[j] / trans_tpm[i];
551 trans_lens[i] += tlens[j] * ta_pct[j];
552 trans_eels[i] += eel[j] * ta_pct[j];
557 gt_pct.assign(m_trans, 0.0);
558 for (int i = 0; i < m; i++)
559 if (gene_tpm[i] >= EPSILON) {
560 int b = gt.spAt(i), e = gt.spAt(i + 1);
561 for (int j = b; j < e; j++) gt_pct[j] = trans_tpm[j] / gene_tpm[i];
567 sprintf(outF, "%s.sim.alleles.results", outFN);
568 fo = fopen(outF, "w");
569 fprintf(fo, "allele_id\ttranscript_id\tgene_id\tlength\teffective_length\tcount\tTPM\tFPKM\tAlleleIsoPct\tAlleleGenePct\n");
570 for (int i = 1; i <= M; i++) {
571 const Transcript& transcript = transcripts.getTranscriptAt(i);
572 fprintf(fo, "%s\t%s\t%s\t%d\t%.2f\t%.2f\t%.2f\t%.2f\t%.2f\t%.2f\n", transcript.getSeqName().c_str(), transcript.getTranscriptID().c_str(), transcript.getGeneID().c_str(), tlens[i],
573 eel[i], counts[i], tpm[i], fpkm[i], ta_pct[i] * 1e2, isopct[i] * 1e2);
579 sprintf(outF, "%s.sim.isoforms.results", outFN);
580 fo = fopen(outF, "w");
581 fprintf(fo, "transcript_id\tgene_id\tlength\teffective_length\tcount\tTPM\tFPKM\tIsoPct\n");
583 for (int i = 1; i <= M; i++) {
584 const Transcript& transcript = transcripts.getTranscriptAt(i);
585 fprintf(fo, "%s\t%s\t%d\t%.2f\t%.2f\t%.2f\t%.2f\t%.2f\n", transcript.getTranscriptID().c_str(), transcript.getGeneID().c_str(), tlens[i],
586 eel[i], counts[i], tpm[i], fpkm[i], isopct[i] * 1e2);
590 for (int i = 0; i < m_trans; i++) {
591 const Transcript& transcript = transcripts.getTranscriptAt(ta.spAt(i));
592 fprintf(fo, "%s\t%s\t%.2f\t%.2f\t%.2f\t%.2f\t%.2f\t%.2f\n", transcript.getTranscriptID().c_str(), transcript.getGeneID().c_str(), trans_lens[i],
593 trans_eels[i], trans_counts[i], trans_tpm[i], trans_fpkm[i], gt_pct[i] * 1e2);
599 sprintf(outF, "%s.sim.genes.results", outFN);
600 fo = fopen(outF, "w");
601 fprintf(fo, "gene_id\ttranscript_id(s)\tlength\teffective_length\tcount\tTPM\tFPKM\n");
602 for (int i = 0; i < m; i++) {
603 int b = gi.spAt(i), e = gi.spAt(i + 1);
604 const std::string& gene_id = transcripts.getTranscriptAt(b).getGeneID();
605 fprintf(fo, "%s\t", gene_id.c_str());
606 std::string curtid = "", tid;
607 for (int j = b; j < e; j++) {
608 tid = transcripts.getTranscriptAt(j).getTranscriptID();
610 if (curtid != "") fprintf(fo, ",");
611 fprintf(fo, "%s", tid.c_str());
615 fprintf(fo, "\t%.2f\t%.2f\t%.2f\t%.2f\t%.2f\n", glens[i], gene_eels[i], gene_counts[i], gene_tpm[i], gene_fpkm[i]);