X-Git-Url: https://git.donarmstrong.com/?a=blobdiff_plain;f=EM.cpp;h=59783f3089b5a96c7ea67db0113f9b558818d861;hb=58d504aaf36ae486b1dba6d03e0e9f1c25855037;hp=09b85be234bb432724524574519ef6affda7f2d7;hpb=9c2e46183a19d661f0a618a8eabe8ce1f6a8e2d6;p=rsem.git

diff --git a/EM.cpp b/EM.cpp
index 09b85be..59783f3 100644
--- 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