From d484af5f8ab839e296f01d01ff713c37a8db6e02 Mon Sep 17 00:00:00 2001 From: Bo Li Date: Sat, 15 Sep 2012 15:27:44 -0500 Subject: [PATCH] If a gene is not expressed, its length/effective length is defined as the unweighted mean of its isoforms' --- EM.cpp | 22 ++++++++++++++-------- 1 file changed, 14 insertions(+), 8 deletions(-) diff --git a/EM.cpp b/EM.cpp index 09b85be..fdc75cc 100644 --- a/EM.cpp +++ b/EM.cpp @@ -378,19 +378,25 @@ void writeResults(ModelType& model, double* counts) { 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]; + if (gene_tpm[i] < EPSILON) { + double frac = 1.0 / (e - b); + for (int j = b; j < e; j++) { + glens[i] += tlens[j] * frac; + gene_eels[i] += eel[j] * frac; + } + } + else { + for (int j = b; j < e; j++) { + isopct[j] = tpm[j] / gene_tpm[i]; + glens[i] += tlens[j] * isopct[j]; + gene_eels[i] += eel[j] * isopct[j]; + } + } } //isoform level results -- 2.39.2