X-Git-Url: https://git.donarmstrong.com/?a=blobdiff_plain;f=SingleModel.h;h=59db6ec69d30ce228535312f9cd91da476be2cae;hb=cb94fd597b180aa7cb01ae84c9d1025201b98d8e;hp=756e3d5cdecf386339eb90c303e9aa6b908a26e5;hpb=68a2be089a876aba126e384837559aaab40431bf;p=rsem.git

diff --git a/SingleModel.h b/SingleModel.h
index 756e3d5..59db6ec 100644
--- a/SingleModel.h
+++ b/SingleModel.h
@@ -269,12 +269,17 @@ void SingleModel::estimateFromReads(const char* readFN) {
 	for (int i = 0; i < 3; i++)
 		if (N[i] > 0) {
 			genReadFileNames(readFN, i, read_type, s, readFs);
-			ReadReader<SingleRead> reader(s, readFs);
+			ReadReader<SingleRead> reader(s, readFs, refs->hasPolyA(), seedLen); // allow calculation of calc_lq() function
 
 			int cnt = 0;
 			while (reader.next(read)) {
-				mld != NULL ? mld->update(read.getReadLength(), 1.0) : gld->update(read.getReadLength(), 1.0);
-				if (i == 0) { npro->updateC(read.getReadSeq()); }
+				if (!read.isLowQuality()) {
+					mld != NULL ? mld->update(read.getReadLength(), 1.0) : gld->update(read.getReadLength(), 1.0);
+					if (i == 0) { npro->updateC(read.getReadSeq()); }
+				}
+				else if (verbose && read.getReadLength() < seedLen) {
+					printf("Warning: Read %s is ignored due to read length %d < seed length %d!\n", read.getName().c_str(), read.getReadLength(), seedLen);
+				}
 
 				++cnt;
 				if (verbose && cnt % 1000000 == 0) { printf("%d READS PROCESSED\n", cnt); }
@@ -321,12 +326,12 @@ void SingleModel::read(const char* inpF) {
 	FILE *fi = fopen(inpF, "r");
 	if (fi == NULL) { fprintf(stderr, "Cannot open %s! It may not exist.\n", inpF); exit(-1); }
 
-	fscanf(fi, "%d", &val);
+	assert(fscanf(fi, "%d", &val) == 1);
 	assert(val == model_type);
 
 	ori->read(fi);
 	gld->read(fi);
-	fscanf(fi, "%d", &val);
+	assert(fscanf(fi, "%d", &val) == 1);
 	if (val > 0) {
 		if (mld == NULL) mld = new LenDist();
 		mld->read(fi);
@@ -339,7 +344,7 @@ void SingleModel::read(const char* inpF) {
 		if (M == 0) M = val;
 		if (M == val) {
 			mw = new double[M + 1];
-			for (int i = 0; i <= M; i++) fscanf(fi, "%lf", &mw[i]);
+			for (int i = 0; i <= M; i++) assert(fscanf(fi, "%lf", &mw[i]) == 1);
 		}
 	}
 
@@ -440,68 +445,67 @@ void SingleModel::finishSimulation() {
 }
 
 void SingleModel::calcMW() {
-  double probF, probR;
-  
-  assert(seedLen >= OLEN && (mld == NULL ? gld->getMinL() : mld->getMinL()) >= seedLen);
-  
-  memset(mw, 0, sizeof(double) * (M + 1));
-  mw[0] = 1.0;
-  
-  
-  probF = ori->getProb(0);
-  probR = ori->getProb(1);
+	double probF, probR;
+
+	assert((mld == NULL ? gld->getMinL() : mld->getMinL()) >= seedLen);
   
-  for (int i = 1; i <= M; i++) { 
-    RefSeq& ref = refs->getRef(i);
-    int totLen = ref.getTotLen();
-    int fullLen = ref.getFullLen();
-    double value = 0.0;
-    int minL, maxL;
-    int effL, pfpos;
-    int end = std::min(fullLen, totLen - seedLen + 1);
-    double factor;
-
-    for (int seedPos = 0; seedPos < end; seedPos++) 
-      if (ref.getMask(seedPos)) {
-	//forward
-	minL = gld->getMinL();
-	maxL = std::min(gld->getMaxL(), totLen - seedPos);
-	pfpos = seedPos;
-	for (int fragLen = minL; fragLen <= maxL; fragLen++) {
-	  effL = std::min(fullLen, totLen - fragLen + 1); 
-	  factor = (mld == NULL ? 1.0 : mld->getAdjustedCumulativeProb(std::min(mld->getMaxL(), fragLen), fragLen)); 
-	  value += probF * gld->getAdjustedProb(fragLen, totLen) * rspd->getAdjustedProb(pfpos, effL, fullLen) * factor; 
-	}
-	//reverse
-	minL = gld->getMinL();
-	maxL = std::min(gld->getMaxL(), seedPos + seedLen);
-	for (int fragLen = minL; fragLen <= maxL; fragLen++) {
-	  pfpos = seedPos - (fragLen - seedLen);
-	  effL = std::min(fullLen, totLen - fragLen + 1);
-	  factor = (mld == NULL ? 1.0 : mld->getAdjustedCumulativeProb(std::min(mld->getMaxL(), fragLen), fragLen)); 
-	  value += probR * gld->getAdjustedProb(fragLen, totLen) * rspd->getAdjustedProb(pfpos, effL, fullLen) * factor;
-	}
-      }
-    
-    //for reverse strand masking
-    for (int seedPos = end; seedPos <= totLen - seedLen; seedPos++) {
-      minL = std::max(gld->getMinL(), seedPos + seedLen - fullLen + 1);
-      maxL = std::min(gld->getMaxL(), seedPos + seedLen);
-      for (int fragLen = minL; fragLen <= maxL; fragLen++) {
-	pfpos = seedPos - (fragLen - seedLen);
-	effL = std::min(fullLen, totLen - fragLen + 1);
-	factor = (mld == NULL ? 1.0 : mld->getAdjustedCumulativeProb(std::min(mld->getMaxL(), fragLen), fragLen)); 
-	value += probR * gld->getAdjustedProb(fragLen, totLen) * rspd->getAdjustedProb(pfpos, effL, fullLen) * factor;
-      }
-    }
+	memset(mw, 0, sizeof(double) * (M + 1));
+	mw[0] = 1.0;
+
+	probF = ori->getProb(0);
+	probR = ori->getProb(1);
+
+	for (int i = 1; i <= M; i++) {
+		RefSeq& ref = refs->getRef(i);
+		int totLen = ref.getTotLen();
+		int fullLen = ref.getFullLen();
+		double value = 0.0;
+		int minL, maxL;
+		int effL, pfpos;
+		int end = std::min(fullLen, totLen - seedLen + 1);
+		double factor;
+
+		for (int seedPos = 0; seedPos < end; seedPos++)
+			if (ref.getMask(seedPos)) {
+				//forward
+				minL = gld->getMinL();
+				maxL = std::min(gld->getMaxL(), totLen - seedPos);
+				pfpos = seedPos;
+				for (int fragLen = minL; fragLen <= maxL; fragLen++) {
+					effL = std::min(fullLen, totLen - fragLen + 1);
+					factor = (mld == NULL ? 1.0 : mld->getAdjustedCumulativeProb(std::min(mld->getMaxL(), fragLen), fragLen));
+					value += probF * gld->getAdjustedProb(fragLen, totLen) * rspd->getAdjustedProb(pfpos, effL, fullLen) * factor;
+				}
+				//reverse
+				minL = gld->getMinL();
+				maxL = std::min(gld->getMaxL(), seedPos + seedLen);
+				for (int fragLen = minL; fragLen <= maxL; fragLen++) {
+					pfpos = seedPos - (fragLen - seedLen);
+					effL = std::min(fullLen, totLen - fragLen + 1);
+					factor = (mld == NULL ? 1.0 : mld->getAdjustedCumulativeProb(std::min(mld->getMaxL(), fragLen), fragLen));
+					value += probR * gld->getAdjustedProb(fragLen, totLen) * rspd->getAdjustedProb(pfpos, effL, fullLen) * factor;
+				}
+			}
     
-    mw[i] = 1.0 - value;
+		//for reverse strand masking
+		for (int seedPos = end; seedPos <= totLen - seedLen; seedPos++) {
+			minL = std::max(gld->getMinL(), seedPos + seedLen - fullLen + 1);
+			maxL = std::min(gld->getMaxL(), seedPos + seedLen);
+			for (int fragLen = minL; fragLen <= maxL; fragLen++) {
+				pfpos = seedPos - (fragLen - seedLen);
+				effL = std::min(fullLen, totLen - fragLen + 1);
+				factor = (mld == NULL ? 1.0 : mld->getAdjustedCumulativeProb(std::min(mld->getMaxL(), fragLen), fragLen));
+				value += probR * gld->getAdjustedProb(fragLen, totLen) * rspd->getAdjustedProb(pfpos, effL, fullLen) * factor;
+			}
+		}
 
-    if (mw[i] < 1e-8) { 
-      //      fprintf(stderr, "Warning: %dth reference sequence is masked for almost all positions!\n", i);
-      mw[i] = 0.0;
-    }
-  }
+		mw[i] = 1.0 - value;
+
+		if (mw[i] < 1e-8) {
+			//      fprintf(stderr, "Warning: %dth reference sequence is masked for almost all positions!\n", i);
+			mw[i] = 0.0;
+		}
+	}
 }
 
 #endif /* SINGLEMODEL_H_ */