BUGFIX : dangling-end probability computation using maxLength+1 ED values

Author: martin-raden

ChangeLog

@@ -22,6 +22,7 @@
 - BUGFIX : target accessibility limits and seed-free outMinPu filtering
 - BUGFIX : accessibility range splitting and intramolecular ensemble energies
 - BUGFIX : interaction assignment, seed comparison and output handling
+- BUGFIX : dangling-end probability computation using maxLength+1 ED values
 - add exhaustive interaction and seed oracles for short-sequence tests
 
 - reduced memory use for exact non-seed ensemble prediction without trackers
@@ -44,6 +45,12 @@
  * src/IntaRNA/Makefile.am :
  * src/bin/Makefile.am :
    * undo some diagnostic flags to allow compilation on test system
+ * IntaRNA/Accessibility* :
+   * constructors: generate/read/keep data for maxLength+1 to allow for 
+     dangling-end probability computation
+ * tests/AccessibilityFromStream_test.cpp :
+ * tests/AccessibilityVrna_test.cpp :
+   + check for maxLength+1 ED values
 
 260919 Alexander Mitrofanov
  * IntaRNA/PredictorMfe2d* :

src/IntaRNA/Accessibility.cpp

@@ -30,7 +30,7 @@ operator<<(std::ostream& out, const Accessibility& acc)
 		// print first (without delimiter)
 		out <<acc.getED(i, i);
 		// print remaining with delimiter
-		for (size_t j=i+1; j<std::min(acc.getMaxLength()+i,acc.getSequence().size());j++) {
+		for (size_t j=i+1; j<std::min(acc.getMaxLength()+1+i,acc.getSequence().size());j++) { // "getMaxLength()+1" to allow for dangling-end probability computation
 			out <<delimiter <<(acc.getED(i, j)>=0?" ":"") <<acc.getED(i, j);
 		}
 		out <<"\n";
@@ -67,7 +67,7 @@ writeRNAplfold_text( std::ostream& out, const Z_type RT, const bool writeProbs )
 
 	// length information
 	out <<" #i$	l=";
-	for (size_t l=1; l<=std::min(getSequence().size(),getMaxLength()); l++) {
+	for (size_t l=1; l<=std::min(getSequence().size(),getMaxLength()+1); l++) { // "getMaxLength()+1" to allow for dangling-end probability computation
 		out <<l <<'\t';
 	}
 	out <<'\n';
@@ -76,7 +76,7 @@ writeRNAplfold_text( std::ostream& out, const Z_type RT, const bool writeProbs )
 	for (size_t j=0; j<getSequence().size(); j++) {
 		// print end of window = j
 		out <<(j+1) <<'\t';
-		size_t maxL = std::min(j+1,getMaxLength());
+		size_t maxL = std::min(j+1,std::min(getSequence().size(),getMaxLength()+1)); // "getMaxLength()+1" to allow for dangling-end probability computation
 		// for each (increasing) window length, print value
 		for ( size_t l = 1; l <= maxL; l++ ) {
 			if (writeProbs) {
@@ -100,7 +100,7 @@ writeRNAplfold_text( std::ostream& out, const Z_type RT, const bool writeProbs )
 			}
 		}
 		// print NA for remaining entries
-		for ( size_t l = maxL+1; l<=getMaxLength(); l++ ) {
+		for ( size_t l = maxL+1; l<=std::min(getSequence().size(),getMaxLength()+1); l++ ) { // "getMaxLength()+1" to allow for dangling-end probability computation
 			out <<"NA\t";
 		}
 		// line break

src/IntaRNA/AccessibilityBasePair.cpp

@@ -70,10 +70,7 @@ AccessibilityBasePair::~AccessibilityBasePair()
 E_type
 AccessibilityBasePair::getED( const size_t from, const size_t to ) const
 {
-  if (from > to || to < 0 || from >= seq.size()) {
-    throw std::runtime_error( "AccessibilityBasePair::getED("+toString(from)+","+toString(to)
-                             +") : Arguments must satisfy 0 <= from <= to < seq.length" );
-  }
+  checkIndices(from,to);
   return logPu(from, to);
 };
 

src/IntaRNA/AccessibilityDisabled.h

@@ -46,7 +46,7 @@ class AccessibilityDisabled: public Accessibility {
 	 * @param from the start index of the regions (from <= to)
 	 * @param to the end index of the regions (to <= seq.length())
 	 *
-	 * @return 0  if (j-1+1) <= maxLength or ED_UPPER_BOUND otherwise
+	 * @return 0  if (j-i+1) <= maxLength or ED_UPPER_BOUND otherwise
 	 */
 	virtual
 	E_type
@@ -88,7 +88,7 @@ getED( const size_t from, const size_t to ) const
 	// input check
 	checkIndices(from,to);
 
-	if ((to-from+1) <= getMaxLength()) {
+	if ((to-from+1) <= getMaxLength()+1) { // "getMaxLength()+1" to allow for dangling-end probability computation
 		// check for constrained end positions
 		if (!getAccConstraint().isAccessible(from) || !getAccConstraint().isAccessible(to)) {
 			// end position blocked --> omit accessibility

src/IntaRNA/AccessibilityFromStream.h

@@ -139,7 +139,7 @@ getED( const size_t from, const size_t to ) const
 	// input range check
 	checkIndices(from,to);
 
-	if ((to-from+1) <= getMaxLength()) {
+	if ((to-from+1) <= getMaxLength()+1) { // "getMaxLength()+1" to allow for dangling-end probability computation
 		// check for constrained end positions
 		if (!getAccConstraint().isAccessible(from) || !getAccConstraint().isAccessible(to)) {
 			// end position blocked --> omit accessibility

src/IntaRNA/AccessibilityVrna.cpp

@@ -83,7 +83,8 @@ AccessibilityVrna::AccessibilityVrna(
 		)
  :
 	Accessibility( seq, maxLength, accConstraint ),
-	edValues( getSequence().size(), getSequence().size(), 0, getMaxLength() )
+	edValues( getSequence().size(), getSequence().size(), 
+				0, std::min(getSequence().size(),getMaxLength()+1) ) // +1 to allow for dangling-end probability computation
 {
 	// if sequence shows minimal length
 	if (seq.size() > 4) {
@@ -134,7 +135,7 @@ callbackForStorage(FLT_OR_DBL   *pr,
 	    // copy unpaired data for all available interval lengths
 	    // but ensure interval does not contain blocked positions
 	    const bool rightEndBlocked = accConstr.isMarkedBlocked(j-1);
-	    for (int l = std::min(j,std::min(pr_size,std::min(max,(int)storageRT.first->getMaxLength()))); l>=1; l--) {
+	    for (int l = std::min(j,std::min(pr_size,std::min(max,(int)storageRT.first->getMaxLength()+1))); l>=1; l--) { // "getMaxLength()+1" to allow for dangling-end probability computation
 			// get unpaired probability
 			FLT_OR_DBL prob_unpaired = pr[l];
 //			TODO: check for [0,1] range and correct if needed (print WARNING)

src/IntaRNA/AccessibilityVrna.h

@@ -149,7 +149,7 @@ getED( const size_t from, const size_t to ) const
 	// input range check
 	checkIndices(from,to);
 
-	if ((to-from+1) <= getMaxLength()) {
+	if ((to-from+1) <= getMaxLength() +1 ) { // "getMaxLength()+1" to allow for dangling-end probability computation
 		// check for constrained end positions
 		if (!getAccConstraint().isAccessible(from) || !getAccConstraint().isAccessible(to)) {
 			// end position blocked --> omit accessibility

tests/.gitignore

@@ -1,2 +1,3 @@
 /*.log
 /*.trs
+runApiTests

tests/AccessibilityFromStream_test.cpp

@@ -67,9 +67,13 @@ TEST_CASE( "AccessibilityFromStream", "[AccessibilityFromStream]" ) {
 //		std::cerr <<"orig data:\n" <<accString;
 //		std::cerr <<"ED data:\n" <<acc;
 
+		// check max length
+		REQUIRE( acc.getMaxLength() == 9 );
+
 		// check elements
-		REQUIRE( acc.getED(29, 29) == 0 );
-		REQUIRE( acc.getED(21, 29) == 690 );
+		REQUIRE( acc.getED(29, 29) == 0 );   // length 1 (0.9980056 -> -log(0.9980056)*100 = 0.19944 -> 0)
+		REQUIRE( acc.getED(21, 29) == 690 ); // length 9 (0.001002708 -> -log(0.001002708)*100 = 690.03 -> 690)
+		REQUIRE( acc.getED(20, 29) == 732 ); // check for length exceeding maxLength+1 used in dangling-end probability computation
 //		// old checks not working for integer-based ED type
 //		REQUIRE( std::exp( - acc.getED(29, 29) ) > 0.998 );
 //		REQUIRE( std::exp( - acc.getED(29, 29) ) < 0.999 );

tests/AccessibilityVrna_test.cpp

@@ -43,6 +43,16 @@ TEST_CASE("AccessibilityVrna", "[AccessibilityVrna]") {
 		// 'x' is passed at both C-string boundaries and forces unpaired bases.
 		REQUIRE( E_equal(acc.getED(0, 0), 0) );
 		REQUIRE( E_equal(acc.getED(rna.size()-1, rna.size()-1), 0) );
+		// 'p' is passed at both C-string boundaries and forces paired bases.
+		REQUIRE( E_equal(acc.getED(1, 1), Accessibility::ED_UPPER_BOUND) );
+		REQUIRE( E_equal(acc.getED(rna.size()-2, rna.size()-2), Accessibility::ED_UPPER_BOUND) );
+		
+		// check maximal length
+		REQUIRE( acc.getMaxLength() == 1 );
+		// check for ED values of length+1 regions required for dangling-end probability computation
+		REQUIRE( E_equal(acc.getED(0, 1), Accessibility::ED_UPPER_BOUND) ); // not accessible due to 'p' at position 1
+		//std::cerr << "acc.getED(2,3) = " << acc.getED(2,3) << std::endl;
+		REQUIRE( E_equal(acc.getED(2,3), 610) ); // should be accessible for dangling-end probability computation of position 2 or 3
 
 		// Exercise the same x|........|x bytes through computeES() and
 		// computeIntraEall(), including the reversed accessibility path.