WriteResults.h

   1 #ifndef WRITERESULTS_H_
   2 #define WRITERESULTS_H_
   3
   4 #include<cmath>
   5 #include<cstdio>
   6 #include<vector>
   7 #include<string>
   8 #include<fstream>
   9 #include<algorithm>
  10
  11 #include "utils.h"
  12 #include "my_assert.h"
  13 #include "GroupInfo.h"
  14 #include "Transcript.h"
  15 #include "Transcripts.h"
  16 #include "Refs.h"
  17
  18 #include "Model.h"
  19 #include "SingleModel.h"
  20 #include "SingleQModel.h"
  21 #include "PairedEndModel.h"
  22 #include "PairedEndQModel.h"
  23
  24 template<class ModelType>
  25 void calcExpectedEffectiveLengths(int M, Refs& refs, ModelType& model, std::vector<double>& eel) {
  26         int lb, ub, span;
  27         double *pdf = NULL, *cdf = NULL, *clen = NULL; // clen[i] = sigma_{j=1}^{i}pdf[i]*(lb+i)
  28
  29         model.getGLD().copyTo(pdf, cdf, lb, ub, span);
  30         clen = new double[span + 1];
  31         clen[0] = 0.0;
  32         for (int i = 1; i <= span; i++) {
  33                 clen[i] = clen[i - 1] + pdf[i] * (lb + i);
  34         }
  35
  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);
  42
  43                 if (pos2 == 0) { eel[i] = 0.0; continue; }
  44
  45                 eel[i] = fullLen * cdf[pos1] + ((cdf[pos2] - cdf[pos1]) * (totLen + 1) - (clen[pos2] - clen[pos1]));
  46                 assert(eel[i] >= 0);
  47                 if (eel[i] < MINEEL) { eel[i] = 0.0; }
  48         }
  49
  50         delete[] pdf;
  51         delete[] cdf;
  52         delete[] clen;
  53 }
  54
  55 void polishTheta(int M, std::vector<double>& theta, const std::vector<double>& eel, const double* mw) {
  56         double sum = 0.0;
  57
  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
  61          */
  62
  63         for (int i = 0; i <= M; i++) {
  64                 // i == 0, mw[i] == 1
  65                 if (i > 0 && (mw[i] < EPSILON || eel[i] < EPSILON)) {
  66                         theta[i] = 0.0;
  67                         continue;
  68                 }
  69                 theta[i] = theta[i] / mw[i];
  70                 sum += theta[i];
  71         }
  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;
  75 }
  76
  77 void calcExpressionValues(int M, const std::vector<double>& theta, const std::vector<double>& eel, std::vector<double>& tpm, std::vector<double>& fpkm) {
  78         double denom;
  79         std::vector<double> frac;
  80
  81         //calculate fraction of count over all mappabile reads
  82         denom = 0.0;
  83         frac.assign(M + 1, 0.0);
  84         for (int i = 1; i <= M; i++)
  85           if (eel[i] >= EPSILON) {
  86             frac[i] = theta[i];
  87             denom += frac[i];
  88           }
  89         general_assert(denom >= EPSILON, "No alignable reads?!");
  90         for (int i = 1; i <= M; i++) frac[i] /= denom;
  91
  92         //calculate FPKM
  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];
  96
  97         //calculate TPM
  98         tpm.assign(M + 1, 0.0);
  99         denom = 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;
 102 }
 103
 104 inline bool isAlleleSpecific(const char* refName, GroupInfo* gt = NULL, GroupInfo* ta = NULL) {
 105   bool alleleS;
 106   char gtF[STRLEN], taF[STRLEN];
 107
 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();
 114
 115   if (alleleS) {
 116     if (gt != NULL) gt->load(gtF);
 117     if (ta != NULL) ta->load(taF);
 118   }
 119
 120   return alleleS;
 121 }
 122
 123 void writeResultsEM(int M, const char* refName, const char* imdName, Transcripts& transcripts, std::vector<double>& theta, std::vector<double>& eel, double* counts) {
 124         char outF[STRLEN];
 125         FILE *fo;
 126
 127         int m;
 128         GroupInfo gi;
 129         char groupF[STRLEN];
 130
 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;
 134
 135         // Load group info
 136         sprintf(groupF, "%s.grp", refName);
 137         gi.load(groupF);
 138         m = gi.getm();
 139
 140         // For allele-specific expression
 141         int m_trans = 0;
 142         GroupInfo gt, ta;
 143         std::vector<double> trans_lens, trans_eels, trans_counts, trans_tpm, trans_fpkm, ta_pct, gt_pct;
 144
 145         bool alleleS = isAlleleSpecific(refName, &gt, &ta); // if allele-specific
 146
 147         calcExpressionValues(M, theta, eel, tpm, fpkm);
 148
 149         //calculate IsoPct, etc.
 150         isopct.assign(M + 1, 0.0);
 151         tlens.assign(M + 1, 0);
 152
 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);
 155
 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();
 161
 162                         gene_counts[i] += counts[j];
 163                         gene_tpm[i] += tpm[j];
 164                         gene_fpkm[i] += fpkm[j];
 165                 }
 166
 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;
 172                         }
 173                 }
 174                 else {
 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];
 179                         }
 180                 }
 181         }
 182
 183         if (alleleS) {
 184           m_trans = ta.getm();
 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);
 188
 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];
 195                 }
 196
 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;
 202                         }
 203                 }
 204                 else {
 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];
 209                         }
 210                 }
 211           }
 212
 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];
 218             }
 219         }
 220
 221         if (!alleleS) {
 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'));
 228           }
 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'));
 232           }
 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'));
 245           fclose(fo);
 246         }
 247         else {
 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'));
 254           }
 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'));
 258           }
 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'));
 262           }
 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'));
 277           fclose(fo);
 278
 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'));
 285           }
 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'));
 289           }
 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'));
 302           fclose(fo);
 303         }
 304
 305         //gene level results
 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'));
 311         }
 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();
 317                         if (curtid != tid) {
 318                           if (curtid != "") fprintf(fo, ",");
 319                           fprintf(fo, "%s", tid.c_str());
 320                           curtid = tid;
 321                         }
 322                 }
 323                 fprintf(fo, "%c", (i < m - 1 ? '\t' : '\n'));
 324         }
 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'));
 335         fclose(fo);
 336
 337         if (verbose) { printf("Expression Results are written!\n"); }
 338 }
 339
 340 void writeResultsGibbs(int M, int m, int m_trans, GroupInfo& gi, GroupInfo &gt, 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) {
 341         char outF[STRLEN];
 342         FILE *fo;
 343
 344         std::vector<double> isopct;
 345         std::vector<double> gene_counts, gene_tpm, gene_fpkm;
 346
 347         // For allele-specific expression
 348         std::vector<double> trans_counts, trans_tpm, trans_fpkm, ta_pct, gt_pct;
 349
 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);
 353
 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];
 360                 }
 361                 if (gene_tpm[i] < EPSILON) continue;
 362                 for (int j = b; j < e; j++)
 363                         isopct[j] = pme_tpm[j] / gene_tpm[i];
 364         }
 365
 366         if (alleleS) {
 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);
 369
 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];
 376                 }
 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];
 380           }
 381
 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];
 387             }
 388         }
 389
 390         if (!alleleS) {
 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) + "!");
 395
 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'));
 406           fclose(fo);
 407         }
 408         else {
 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) + "!");
 413
 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'));
 426           fclose(fo);
 427
 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) + "!");
 432
 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'));
 443           fclose(fo);
 444         }
 445
 446         //gene level results
 447         sprintf(outF, "%s.gene_res", imdName);
 448         fo = fopen(outF, "a");
 449         general_assert(fo != NULL, "Cannot open " + cstrtos(outF) + "!");
 450
 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'));
 459         fclose(fo);
 460
 461         if (verbose) { printf("Gibbs based expression values are written!\n"); }
 462 }
 463
 464 void writeResultsSimulation(int M, char* refName, char* outFN, Transcripts& transcripts, std::vector<double>& eel, std::vector<double>& counts) {
 465         char outF[STRLEN];
 466         FILE *fo;
 467
 468         int m;
 469         GroupInfo gi;
 470         char groupF[STRLEN];
 471
 472         // Load group info
 473         sprintf(groupF, "%s.grp", refName);
 474         gi.load(groupF);
 475         m = gi.getm();
 476
 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;
 480
 481         // For allele-specific expression
 482         int m_trans = 0;
 483         GroupInfo gt, ta;
 484         std::vector<double> trans_lens, trans_eels, trans_counts, trans_tpm, trans_fpkm, ta_pct, gt_pct;
 485
 486         bool alleleS = isAlleleSpecific(refName, &gt, &ta); // if allele-specific
 487
 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!");
 490
 491         calcExpressionValues(M, counts, eel, tpm, fpkm);
 492
 493         //calculate IsoPct, etc.
 494         isopct.assign(M + 1, 0.0);
 495         tlens.assign(M + 1, 0);
 496
 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);
 499
 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();
 505
 506             gene_counts[i] += counts[j];
 507             gene_tpm[i] += tpm[j];
 508             gene_fpkm[i] += fpkm[j];
 509           }
 510
 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;
 516             }
 517           }
 518           else {
 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];
 523             }
 524           }
 525         }
 526
 527         if (alleleS) {
 528           m_trans = ta.getm();
 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);
 532
 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];
 539             }
 540
 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;
 546               }
 547             }
 548             else {
 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];
 553               }
 554             }
 555           }
 556
 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];
 562             }
 563         }
 564
 565         //allele level
 566         if (alleleS) {
 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);
 574           }
 575           fclose(fo);
 576         }
 577
 578         //isoform level
 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");
 582         if (!alleleS) {
 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);
 587           }
 588         }
 589         else {
 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);
 594           }
 595         }
 596         fclose(fo);
 597
 598         //gene level
 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();
 609                         if (curtid != tid) {
 610                           if (curtid != "") fprintf(fo, ",");
 611                           fprintf(fo, "%s", tid.c_str());
 612                           curtid = tid;
 613                         }
 614                 }
 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]);
 616         }
 617         fclose(fo);
 618 }
 619
 620 #endif