1 #ifndef WRITERESULTS_H_
2 #define WRITERESULTS_H_
11 #include "my_assert.h"
12 #include "GroupInfo.h"
13 #include "Transcript.h"
14 #include "Transcripts.h"
18 #include "SingleModel.h"
19 #include "SingleQModel.h"
20 #include "PairedEndModel.h"
21 #include "PairedEndQModel.h"
23 template<class ModelType>
24 void calcExpectedEffectiveLengths(int M, Refs& refs, ModelType& model, std::vector<double>& eel) {
26 double *pdf = NULL, *cdf = NULL, *clen = NULL; // clen[i] = sigma_{j=1}^{i}pdf[i]*(lb+i)
28 model.getGLD().copyTo(pdf, cdf, lb, ub, span);
29 clen = new double[span + 1];
31 for (int i = 1; i <= span; i++) {
32 clen[i] = clen[i - 1] + pdf[i] * (lb + i);
35 eel.assign(M + 1, 0.0);
36 for (int i = 1; i <= M; i++) {
37 int totLen = refs.getRef(i).getTotLen();
38 int fullLen = refs.getRef(i).getFullLen();
39 int pos1 = std::max(std::min(totLen - fullLen + 1, ub) - lb, 0);
40 int pos2 = std::max(std::min(totLen, ub) - lb, 0);
42 if (pos2 == 0) { eel[i] = 0.0; continue; }
44 eel[i] = fullLen * cdf[pos1] + ((cdf[pos2] - cdf[pos1]) * (totLen + 1) - (clen[pos2] - clen[pos1]));
46 if (eel[i] < MINEEL) { eel[i] = 0.0; }
54 void polishTheta(int M, std::vector<double>& theta, const std::vector<double>& eel, const double* mw) {
57 /* The reason that for noise gene, mw value is 1 is :
58 * currently, all masked positions are for poly(A) sites, which in theory should be filtered out.
59 * So the theta0 does not containing reads from any masked position
62 for (int i = 0; i <= M; i++) {
64 if (i > 0 && (mw[i] < EPSILON || eel[i] < EPSILON)) {
68 theta[i] = theta[i] / mw[i];
71 // currently is OK, since no transcript should be masked totally, only the poly(A) tail related part will be masked
72 general_assert(sum >= EPSILON, "No effective length is no less than" + ftos(MINEEL, 6) + " !");
73 for (int i = 0; i <= M; i++) theta[i] /= sum;
76 void calcExpressionValues(int M, const std::vector<double>& theta, const std::vector<double>& eel, std::vector<double>& tpm, std::vector<double>& fpkm) {
78 std::vector<double> frac;
80 //calculate fraction of count over all mappabile reads
82 frac.assign(M + 1, 0.0);
83 for (int i = 1; i <= M; i++)
84 if (eel[i] >= EPSILON) {
88 general_assert(denom >= EPSILON, "No alignable reads?!");
89 for (int i = 1; i <= M; i++) frac[i] /= denom;
92 fpkm.assign(M + 1, 0.0);
93 for (int i = 1; i <= M; i++)
94 if (eel[i] >= EPSILON) fpkm[i] = frac[i] * 1e9 / eel[i];
97 tpm.assign(M + 1, 0.0);
99 for (int i = 1; i <= M; i++) denom += fpkm[i];
100 for (int i = 1; i <= M; i++) tpm[i] = fpkm[i] / denom * 1e6;
103 inline bool isAlleleSpecific(const char* refName, GroupInfo* gt = NULL, GroupInfo* ta = NULL) {
105 char gtF[STRLEN], taF[STRLEN];
107 sprintf(gtF, "%s.gt", refName);
108 sprintf(taF, "%s.ta", refName);
109 std::ifstream gtIF(gtF), taIF(taF);
110 alleleS = gtIF.is_open() && taIF.is_open();
111 if (gtIF.is_open()) gtIF.close();
112 if (taIF.is_open()) taIF.close();
115 if (gt != NULL) gt->load(gtF);
116 if (ta != NULL) ta->load(taF);
122 void writeResultsEM(int M, const char* refName, const char* imdName, Transcripts& transcripts, std::vector<double>& theta, std::vector<double>& eel, double* counts) {
130 std::vector<int> tlens;
131 std::vector<double> fpkm, tpm, isopct;
132 std::vector<double> glens, gene_eels, gene_counts, gene_tpm, gene_fpkm;
135 sprintf(groupF, "%s.grp", refName);
139 // For allele-specific expression
142 std::vector<double> trans_lens, trans_eels, trans_counts, trans_tpm, trans_fpkm, ta_pct, gt_pct;
144 bool alleleS = isAlleleSpecific(refName, >, &ta); // if allele-specific
146 calcExpressionValues(M, theta, eel, tpm, fpkm);
148 //calculate IsoPct, etc.
149 isopct.assign(M + 1, 0.0);
150 tlens.assign(M + 1, 0);
152 glens.assign(m, 0.0); gene_eels.assign(m, 0.0);
153 gene_counts.assign(m, 0.0); gene_tpm.assign(m, 0.0); gene_fpkm.assign(m, 0.0);
155 for (int i = 0; i < m; i++) {
156 int b = gi.spAt(i), e = gi.spAt(i + 1);
157 for (int j = b; j < e; j++) {
158 const Transcript& transcript = transcripts.getTranscriptAt(j);
159 tlens[j] = transcript.getLength();
161 gene_counts[i] += counts[j];
162 gene_tpm[i] += tpm[j];
163 gene_fpkm[i] += fpkm[j];
166 if (gene_tpm[i] < EPSILON) {
167 double frac = 1.0 / (e - b);
168 for (int j = b; j < e; j++) {
169 glens[i] += tlens[j] * frac;
170 gene_eels[i] += eel[j] * frac;
174 for (int j = b; j < e; j++) {
175 isopct[j] = tpm[j] / gene_tpm[i];
176 glens[i] += tlens[j] * isopct[j];
177 gene_eels[i] += eel[j] * isopct[j];
184 ta_pct.assign(M + 1, 0.0);
185 trans_lens.assign(m_trans, 0.0); trans_eels.assign(m_trans, 0.0);
186 trans_counts.assign(m_trans, 0.0); trans_tpm.assign(m_trans, 0.0); trans_fpkm.assign(m_trans, 0.0);
188 for (int i = 0; i < m_trans; i++) {
189 int b = ta.spAt(i), e = ta.spAt(i + 1);
190 for (int j = b; j < e; j++) {
191 trans_counts[i] += counts[j];
192 trans_tpm[i] += tpm[j];
193 trans_fpkm[i] += fpkm[j];
196 if (trans_tpm[i] < EPSILON) {
197 double frac = 1.0 / (e - b);
198 for (int j = b; j < e; j++) {
199 trans_lens[i] += tlens[j] * frac;
200 trans_eels[i] += eel[j] * frac;
204 for (int j = b; j < e; j++) {
205 ta_pct[j] = tpm[j] / trans_tpm[i];
206 trans_lens[i] += tlens[j] * ta_pct[j];
207 trans_eels[i] += eel[j] * ta_pct[j];
212 gt_pct.assign(m_trans, 0.0);
213 for (int i = 0; i < m; i++)
214 if (gene_tpm[i] >= EPSILON) {
215 int b = gt.spAt(i), e = gt.spAt(i + 1);
216 for (int j = b; j < e; j++) gt_pct[j] = trans_tpm[j] / gene_tpm[i];
221 //isoform level results
222 sprintf(outF, "%s.iso_res", imdName);
223 fo = fopen(outF, "w");
224 for (int i = 1; i <= M; i++) {
225 const Transcript& transcript = transcripts.getTranscriptAt(i);
226 fprintf(fo, "%s%c", transcript.getTranscriptID().c_str(), (i < M ? '\t' : '\n'));
228 for (int i = 1; i <= M; i++) {
229 const Transcript& transcript = transcripts.getTranscriptAt(i);
230 fprintf(fo, "%s%c", transcript.getGeneID().c_str(), (i < M ? '\t' : '\n'));
232 for (int i = 1; i <= M; i++)
233 fprintf(fo, "%d%c", tlens[i], (i < M ? '\t' : '\n'));
234 for (int i = 1; i <= M; i++)
235 fprintf(fo, "%.2f%c", eel[i], (i < M ? '\t' : '\n'));
236 for (int i = 1; i <= M; i++)
237 fprintf(fo, "%.2f%c", counts[i], (i < M ? '\t' : '\n'));
238 for (int i = 1; i <= M; i++)
239 fprintf(fo, "%.2f%c", tpm[i], (i < M ? '\t' : '\n'));
240 for (int i = 1; i <= M; i++)
241 fprintf(fo, "%.2f%c", fpkm[i], (i < M ? '\t' : '\n'));
242 for (int i = 1; i <= M; i++)
243 fprintf(fo, "%.2f%c", isopct[i] * 1e2, (i < M ? '\t' : '\n'));
247 // allele level results
248 sprintf(outF, "%s.allele_res", imdName);
249 fo = fopen(outF, "w");
250 for (int i = 1; i <= M; i++) {
251 const Transcript& transcript = transcripts.getTranscriptAt(i);
252 fprintf(fo, "%s%c", transcript.getSeqName().c_str(), (i < M ? '\t' : '\n'));
254 for (int i = 1; i <= M; i++) {
255 const Transcript& transcript = transcripts.getTranscriptAt(i);
256 fprintf(fo, "%s%c", transcript.getTranscriptID().c_str(), (i < M ? '\t' : '\n'));
258 for (int i = 1; i <= M; i++) {
259 const Transcript& transcript = transcripts.getTranscriptAt(i);
260 fprintf(fo, "%s%c", transcript.getGeneID().c_str(), (i < M ? '\t' : '\n'));
262 for (int i = 1; i <= M; i++)
263 fprintf(fo, "%d%c", tlens[i], (i < M ? '\t' : '\n'));
264 for (int i = 1; i <= M; i++)
265 fprintf(fo, "%.2f%c", eel[i], (i < M ? '\t' : '\n'));
266 for (int i = 1; i <= M; i++)
267 fprintf(fo, "%.2f%c", counts[i], (i < M ? '\t' : '\n'));
268 for (int i = 1; i <= M; i++)
269 fprintf(fo, "%.2f%c", tpm[i], (i < M ? '\t' : '\n'));
270 for (int i = 1; i <= M; i++)
271 fprintf(fo, "%.2f%c", fpkm[i], (i < M ? '\t' : '\n'));
272 for (int i = 1; i <= M; i++)
273 fprintf(fo, "%.2f%c", ta_pct[i] * 1e2, (i < M ? '\t' : '\n'));
274 for (int i = 1; i <= M; i++)
275 fprintf(fo, "%.2f%c", isopct[i] * 1e2, (i < M ? '\t' : '\n'));
278 // isoform level results
279 sprintf(outF, "%s.iso_res", imdName);
280 fo = fopen(outF, "w");
281 for (int i = 0; i < m_trans; i++) {
282 const Transcript& transcript = transcripts.getTranscriptAt(ta.spAt(i));
283 fprintf(fo, "%s%c", transcript.getTranscriptID().c_str(), (i < m_trans - 1 ? '\t' : '\n'));
285 for (int i = 0; i < m_trans; i++) {
286 const Transcript& transcript = transcripts.getTranscriptAt(ta.spAt(i));
287 fprintf(fo, "%s%c", transcript.getGeneID().c_str(), (i < m_trans - 1 ? '\t' : '\n'));
289 for (int i = 0; i < m_trans; i++)
290 fprintf(fo, "%.2f%c", trans_lens[i], (i < m_trans - 1 ? '\t' : '\n'));
291 for (int i = 0; i < m_trans; i++)
292 fprintf(fo, "%.2f%c", trans_eels[i], (i < m_trans - 1 ? '\t' : '\n'));
293 for (int i = 0; i < m_trans; i++)
294 fprintf(fo, "%.2f%c", trans_counts[i], (i < m_trans - 1 ? '\t' : '\n'));
295 for (int i = 0; i < m_trans; i++)
296 fprintf(fo, "%.2f%c", trans_tpm[i], (i < m_trans - 1 ? '\t' : '\n'));
297 for (int i = 0; i < m_trans; i++)
298 fprintf(fo, "%.2f%c", trans_fpkm[i], (i < m_trans - 1 ? '\t' : '\n'));
299 for (int i = 0; i < m_trans; i++)
300 fprintf(fo, "%.2f%c", gt_pct[i] * 1e2, (i < m_trans - 1 ? '\t' : '\n'));
305 sprintf(outF, "%s.gene_res", imdName);
306 fo = fopen(outF, "w");
307 for (int i = 0; i < m; i++) {
308 const Transcript& transcript = transcripts.getTranscriptAt(gi.spAt(i));
309 fprintf(fo, "%s%c", transcript.getGeneID().c_str(), (i < m - 1 ? '\t' : '\n'));
311 for (int i = 0; i < m; i++) {
312 int b = gi.spAt(i), e = gi.spAt(i + 1);
313 std::string curtid = "", tid;
314 for (int j = b; j < e; j++) {
315 tid = transcripts.getTranscriptAt(j).getTranscriptID();
317 if (curtid != "") fprintf(fo, ",");
318 fprintf(fo, "%s", tid.c_str());
322 fprintf(fo, "%c", (i < m - 1 ? '\t' : '\n'));
324 for (int i = 0; i < m; i++)
325 fprintf(fo, "%.2f%c", glens[i], (i < m - 1 ? '\t' : '\n'));
326 for (int i = 0; i < m; i++)
327 fprintf(fo, "%.2f%c", gene_eels[i], (i < m - 1 ? '\t' : '\n'));
328 for (int i = 0; i < m; i++)
329 fprintf(fo, "%.2f%c", gene_counts[i], (i < m - 1 ? '\t' : '\n'));
330 for (int i = 0; i < m; i++)
331 fprintf(fo, "%.2f%c", gene_tpm[i], (i < m - 1 ? '\t' : '\n'));
332 for (int i = 0; i < m; i++)
333 fprintf(fo, "%.2f%c", gene_fpkm[i], (i < m - 1 ? '\t' : '\n'));
336 if (verbose) { printf("Expression Results are written!\n"); }
339 void writeResultsGibbs(int M, char* refName, char* imdName, std::vector<double>& pme_c, std::vector<double>& pme_fpkm, std::vector<double>& pme_tpm) {
346 std::vector<double> isopct;
347 std::vector<double> gene_counts, gene_tpm, gene_fpkm;
350 sprintf(groupF, "%s.grp", refName);
354 // For allele-specific expression
357 std::vector<double> trans_counts, trans_tpm, trans_fpkm, ta_pct, gt_pct;
359 bool alleleS = isAlleleSpecific(refName, >, &ta); // if allele-specific
361 //calculate IsoPct, etc.
362 isopct.assign(M + 1, 0.0);
363 gene_counts.assign(m, 0.0); gene_tpm.assign(m, 0.0); gene_fpkm.assign(m, 0.0);
365 for (int i = 0; i < m; i++) {
366 int b = gi.spAt(i), e = gi.spAt(i + 1);
367 for (int j = b; j < e; j++) {
368 gene_counts[i] += pme_c[j];
369 gene_tpm[i] += pme_tpm[j];
370 gene_fpkm[i] += pme_fpkm[j];
372 if (gene_tpm[i] < EPSILON) continue;
373 for (int j = b; j < e; j++)
374 isopct[j] = pme_tpm[j] / gene_tpm[i];
379 ta_pct.assign(M + 1, 0.0);
380 trans_counts.assign(m_trans, 0.0); trans_tpm.assign(m_trans, 0.0); trans_fpkm.assign(m_trans, 0.0);
382 for (int i = 0; i < m_trans; i++) {
383 int b = ta.spAt(i), e = ta.spAt(i + 1);
384 for (int j = b; j < e; j++) {
385 trans_counts[i] += pme_c[j];
386 trans_tpm[i] += pme_tpm[j];
387 trans_fpkm[i] += pme_fpkm[j];
389 if (trans_tpm[i] < EPSILON) continue;
390 for (int j = b; j < e; j++)
391 ta_pct[j] = pme_tpm[j] / trans_tpm[i];
394 gt_pct.assign(m_trans, 0.0);
395 for (int i = 0; i < m; i++)
396 if (gene_tpm[i] >= EPSILON) {
397 int b = gt.spAt(i), e = gt.spAt(i + 1);
398 for (int j = b; j < e; j++) gt_pct[j] = trans_tpm[j] / gene_tpm[i];
403 //isoform level results
404 sprintf(outF, "%s.iso_res", imdName);
405 fo = fopen(outF, "a");
406 general_assert(fo != NULL, "Cannot open " + cstrtos(outF) + "!");
408 for (int i = 1; i <= M; i++)
409 fprintf(fo, "%.2f%c", pme_c[i], (i < M ? '\t' : '\n'));
410 for (int i = 1; i <= M; i++)
411 fprintf(fo, "%.2f%c", pme_tpm[i], (i < M ? '\t' : '\n'));
412 for (int i = 1; i <= M; i++)
413 fprintf(fo, "%.2f%c", pme_fpkm[i], (i < M ? '\t' : '\n'));
414 for (int i = 1; i <= M; i++)
415 fprintf(fo, "%.2f%c", isopct[i] * 1e2, (i < M ? '\t' : '\n'));
419 //allele level results
420 sprintf(outF, "%s.allele_res", imdName);
421 fo = fopen(outF, "a");
422 general_assert(fo != NULL, "Cannot open " + cstrtos(outF) + "!");
424 for (int i = 1; i <= M; i++)
425 fprintf(fo, "%.2f%c", pme_c[i], (i < M ? '\t' : '\n'));
426 for (int i = 1; i <= M; i++)
427 fprintf(fo, "%.2f%c", pme_tpm[i], (i < M ? '\t' : '\n'));
428 for (int i = 1; i <= M; i++)
429 fprintf(fo, "%.2f%c", pme_fpkm[i], (i < M ? '\t' : '\n'));
430 for (int i = 1; i <= M; i++)
431 fprintf(fo, "%.2f%c", ta_pct[i] * 1e2, (i < M ? '\t' : '\n'));
432 for (int i = 1; i <= M; i++)
433 fprintf(fo, "%.2f%c", isopct[i] * 1e2, (i < M ? '\t' : '\n'));
436 //isoform level results
437 sprintf(outF, "%s.iso_res", imdName);
438 fo = fopen(outF, "a");
439 general_assert(fo != NULL, "Cannot open " + cstrtos(outF) + "!");
441 for (int i = 0; i < m_trans; i++)
442 fprintf(fo, "%.2f%c", trans_counts[i], (i < m_trans - 1 ? '\t' : '\n'));
443 for (int i = 0; i < m_trans; i++)
444 fprintf(fo, "%.2f%c", trans_tpm[i], (i < m_trans - 1 ? '\t' : '\n'));
445 for (int i = 0; i < m_trans; i++)
446 fprintf(fo, "%.2f%c", trans_fpkm[i], (i < m_trans - 1 ? '\t' : '\n'));
447 for (int i = 0; i < m_trans; i++)
448 fprintf(fo, "%.2f%c", gt_pct[i] * 1e2, (i < m_trans - 1 ? '\t' : '\n'));
453 sprintf(outF, "%s.gene_res", imdName);
454 fo = fopen(outF, "a");
455 general_assert(fo != NULL, "Cannot open " + cstrtos(outF) + "!");
457 for (int i = 0; i < m; i++)
458 fprintf(fo, "%.2f%c", gene_counts[i], (i < m - 1 ? '\t' : '\n'));
459 for (int i = 0; i < m; i++)
460 fprintf(fo, "%.2f%c", gene_tpm[i], (i < m - 1 ? '\t' : '\n'));
461 for (int i = 0; i < m; i++)
462 fprintf(fo, "%.2f%c", gene_fpkm[i], (i < m - 1 ? '\t' : '\n'));
465 if (verbose) { printf("Gibbs based expression values are written!\n"); }
468 void writeResultsSimulation(int M, char* refName, char* outFN, Transcripts& transcripts, std::vector<double>& eel, std::vector<double>& counts) {
477 sprintf(groupF, "%s.grp", refName);
481 std::vector<int> tlens;
482 std::vector<double> tpm, fpkm, isopct;
483 std::vector<double> glens, gene_eels, gene_counts, gene_tpm, gene_fpkm;
485 // For allele-specific expression
488 std::vector<double> trans_lens, trans_eels, trans_counts, trans_tpm, trans_fpkm, ta_pct, gt_pct;
490 bool alleleS = isAlleleSpecific(refName, >, &ta); // if allele-specific
492 for (int i = 1; i <= M; i++)
493 general_assert(eel[i] > EPSILON || counts[i] <= EPSILON, "An isoform whose effecitve length < " + ftos(MINEEL, 6) + " got sampled!");
495 calcExpressionValues(M, counts, eel, tpm, fpkm);
497 //calculate IsoPct, etc.
498 isopct.assign(M + 1, 0.0);
499 tlens.assign(M + 1, 0);
501 glens.assign(m, 0.0); gene_eels.assign(m, 0.0);
502 gene_counts.assign(m, 0.0); gene_tpm.assign(m, 0.0); gene_fpkm.assign(m, 0.0);
504 for (int i = 0; i < m; i++) {
505 int b = gi.spAt(i), e = gi.spAt(i + 1);
506 for (int j = b; j < e; j++) {
507 const Transcript& transcript = transcripts.getTranscriptAt(j);
508 tlens[j] = transcript.getLength();
510 gene_counts[i] += counts[j];
511 gene_tpm[i] += tpm[j];
512 gene_fpkm[i] += fpkm[j];
515 if (gene_tpm[i] < EPSILON) {
516 double frac = 1.0 / (e - b);
517 for (int j = b; j < e; j++) {
518 glens[i] += tlens[j] * frac;
519 gene_eels[i] += eel[j] * frac;
523 for (int j = b; j < e; j++) {
524 isopct[j] = tpm[j] / gene_tpm[i];
525 glens[i] += tlens[j] * isopct[j];
526 gene_eels[i] += eel[j] * isopct[j];
533 ta_pct.assign(M + 1, 0.0);
534 trans_lens.assign(m_trans, 0.0); trans_eels.assign(m_trans, 0.0);
535 trans_counts.assign(m_trans, 0.0); trans_tpm.assign(m_trans, 0.0); trans_fpkm.assign(m_trans, 0.0);
537 for (int i = 0; i < m_trans; i++) {
538 int b = ta.spAt(i), e = ta.spAt(i + 1);
539 for (int j = b; j < e; j++) {
540 trans_counts[i] += counts[j];
541 trans_tpm[i] += tpm[j];
542 trans_fpkm[i] += fpkm[j];
545 if (trans_tpm[i] < EPSILON) {
546 double frac = 1.0 / (e - b);
547 for (int j = b; j < e; j++) {
548 trans_lens[i] += tlens[j] * frac;
549 trans_eels[i] += eel[j] * frac;
553 for (int j = b; j < e; j++) {
554 ta_pct[j] = tpm[j] / trans_tpm[i];
555 trans_lens[i] += tlens[j] * ta_pct[j];
556 trans_eels[i] += eel[j] * ta_pct[j];
561 gt_pct.assign(m_trans, 0.0);
562 for (int i = 0; i < m; i++)
563 if (gene_tpm[i] >= EPSILON) {
564 int b = gt.spAt(i), e = gt.spAt(i + 1);
565 for (int j = b; j < e; j++) gt_pct[j] = trans_tpm[j] / gene_tpm[i];
571 sprintf(outF, "%s.sim.alleles.results", outFN);
572 fo = fopen(outF, "w");
573 fprintf(fo, "allele_id\ttranscript_id\tgene_id\tlength\teffective_length\tcount\tTPM\tFPKM\tAlleleIsoPct\tAlleleGenePct\n");
574 for (int i = 1; i <= M; i++) {
575 const Transcript& transcript = transcripts.getTranscriptAt(i);
576 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],
577 eel[i], counts[i], tpm[i], fpkm[i], ta_pct[i] * 1e2, isopct[i] * 1e2);
583 sprintf(outF, "%s.sim.isoforms.results", outFN);
584 fo = fopen(outF, "w");
585 fprintf(fo, "transcript_id\tgene_id\tlength\teffective_length\tcount\tTPM\tFPKM\tIsoPct\n");
587 for (int i = 1; i <= M; i++) {
588 const Transcript& transcript = transcripts.getTranscriptAt(i);
589 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],
590 eel[i], counts[i], tpm[i], fpkm[i], isopct[i] * 1e2);
594 for (int i = 0; i < m_trans; i++) {
595 const Transcript& transcript = transcripts.getTranscriptAt(ta.spAt(i));
596 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],
597 trans_eels[i], trans_counts[i], trans_tpm[i], trans_fpkm[i], gt_pct[i] * 1e2);
603 sprintf(outF, "%s.sim.genes.results", outFN);
604 fo = fopen(outF, "w");
605 fprintf(fo, "gene_id\ttranscript_id(s)\tlength\teffective_length\tcount\tTPM\tFPKM\n");
606 for (int i = 0; i < m; i++) {
607 int b = gi.spAt(i), e = gi.spAt(i + 1);
608 const std::string& gene_id = transcripts.getTranscriptAt(b).getGeneID();
609 fprintf(fo, "%s\t", gene_id.c_str());
610 std::string curtid = "", tid;
611 for (int j = b; j < e; j++) {
612 tid = transcripts.getTranscriptAt(j).getTranscriptID();
614 if (curtid != "") fprintf(fo, ",");
615 fprintf(fo, "%s", tid.c_str());
619 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]);