WriteResults.h

   1 #ifndef WRITERESULTS_H_
   2 #define WRITERESULTS_H_
   3
   4 #include<cstdio>
   5 #include<vector>
   6 #include<string>
   7 #include<fstream>
   8 #include<algorithm>
   9
  10 #include "utils.h"
  11 #include "my_assert.h"
  12 #include "GroupInfo.h"
  13 #include "Transcript.h"
  14 #include "Transcripts.h"
  15 #include "Refs.h"
  16
  17 #include "Model.h"
  18 #include "SingleModel.h"
  19 #include "SingleQModel.h"
  20 #include "PairedEndModel.h"
  21 #include "PairedEndQModel.h"
  22
  23 template<class ModelType>
  24 void calcExpectedEffectiveLengths(int M, Refs& refs, ModelType& model, std::vector<double>& eel) {
  25         int lb, ub, span;
  26         double *pdf = NULL, *cdf = NULL, *clen = NULL; // clen[i] = sigma_{j=1}^{i}pdf[i]*(lb+i)
  27
  28         model.getGLD().copyTo(pdf, cdf, lb, ub, span);
  29         clen = new double[span + 1];
  30         clen[0] = 0.0;
  31         for (int i = 1; i <= span; i++) {
  32                 clen[i] = clen[i - 1] + pdf[i] * (lb + i);
  33         }
  34
  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);
  41
  42                 if (pos2 == 0) { eel[i] = 0.0; continue; }
  43
  44                 eel[i] = fullLen * cdf[pos1] + ((cdf[pos2] - cdf[pos1]) * (totLen + 1) - (clen[pos2] - clen[pos1]));
  45                 assert(eel[i] >= 0);
  46                 if (eel[i] < MINEEL) { eel[i] = 0.0; }
  47         }
  48
  49         delete[] pdf;
  50         delete[] cdf;
  51         delete[] clen;
  52 }
  53
  54 void polishTheta(int M, std::vector<double>& theta, const std::vector<double>& eel, const double* mw) {
  55         double sum = 0.0;
  56
  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
  60          */
  61
  62         for (int i = 0; i <= M; i++) {
  63                 // i == 0, mw[i] == 1
  64                 if (i > 0 && (mw[i] < EPSILON || eel[i] < EPSILON)) {
  65                         theta[i] = 0.0;
  66                         continue;
  67                 }
  68                 theta[i] = theta[i] / mw[i];
  69                 sum += theta[i];
  70         }
  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;
  74 }
  75
  76 void calcExpressionValues(int M, const std::vector<double>& theta, const std::vector<double>& eel, std::vector<double>& tpm, std::vector<double>& fpkm) {
  77         double denom;
  78         std::vector<double> frac;
  79
  80         //calculate fraction of count over all mappabile reads
  81         denom = 0.0;
  82         frac.assign(M + 1, 0.0);
  83         for (int i = 1; i <= M; i++)
  84           if (eel[i] >= EPSILON) {
  85             frac[i] = theta[i];
  86             denom += frac[i];
  87           }
  88         general_assert(denom >= EPSILON, "No alignable reads?!");
  89         for (int i = 1; i <= M; i++) frac[i] /= denom;
  90
  91         //calculate FPKM
  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];
  95
  96         //calculate TPM
  97         tpm.assign(M + 1, 0.0);
  98         denom = 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;
 101 }
 102
 103 inline bool isAlleleSpecific(const char* refName, GroupInfo* gt = NULL, GroupInfo* ta = NULL) {
 104   bool alleleS;
 105   char gtF[STRLEN], taF[STRLEN];
 106
 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();
 113
 114   if (alleleS) {
 115     if (gt != NULL) gt->load(gtF);
 116     if (ta != NULL) ta->load(taF);
 117   }
 118
 119   return alleleS;
 120 }
 121
 122 void writeResultsEM(int M, const char* refName, const char* imdName, Transcripts& transcripts, std::vector<double>& theta, std::vector<double>& eel, double* counts) {
 123         char outF[STRLEN];
 124         FILE *fo;
 125
 126         int m;
 127         GroupInfo gi;
 128         char groupF[STRLEN];
 129
 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;
 133
 134         // Load group info
 135         sprintf(groupF, "%s.grp", refName);
 136         gi.load(groupF);
 137         m = gi.getm();
 138
 139         // For allele-specific expression
 140         int m_trans = 0;
 141         GroupInfo gt, ta;
 142         std::vector<double> trans_lens, trans_eels, trans_counts, trans_tpm, trans_fpkm, ta_pct, gt_pct;
 143
 144         bool alleleS = isAlleleSpecific(refName, &gt, &ta); // if allele-specific
 145
 146         calcExpressionValues(M, theta, eel, tpm, fpkm);
 147
 148         //calculate IsoPct, etc.
 149         isopct.assign(M + 1, 0.0);
 150         tlens.assign(M + 1, 0);
 151
 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);
 154
 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();
 160
 161                         gene_counts[i] += counts[j];
 162                         gene_tpm[i] += tpm[j];
 163                         gene_fpkm[i] += fpkm[j];
 164                 }
 165
 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;
 171                         }
 172                 }
 173                 else {
 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];
 178                         }
 179                 }
 180         }
 181
 182         if (alleleS) {
 183           m_trans = ta.getm();
 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);
 187
 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];
 194                 }
 195
 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;
 201                         }
 202                 }
 203                 else {
 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];
 208                         }
 209                 }
 210           }
 211
 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];
 217             }
 218         }
 219
 220         if (!alleleS) {
 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'));
 227           }
 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'));
 231           }
 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'));
 244           fclose(fo);
 245         }
 246         else {
 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'));
 253           }
 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'));
 257           }
 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'));
 261           }
 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'));
 276           fclose(fo);
 277
 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'));
 284           }
 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'));
 288           }
 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'));
 301           fclose(fo);
 302         }
 303
 304         //gene level results
 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'));
 310         }
 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();
 316                         if (curtid != tid) {
 317                           if (curtid != "") fprintf(fo, ",");
 318                           fprintf(fo, "%s", tid.c_str());
 319                           curtid = tid;
 320                         }
 321                 }
 322                 fprintf(fo, "%c", (i < m - 1 ? '\t' : '\n'));
 323         }
 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'));
 334         fclose(fo);
 335
 336         if (verbose) { printf("Expression Results are written!\n"); }
 337 }
 338
 339 void writeResultsGibbs(int M, char* refName, char* imdName, std::vector<double>& pme_c, std::vector<double>& pme_fpkm, std::vector<double>& pme_tpm) {
 340         char outF[STRLEN];
 341         FILE *fo;
 342
 343         int m;
 344         GroupInfo gi;
 345         char groupF[STRLEN];
 346         std::vector<double> isopct;
 347         std::vector<double> gene_counts, gene_tpm, gene_fpkm;
 348
 349         // Load group info
 350         sprintf(groupF, "%s.grp", refName);
 351         gi.load(groupF);
 352         m = gi.getm();
 353
 354         // For allele-specific expression
 355         int m_trans = 0;
 356         GroupInfo gt, ta;
 357         std::vector<double> trans_counts, trans_tpm, trans_fpkm, ta_pct, gt_pct;
 358
 359         bool alleleS = isAlleleSpecific(refName, &gt, &ta); // if allele-specific
 360
 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);
 364
 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];
 371                 }
 372                 if (gene_tpm[i] < EPSILON) continue;
 373                 for (int j = b; j < e; j++)
 374                         isopct[j] = pme_tpm[j] / gene_tpm[i];
 375         }
 376
 377         if (alleleS) {
 378           m_trans = ta.getm();
 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);
 381
 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];
 388                 }
 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];
 392           }
 393
 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];
 399             }
 400         }
 401
 402         if (!alleleS) {
 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) + "!");
 407
 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'));
 416           fclose(fo);
 417         }
 418         else {
 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) + "!");
 423
 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'));
 434           fclose(fo);
 435
 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) + "!");
 440
 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'));
 449           fclose(fo);
 450         }
 451
 452         //gene level results
 453         sprintf(outF, "%s.gene_res", imdName);
 454         fo = fopen(outF, "a");
 455         general_assert(fo != NULL, "Cannot open " + cstrtos(outF) + "!");
 456
 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'));
 463         fclose(fo);
 464
 465         if (verbose) { printf("Gibbs based expression values are written!\n"); }
 466 }
 467
 468 void writeResultsSimulation(int M, char* refName, char* outFN, Transcripts& transcripts, std::vector<double>& eel, std::vector<double>& counts) {
 469         char outF[STRLEN];
 470         FILE *fo;
 471
 472         int m;
 473         GroupInfo gi;
 474         char groupF[STRLEN];
 475
 476         // Load group info
 477         sprintf(groupF, "%s.grp", refName);
 478         gi.load(groupF);
 479         m = gi.getm();
 480
 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;
 484
 485         // For allele-specific expression
 486         int m_trans = 0;
 487         GroupInfo gt, ta;
 488         std::vector<double> trans_lens, trans_eels, trans_counts, trans_tpm, trans_fpkm, ta_pct, gt_pct;
 489
 490         bool alleleS = isAlleleSpecific(refName, &gt, &ta); // if allele-specific
 491
 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!");
 494
 495         calcExpressionValues(M, counts, eel, tpm, fpkm);
 496
 497         //calculate IsoPct, etc.
 498         isopct.assign(M + 1, 0.0);
 499         tlens.assign(M + 1, 0);
 500
 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);
 503
 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();
 509
 510             gene_counts[i] += counts[j];
 511             gene_tpm[i] += tpm[j];
 512             gene_fpkm[i] += fpkm[j];
 513           }
 514
 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;
 520             }
 521           }
 522           else {
 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];
 527             }
 528           }
 529         }
 530
 531         if (alleleS) {
 532           m_trans = ta.getm();
 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);
 536
 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];
 543             }
 544
 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;
 550               }
 551             }
 552             else {
 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];
 557               }
 558             }
 559           }
 560
 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];
 566             }
 567         }
 568
 569         //allele level
 570         if (alleleS) {
 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);
 578           }
 579           fclose(fo);
 580         }
 581
 582         //isoform level
 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");
 586         if (!alleleS) {
 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);
 591           }
 592         }
 593         else {
 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);
 598           }
 599         }
 600         fclose(fo);
 601
 602         //gene level
 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();
 613                         if (curtid != tid) {
 614                           if (curtid != "") fprintf(fo, ",");
 615                           fprintf(fo, "%s", tid.c_str());
 616                           curtid = tid;
 617                         }
 618                 }
 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]);
 620         }
 621         fclose(fo);
 622 }
 623
 624 #endif