-void calcExpressionValues(const vector<double>& theta, const vector<double>& eel, vector<double>& tpm, vector<double>& fpkm) {
- double denom;
- vector<double> frac;
-
- //calculate fraction of count over all mappabile reads
- denom = 0.0;
- frac.assign(M + 1, 0.0);
- for (int i = 1; i <= M; i++)
- if (eel[i] >= EPSILON) {
- frac[i] = theta[i];
- denom += frac[i];
- }
- general_assert(denom > 0, "No alignable reads?!");
- for (int i = 1; i <= M; i++) frac[i] /= denom;
-
- //calculate FPKM
- fpkm.assign(M + 1, 0.0);
- for (int i = 1; i <= M; i++)
- if (eel[i] >= EPSILON) fpkm[i] = frac[i] * 1e9 / eel[i];
-
- //calculate TPM
- tpm.assign(M + 1, 0.0);
- denom = 0.0;
- for (int i = 1; i <= M; i++) denom += fpkm[i];
- for (int i = 1; i <= M; i++) tpm[i] = fpkm[i] / denom * 1e6;
-}
-
-void writeResFiles(char* outFN) {
- FILE *fo;
- vector<int> tlens;
- vector<double> fpkm, tpm, isopct;
- vector<double> glens, gene_eels, gene_counts, gene_tpm, gene_fpkm;
-
- for (int i = 1; i <= M; i++)
- general_assert(eel[i] > EPSILON || counts[i] <= EPSILON, "An isoform whose effecitve length < " + ftos(MINEEL, 6) + " got sampled!");
-
- calcExpressionValues(counts, eel, tpm, fpkm);
-
- //calculate IsoPct, etc.
- isopct.assign(M + 1, 0.0);
- tlens.assign(M + 1, 0);
-
- glens.assign(m, 0.0); gene_eels.assign(m, 0.0);
- gene_counts.assign(m, 0.0); gene_tpm.assign(m, 0.0); gene_fpkm.assign(m, 0.0);
-
- for (int i = 0; i < m; i++) {
- int b = gi.spAt(i), e = gi.spAt(i + 1);
- for (int j = b; j < e; j++) {
- const Transcript& transcript = transcripts.getTranscriptAt(j);
- tlens[j] = transcript.getLength();
-
- glens[i] += tlens[j] * tpm[j];
- gene_eels[i] += eel[j] * tpm[j];
- gene_counts[i] += counts[j];
- gene_tpm[i] += tpm[j];
- gene_fpkm[i] += fpkm[j];
- }
-
- if (gene_tpm[i] < EPSILON) continue;
-
- for (int j = b; j < e; j++)
- isopct[j] = tpm[j] / gene_tpm[i];
- glens[i] /= gene_tpm[i];
- gene_eels[i] /= gene_tpm[i];
- }
-
- //isoform level
- sprintf(isoResF, "%s.sim.isoforms.results", outFN);
- fo = fopen(isoResF, "w");
- fprintf(fo, "transcript_id\tgene_id\tlength\teffective_length\tcount\tTPM\tFPKM\tIsoPct\n");
- for (int i = 1; i <= M; i++) {
- const Transcript& transcript = transcripts.getTranscriptAt(i);
- 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],
- eel[i], counts[i], tpm[i], fpkm[i], isopct[i] * 1e2);
- }
- fclose(fo);
-
- //gene level
- sprintf(geneResF, "%s.sim.genes.results", outFN);
- fo = fopen(geneResF, "w");
- fprintf(fo, "gene_id\ttranscript_id(s)\tlength\teffective_length\tcount\tTPM\tFPKM\n");
- for (int i = 0; i < m; i++) {
- int b = gi.spAt(i), e = gi.spAt(i + 1);
- const string& gene_id = transcripts.getTranscriptAt(b).getGeneID();
- fprintf(fo, "%s\t", gene_id.c_str());
- for (int j = b; j < e; j++) {
- fprintf(fo, "%s%c", transcripts.getTranscriptAt(j).getTranscriptID().c_str(), (j < e - 1 ? ',' : '\t'));
- }
- fprintf(fo, "%.2f\t%.2f\t%.2f\t%.2f\t%.2f\n", glens[i], gene_eels[i], gene_counts[i], gene_tpm[i], gene_fpkm[i]);
- }
- fclose(fo);
-}
-