Fixed a minor bug which only affects paired-end reads for reporting how many alignmen...

[rsem.git] / EM.cpp

diff --git a/EM.cpp b/EM.cpp
index 09b85be234bb432724524574519ef6affda7f2d7..59783f3089b5a96c7ea67db0113f9b558818d861 100644 (file)
--- a/EM.cpp
+++ b/EM.cpp
@@ -336,7 +336,7 @@ void calcExpressionValues(const vector<double>& theta, const vector<double>& eel
            frac[i] = theta[i];
            denom += frac[i];
          }
-       general_assert(denom > 0, "No alignable reads?!");
+       general_assert(denom >= EPSILON, "No alignable reads?!");
        for (int i = 1; i <= M; i++) frac[i] /= denom;
   
        //calculate FPKM
@@ -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