v3.1.5

v 3.1.5
BackofenLab · Jan 31, 2020 · cee06fc · cee06fc
2 parents eff24d0 + 06325c3
commit cee06fc
Show file tree

Hide file tree

Showing 12 changed files with 215 additions and 103 deletions.
diff --git a/ChangeLog b/ChangeLog
@@ -10,10 +10,41 @@
 # changes in development version since last release 
 ################################################################################
 
+################################################################################  
+################################################################################
 
+################################################################################
+### version 3.1.5
 ################################################################################  
+
+# IntaRNA :
+- bugfix traceback of interactions with seed on right-boundary
+- bugfix traceback of seeds with bulges when outNoLP present
+
 ################################################################################
 
+200131 Martin Raden
+ * IntaRNA/PredictorMfe*SeedExtension* :
+   * undo bugfix
+   * traceback() :
+     + additional check if seed exceeds right interaction boundary
+ * IntaRNA/SeedHandler :
+   * isFeasibleSeedBasePair() :
+     * debug checks obsolete since part of the check
+   * updateToNextSeed() :
+     * bugfix : right boundary was exclusive (but has to be inclusive)
+ * IntaRNA/SeedHandlerMfe :
+   * getSeedE()
+   * setSeedE()
+   * traceBackSeed()
+     * now using global indices (offset shift done internally)
+   * fillSeed() :
+     + additional feasibility test for noLP
+   * traceback() :
+     * bugfix noLP energy trace : energy was from wrong cells
+     * trace small to large gaps (should be faster)
+ * test/SeedHandlerMfe :
+   + traceback tests (number of bps in traced seeds)
 
 ################################################################################
 ### version 3.1.4

diff --git a/configure.ac b/configure.ac
@@ -1,7 +1,7 @@
 
 AC_PREREQ([2.65])
 # 5 argument version only available with aclocal >= 2.64
-AC_INIT([IntaRNA], [3.1.4], [], [intaRNA], [http://www.bioinf.uni-freiburg.de] )
+AC_INIT([IntaRNA], [3.1.5], [], [intaRNA], [http://www.bioinf.uni-freiburg.de] )
 
 # minimal required version of the boost library 
 BOOST_REQUIRED_VERSION=1.50.0

diff --git a/doc/how-to-release-new-version.md b/doc/how-to-release-new-version.md
@@ -11,8 +11,12 @@
   - upload source code tar ball
   - (upload win binary)
   - (upload linux binary)
+	- (via cygwin) to get cygwin dlls used for compilation on Windows run
+	  - `for f in `ldd ./IntaRNA | grep "/usr/bin/" | awk '{print $3}'`; do cp $f .; done`
+	  - zip these together with `IntaRNA.exe` (and additional script files etc.)
   - (upload API docu pdf and html.zip)
 
+
 - publish release on github
 
 - update [bioconda recipe](https://github.com/bioconda/bioconda-recipes/tree/master/recipes/intarna)

diff --git a/src/IntaRNA/PredictorMfe2dHeuristicSeedExtension.cpp b/src/IntaRNA/PredictorMfe2dHeuristicSeedExtension.cpp
@@ -73,8 +73,8 @@ predict( const IndexRange & r1, const IndexRange & r2 )
 
 	size_t si1 = RnaSequence::lastPos, si2 = RnaSequence::lastPos;
 	while( seedHandler.updateToNextSeed(si1,si2
-			, 0, range_size1+1-seedHandler.getConstraint().getBasePairs()+seedHandler.getConstraint().getMaxUnpaired1()
-			, 0, range_size2+1-seedHandler.getConstraint().getBasePairs()+seedHandler.getConstraint().getMaxUnpaired2()) )
+			, 0, range_size1+1-seedHandler.getConstraint().getBasePairs()
+			, 0, range_size2+1-seedHandler.getConstraint().getBasePairs()) )
 	{
 		E_type seedE = seedHandler.getSeedE(si1, si2);
 		const size_t sl1 = seedHandler.getSeedLength1(si1, si2);

diff --git a/src/IntaRNA/PredictorMfe2dSeedExtension.cpp b/src/IntaRNA/PredictorMfe2dSeedExtension.cpp
@@ -74,8 +74,8 @@ predict( const IndexRange & r1, const IndexRange & r2 )
 
 	size_t si1 = RnaSequence::lastPos, si2 = RnaSequence::lastPos;
 	while( seedHandler.updateToNextSeed(si1,si2
-			, 0, range_size1+1-seedHandler.getConstraint().getBasePairs()+seedHandler.getConstraint().getMaxUnpaired1()
-			, 0, range_size2+1-seedHandler.getConstraint().getBasePairs()+seedHandler.getConstraint().getMaxUnpaired2()) )
+			, 0, range_size1+1-seedHandler.getConstraint().getBasePairs()
+			, 0, range_size2+1-seedHandler.getConstraint().getBasePairs()) )
 	{
 		// get energy and boundaries of seed
 		E_type seedE = seedHandler.getSeedE(si1, si2);
@@ -327,8 +327,8 @@ traceBack( Interaction & interaction )
 
 	size_t si1 = RnaSequence::lastPos, si2 = RnaSequence::lastPos;
 	while( seedHandler.updateToNextSeed(si1,si2
-			, i1,j1+1-seedHandler.getConstraint().getBasePairs()+seedHandler.getConstraint().getMaxUnpaired1()
-			, i2,j2+1-seedHandler.getConstraint().getBasePairs()+seedHandler.getConstraint().getMaxUnpaired2() ) )
+			, i1,j1+1-seedHandler.getConstraint().getBasePairs()
+			, i2,j2+1-seedHandler.getConstraint().getBasePairs() ) )
 	{
 			E_type seedE = seedHandler.getSeedE(si1, si2);
 			const size_t sl1 = seedHandler.getSeedLength1(si1, si2);
@@ -338,6 +338,11 @@ traceBack( Interaction & interaction )
 			const size_t maxMatrixLen1 = energy.getAccessibility1().getMaxLength()-sl1+1;
 			const size_t maxMatrixLen2 = energy.getAccessibility2().getMaxLength()-sl2+1;
 
+			// check if seed exceeds interaction (eg if with bulge)
+			if ( sj1 > j1 || sj2 > j2 ) {
+				continue;
+			}
+
 			hybridE_left.resize( std::min(si1+1, maxMatrixLen1), std::min(si2+1, maxMatrixLen2) );
 			fillHybridE_left( si1, si2 );
 			hybridE_right.resize( std::min(j1-sj1+1, maxMatrixLen1), std::min(j2-sj2+1, maxMatrixLen2) );

diff --git a/src/IntaRNA/PredictorMfe2dSeedExtensionRIblast.cpp b/src/IntaRNA/PredictorMfe2dSeedExtensionRIblast.cpp
@@ -78,8 +78,8 @@ predict( const IndexRange & r1, const IndexRange & r2  )
 
 	size_t si1 = RnaSequence::lastPos, si2 = RnaSequence::lastPos;
 	while( seedHandler.updateToNextSeed(si1,si2
-			, 0, interaction_size1+1-seedHandler.getConstraint().getBasePairs()+seedHandler.getConstraint().getMaxUnpaired1()
-			, 0, interaction_size2+1-seedHandler.getConstraint().getBasePairs()+seedHandler.getConstraint().getMaxUnpaired2()) )
+			, 0, interaction_size1+1-seedHandler.getConstraint().getBasePairs()
+			, 0, interaction_size2+1-seedHandler.getConstraint().getBasePairs()) )
 	{
 		E_type seedE = seedHandler.getSeedE(si1, si2);
 		size_t sl1 = seedHandler.getSeedLength1(si1, si2);
@@ -373,8 +373,8 @@ traceBack( Interaction & interaction )
 
 	size_t si1 = RnaSequence::lastPos, si2 = RnaSequence::lastPos;
 	while( seedHandler.updateToNextSeed(si1,si2
-			, i1,j1+1-seedHandler.getConstraint().getBasePairs()+seedHandler.getConstraint().getMaxUnpaired1()
-			, i2,j2+1-seedHandler.getConstraint().getBasePairs()+seedHandler.getConstraint().getMaxUnpaired2() ) )
+			, i1,j1+1-seedHandler.getConstraint().getBasePairs()
+			, i2,j2+1-seedHandler.getConstraint().getBasePairs() ) )
 	{
 		E_type seedE = seedHandler.getSeedE(si1, si2);
 
@@ -383,6 +383,11 @@ traceBack( Interaction & interaction )
 		size_t sj1 = si1+sl1-1;
 		size_t sj2 = si2+sl2-1;
 
+		// check if seed exceeds interaction (eg if with bulge)
+		if ( sj1 > j1 || sj2 > j2 ) {
+			continue;
+		}
+
 		ExtendedSeed extension;
 		extension.i1 = si1;
 		extension.i2 = si2;

diff --git a/src/IntaRNA/PredictorMfeEns2dHeuristicSeedExtension.cpp b/src/IntaRNA/PredictorMfeEns2dHeuristicSeedExtension.cpp
@@ -75,8 +75,8 @@ predict( const IndexRange & r1, const IndexRange & r2 )
 
 	size_t si1 = RnaSequence::lastPos, si2 = RnaSequence::lastPos;
 	while( seedHandler.updateToNextSeed(si1,si2
-			, 0, range_size1+1-seedHandler.getConstraint().getBasePairs()+seedHandler.getConstraint().getMaxUnpaired1()
-			, 0, range_size2+1-seedHandler.getConstraint().getBasePairs()+seedHandler.getConstraint().getMaxUnpaired2()) )
+			, 0, range_size1+1-seedHandler.getConstraint().getBasePairs()
+			, 0, range_size2+1-seedHandler.getConstraint().getBasePairs()) )
 	{
 		const Z_type seedZ = energy.getBoltzmannWeight( seedHandler.getSeedE(si1, si2) );
 		const size_t sl1 = seedHandler.getSeedLength1(si1, si2);

diff --git a/src/IntaRNA/PredictorMfeEns2dSeedExtension.cpp b/src/IntaRNA/PredictorMfeEns2dSeedExtension.cpp
@@ -76,8 +76,8 @@ predict( const IndexRange & r1, const IndexRange & r2 )
 
 	size_t si1 = RnaSequence::lastPos, si2 = RnaSequence::lastPos;
 	while( seedHandler.updateToNextSeed(si1,si2
-			, 0, range_size1+1-seedHandler.getConstraint().getBasePairs()+seedHandler.getConstraint().getMaxUnpaired1()
-			, 0, range_size2+1-seedHandler.getConstraint().getBasePairs()+seedHandler.getConstraint().getMaxUnpaired2()) )
+			, 0, range_size1+1-seedHandler.getConstraint().getBasePairs()
+			, 0, range_size2+1-seedHandler.getConstraint().getBasePairs()) )
 	{
 		// get Z and boundaries of seed
 		const Z_type seedZ = energy.getBoltzmannWeight( seedHandler.getSeedE(si1, si2) );

diff --git a/src/IntaRNA/SeedHandler.cpp b/src/IntaRNA/SeedHandler.cpp
@@ -11,10 +11,6 @@ bool
 SeedHandler::
 isFeasibleSeedBasePair( const size_t i1, const size_t i2, const bool atEndOfSeed ) const
 {
-#if INTARNA_IN_DEBUG_MODE
-	if ( i1 >= energy.size1() ) throw std::runtime_error("SeedHandler::isFeasibleSeedBasePair: i1("+toString(i1)+") >= energy.size1("+toString(energy.size1())+")");
-	if ( i2 >= energy.size2() ) throw std::runtime_error("SeedHandler::isFeasibleSeedBasePair: i2("+toString(i2)+") >= energy.size2("+toString(energy.size2())+")");
-#endif
 
 	return		i1 < energy.size1() && i2 < energy.size2()
 			&&	energy.isAccessible1(i1)
@@ -52,23 +48,23 @@ updateToNextSeed( size_t & i1_out, size_t & i2_out
 		i2 = i2min;
 	} else {
 		// update to next potential seed position
-		if (++i1 >= i1maxVal) {
+		if (++i1 > i1maxVal) {
 			i1 = i1min;
 			i2++;
 		}
 	}
 
 	// find next valid seed start within range
-	while( i2 < i2maxVal && !(isSeedBound(i1,i2))) {
+	while( i2 <= i2maxVal && !(isSeedBound(i1,i2))) {
 		// update seed position within range
-		if (++i1 == i1maxVal) {
+		if (++i1 > i1maxVal) {
 			i1 = i1min;
 			i2++;
 		}
 	}
 
 	// check if we found a valid seed in the range
-	if (i1 < i1maxVal && i2< i2maxVal) {
+	if (i1 <= i1maxVal && i2 <= i2maxVal) {
 		i1_out = i1;
 		i2_out = i2;
 		return true;

diff --git a/src/IntaRNA/SeedHandlerMfe.cpp b/src/IntaRNA/SeedHandlerMfe.cpp
@@ -64,6 +64,20 @@ fillSeed( const size_t i1min, const size_t i1max, const size_t i2min, const size
 		// bp=0 encodes 2 base pairs
 		for (bpIn=0; bpIn<seedE_rec.shape()[2] && (i1+bpIn+1-offset1)<seed.size1() && (i2+bpIn+1-offset2)<seed.size2(); bpIn++) {
 
+			bool validLeftEnd = true;
+
+			// if no LP : check for direct right-stacking of i
+			if (noLpShift > 0) {
+				// check if feasible extension
+				if (isFeasibleSeedBasePair(i1+1,i2+1)) {
+					// get stacking energy
+					iStackE = energy.getE_interLeft(i1,i1+1,i2,i2+1);
+				} else {
+					// no valid noLP extension possible
+					validLeftEnd = false;
+				}
+			}
+
 			// for feasible unpaired in seq1 in increasing order
 			for (u1=0; u1<seedE_rec.shape()[3] && (i1+bpIn+1+u1-offset1) < seed.size1(); u1++) {
 
@@ -76,26 +90,16 @@ fillSeed( const size_t i1min, const size_t i1max, const size_t i2min, const size
 				// get right seed boundaries
 				j1 = i1+bpIn+1+u1;
 				j2 = i2+bpIn+1+u2;
-				// check if right boundary is complementary and feasible
-				// check if this index range is to be considered for seed search
-				bool validSeedSite = isFeasibleSeedBasePair(j1,j2,true);
-
-				// if no LP : check for direct right-stacking of i
-				if (noLpShift > 0) {
-					// check if feasible extension
-					if (isFeasibleSeedBasePair(i1+1,i2+1)) {
-						// get stacking energy
-						iStackE = energy.getE_interLeft(i1,i1+1,i2,i2+1);
-					} else {
-						validSeedSite = false;
-					}
-				}
 
 				// init current seed energy
 				curE = E_INF;
 
-				// check if right boundary is complementary
-				if (validSeedSite) {
+				// check if this index range is to be considered for seed search
+				// check if boundaries are complementary
+				if (validLeftEnd
+					&& isFeasibleSeedBasePair(j1,j2,true)
+					&& (noLpShift==0 || isFeasibleSeedBasePair(j1-1,j2-1)))
+				{
 
 					// base case: only left and right base pair present
 					if (bpIn==0) {
@@ -108,8 +112,8 @@ fillSeed( const size_t i1min, const size_t i1max, const size_t i2min, const size
 					} else {
 
 						// explicitly check direct stacking extension in noLP mode
-						if (noLpShift > 0 && E_isNotINF( getSeedE( i1+1-offset1, i2+1-offset2, bpIn-1, u1, u2 ) )) {
-							curE = std::min( curE, iStackE + getSeedE( i1+1-offset1, i2+1-offset2, bpIn-1, u1, u2 ) );
+						if (noLpShift > 0 && E_isNotINF( getSeedE( i1+1, i2+1, bpIn-1, u1, u2 ) )) {
+							curE = std::min( curE, iStackE + getSeedE( i1+1, i2+1, bpIn-1, u1, u2 ) );
 						}
 
 						// if enough interior base pairs left
@@ -127,15 +131,15 @@ fillSeed( const size_t i1min, const size_t i1max, const size_t i2min, const size
 								k2 = i2+u2p+1+noLpShift;
 								// check if split pair is complementary
 								// and recursed entry is < E_INF
-								if (! (isFeasibleSeedBasePair(k1,k2) && E_isNotINF( getSeedE( k1-offset1, k2-offset2, bpIn-1-noLpShift, u1-u1p, u2-u2p ) ) ) ) {
+								if (! (isFeasibleSeedBasePair(k1,k2) && E_isNotINF( getSeedE( k1, k2, bpIn-1-noLpShift, u1-u1p, u2-u2p ) ) ) ) {
 									continue; // not complementary -> skip
 								}
 
 								// update mfe for split at k1,k2
 								curE = std::min( curE,
 										iStackE
 										+ energy.getE_interLeft(i1+noLpShift,k1,i2+noLpShift,k2)
-										+ getSeedE( k1-offset1, k2-offset2, bpIn-1-noLpShift, u1-u1p, u2-u2p )
+										+ getSeedE( k1, k2, bpIn-1-noLpShift, u1-u1p, u2-u2p )
 										);
 							} // u2p
 							} // u1p
@@ -145,8 +149,7 @@ fillSeed( const size_t i1min, const size_t i1max, const size_t i2min, const size
 				} // (j1,j2) complementary
 
 				// store seed energy
-				setSeedE( i1-offset1, i2-offset2, bpIn, u1, u2, curE );
-
+				setSeedE( i1, i2, bpIn, u1, u2, curE );
 			} // u2
 			} // u1
 
@@ -176,7 +179,7 @@ fillSeed( const size_t i1min, const size_t i1max, const size_t i2min, const size
 					}
 
 					// get overall interaction energy
-					E_type curEhyb = getSeedE( i1-offset1, i2-offset2, bpIn, u1, u2 ) + energy.getE_init();
+					E_type curEhyb = getSeedE( i1, i2, bpIn, u1, u2 ) + energy.getE_init();
 					curE = energy.getE( i1, j1, i2, j2, curEhyb );
 
 					// check hybrid energy bound including Einit
@@ -195,7 +198,7 @@ fillSeed( const size_t i1min, const size_t i1max, const size_t i2min, const size
 				// reduce bestE to hybridization energy only (init+loops)
 				if (E_isNotINF( bestE )) {
 					// get seed's hybridization loop energies only
-					bestE = getSeedE( i1-offset1, i2-offset2, bpIn, u1best, u2best );
+					bestE = getSeedE( i1, i2, bpIn, u1best, u2best );
 					// count true seed
 					seedCountNotInf++;
 				}
@@ -256,13 +259,11 @@ traceBackSeed( Interaction & interaction
 	// trace each seed base pair (excluding right most)
 	for( size_t bpIn=1+bpInbetween; bpIn-- > 0; ) {
 
-
-
 		// base case: only left and right base pair present
 		if (bpIn==0) {
 			// add left base pair if not left seed boundary
 			if (i1 != i1_) {
-				interaction.basePairs.push_back( energy.getBasePair(i1+offset1,i2+offset2) );
+				interaction.basePairs.push_back( energy.getBasePair(i1,i2) );
 			}
 
 		} else {
@@ -272,12 +273,16 @@ traceBackSeed( Interaction & interaction
 				// check if feasible extension
 				assert(isFeasibleSeedBasePair(i1+1,i2+1));
 				// get stacking energy
-				iStackE = energy.getE_interLeft(i1+offset1,i1+1+offset1,i2+offset2,i2+1+offset2);
+				iStackE = energy.getE_interLeft(i1,i1+1,i2,i2+1);
 				// noLP : check stacking of i
 				if ( E_equal( curE, iStackE + getSeedE( i1+1, i2+1, bpIn-1, u1max, u2max )) ) {
+					// store left base pair if not left seed boundary
+					if (i1 != i1_) {
+						interaction.basePairs.push_back( energy.getBasePair(i1,i2) );
+					}
 					i1++;
 					i2++;
-					curE = getSeedE( i1+1, i2+1, bpIn-1, u1max, u2max );
+					curE = getSeedE( i1, i2, bpIn-1, u1max, u2max );
 					continue;
 				}
 				// sanity check for noLP mode
@@ -288,36 +293,36 @@ traceBackSeed( Interaction & interaction
 			// i1 .. i1+u1p+1 .. j1
 			// i2 .. i2+u2p+1 .. j2
 			bool traceNotFound = true;
-			for (u1=1+u1max; traceNotFound && u1-- > 0;) {
-			for (u2=1+u2max; traceNotFound && u2-- > 0;) {
+			for (u1=0; traceNotFound && u1<=u1max ; u1++) {
+			for (u2=0; traceNotFound && u2<=u2max && (u1+u2)<=uMax ; u2++) {
 				// check if overall number of unpaired is not exceeded
 				// or skip stacked extension since covered above
-				if (u1+u2 > uMax || u1+u2 < noLpShift) {
+				if (u1+u2 < noLpShift) {
 					continue;
 				}
 
 				k1 = i1+u1+1+noLpShift;
 				k2 = i2+u2+1+noLpShift;
 
 				// check if valid trace
-				if ( isFeasibleSeedBasePair(k1+offset1, k2+offset2) && E_isNotINF( getSeedE( k1, k2, bpIn-1, u1max-u1, u2max-u2 ) ) ) {
+				if ( isFeasibleSeedBasePair(k1, k2) && E_isNotINF( getSeedE( k1, k2, bpIn-1-noLpShift, u1max-u1, u2max-u2 ) ) ) {
 
 					// check if correct trace
 					if ( E_equal( curE, iStackE
-										+ energy.getE_interLeft(i1+noLpShift+offset1,k1+offset1,i2+noLpShift+offset2,k2+offset2)
+										+ energy.getE_interLeft(i1+noLpShift,k1,i2+noLpShift,k2)
 										+ getSeedE( k1, k2, bpIn-1-noLpShift, u1max-u1, u2max-u2 )) )
 					{
+						// store next energy value to trace
+						curE = getSeedE( k1, k2, bpIn-1-noLpShift, u1max-u1, u2max-u2 );
 						// store left base pair if not left seed boundary
 						if (i1 != i1_) {
-							interaction.basePairs.push_back( energy.getBasePair(i1+offset1,i2+offset2) );
+							interaction.basePairs.push_back( energy.getBasePair(i1,i2) );
 						}
 						if (noLpShift > 0) {
-							interaction.basePairs.push_back( energy.getBasePair(i1+noLpShift+offset1,i2+noLpShift+offset2) );
+							interaction.basePairs.push_back( energy.getBasePair(i1+noLpShift,i2+noLpShift) );
 							// reflect additional base pair
 							bpIn--;
 						}
-						// store next energy value to trace
-						curE = getSeedE( k1, k2, bpIn-1-noLpShift, u1max-u1, u2max-u2 );
 						// reset for next trace step
 						i1 = k1;
 						i2 = k2;