From d994e6617e9381237dbdb628fbacc6cd1ddeebc8 Mon Sep 17 00:00:00 2001 From: martin-mann Date: Fri, 31 Jan 2020 10:23:46 +0100 Subject: [PATCH 1/8] win binary zip generation --- doc/how-to-release-new-version.md | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/doc/how-to-release-new-version.md b/doc/how-to-release-new-version.md index 0b0f4063..281df428 100644 --- a/doc/how-to-release-new-version.md +++ b/doc/how-to-release-new-version.md @@ -13,6 +13,10 @@ - (upload linux binary) - (upload API docu pdf and html.zip) +- (win binary) 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.) + - publish release on github - update [bioconda recipe](https://github.com/bioconda/bioconda-recipes/tree/master/recipes/intarna) From 6c166d84dc1dc1c093e8c0851252c9ded3f93aad Mon Sep 17 00:00:00 2001 From: martin-mann Date: Fri, 31 Jan 2020 10:37:20 +0100 Subject: [PATCH 2/8] bump dev version --- configure.ac | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/configure.ac b/configure.ac index 41871bab..1b48235f 100644 --- 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 From d91992a42de9e29739532cc3992170d0e2f475c1 Mon Sep 17 00:00:00 2001 From: martin-mann Date: Fri, 31 Jan 2020 10:39:02 +0100 Subject: [PATCH 3/8] restructure --- doc/how-to-release-new-version.md | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/doc/how-to-release-new-version.md b/doc/how-to-release-new-version.md index 281df428..8bb0cce6 100644 --- a/doc/how-to-release-new-version.md +++ b/doc/how-to-release-new-version.md @@ -11,11 +11,11 @@ - 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) -- (win binary) 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.) - publish release on github From ca3c0664bfaab446ec98cf24bda44a3dc8b8037f Mon Sep 17 00:00:00 2001 From: martin-mann Date: Fri, 31 Jan 2020 16:13:19 +0100 Subject: [PATCH 4/8] SeedHandler.cpp --- src/IntaRNA/SeedHandler.cpp | 12 ++++-------- 1 file changed, 4 insertions(+), 8 deletions(-) diff --git a/src/IntaRNA/SeedHandler.cpp b/src/IntaRNA/SeedHandler.cpp index 20358fe1..ed1196eb 100644 --- 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; From fd839bb22bbc11a9bd9bbf23a21ad4023e84358f Mon Sep 17 00:00:00 2001 From: martin-mann Date: Fri, 31 Jan 2020 16:13:30 +0100 Subject: [PATCH 5/8] * IntaRNA/PredictorMfe*SeedExtension* : * undo bugfix * traceback() : + additional check if seed exceeds right interaction boundary --- .../PredictorMfe2dHeuristicSeedExtension.cpp | 4 ++-- src/IntaRNA/PredictorMfe2dSeedExtension.cpp | 13 +++++++++---- src/IntaRNA/PredictorMfe2dSeedExtensionRIblast.cpp | 13 +++++++++---- .../PredictorMfeEns2dHeuristicSeedExtension.cpp | 4 ++-- src/IntaRNA/PredictorMfeEns2dSeedExtension.cpp | 4 ++-- 5 files changed, 24 insertions(+), 14 deletions(-) diff --git a/src/IntaRNA/PredictorMfe2dHeuristicSeedExtension.cpp b/src/IntaRNA/PredictorMfe2dHeuristicSeedExtension.cpp index 58c2035f..4e124e0f 100644 --- 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 index a178b895..38521f2b 100644 --- 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 index 864e652a..f150ff4c 100644 --- 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 index 21523ea6..bc1fe839 100644 --- 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 index 4e5da9bb..cc7aa459 100644 --- 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) ); From d2681f6a446bfbc007ab340519c9297aadae1d51 Mon Sep 17 00:00:00 2001 From: martin-mann Date: Fri, 31 Jan 2020 16:20:35 +0100 Subject: [PATCH 6/8] + traceback tests (number of bps in traced seeds --- tests/SeedHandlerMfe_test.cpp | 118 ++++++++++++++++++++++++---------- 1 file changed, 85 insertions(+), 33 deletions(-) diff --git a/tests/SeedHandlerMfe_test.cpp b/tests/SeedHandlerMfe_test.cpp index 6f5487cc..f78ffa4a 100644 --- a/tests/SeedHandlerMfe_test.cpp +++ b/tests/SeedHandlerMfe_test.cpp @@ -64,6 +64,22 @@ TEST_CASE( "SeedHandlerMfe", "[SeedHandlerMfe]") { ///////////////////////////////////////////// FILLSEED //////////////////////////////////////////////// ////////////////////////////////////////////////////////////////////////////////////////////////////////// REQUIRE(sHM.fillSeed(0,energy.size1()-1, 0,energy.size2()-1) == 4); // getting only one seed per si, thus 2x2 + + ////////////////////////////////////////////////////////////////////////////////////////////////////////// + ///////////////////////////////////////////// TRACEBACK /////////////////////////////////////////////// + ////////////////////////////////////////////////////////////////////////////////////////////////////////// + size_t si1=9999999, si2=9999999; + while( sHM.updateToNextSeed(si1,si2,0,r1.size()-1,0,r2.size()) ) { + // prepare trace + Interaction i(r1,r2); + i.basePairs.push_back(energy.getBasePair(si1,si2)); + i.basePairs.push_back(energy.getBasePair(si1+sHM.getSeedLength1(si1,si2)-1,si2+sHM.getSeedLength2(si1,si2)-1)); + // trace base pairs + sHM.traceBackSeed(i, si1, si2); + // check number of traced base pairs + REQUIRE(i.basePairs.size() == 4); + } + } SECTION("SeedHandlerMfe: no lp - 1", "[SeedHandlerMfe]") { @@ -75,42 +91,59 @@ TEST_CASE( "SeedHandlerMfe", "[SeedHandlerMfe]") { InteractionEnergyBasePair energy(acc1, racc); { - size_t bp = 2; - // seedBP / seedMaxUP / seedTMaxUP / seedQMaxUP / seedMaxE / seedMaxED / seedTRange / seedQRange / seedTQ - SeedConstraint sC(bp,3,3,0,0 - , AccessibilityDisabled::ED_UPPER_BOUND - , 0 - , IndexRangeList("") - , IndexRangeList("") - , "" - , false, false, true - ); - - SeedHandlerMfe sHM(energy, sC); - - ////////////////////////////////////////////////////////////////////////////////////////////////////////// - ///////////////////////////////////////////// FILLSEED //////////////////////////////////////////////// - ////////////////////////////////////////////////////////////////////////////////////////////////////////// - REQUIRE(sHM.fillSeed(0,energy.size1()-1, 0,energy.size2()-1) == 4); + size_t bp = 2; + // seedBP / seedMaxUP / seedTMaxUP / seedQMaxUP / seedMaxE / seedMaxED / seedTRange / seedQRange / seedTQ + SeedConstraint sC(bp,3,3,0,0 + , AccessibilityDisabled::ED_UPPER_BOUND + , 0 + , IndexRangeList("") + , IndexRangeList("") + , "" + , false, false, true + ); + + SeedHandlerMfe sHM(energy, sC); + + ////////////////////////////////////////////////////////////////////////////////////////////////////////// + ///////////////////////////////////////////// FILLSEED //////////////////////////////////////////////// + ////////////////////////////////////////////////////////////////////////////////////////////////////////// + REQUIRE(sHM.fillSeed(0,energy.size1()-1, 0,energy.size2()-1) == 4); } for (size_t bp = 3; bp <= 4; bp++) { - // seedBP / seedMaxUP / seedTMaxUP / seedQMaxUP / seedMaxE / seedMaxED / seedTRange / seedQRange / seedTQ - SeedConstraint sC(bp,3,3,0,0 - , AccessibilityDisabled::ED_UPPER_BOUND - , 0 - , IndexRangeList("") - , IndexRangeList("") - , "" - , false, false, true - ); + // seedBP / seedMaxUP / seedTMaxUP / seedQMaxUP / seedMaxE / seedMaxED / seedTRange / seedQRange / seedTQ + SeedConstraint sC(bp,3,3,0,0 + , AccessibilityDisabled::ED_UPPER_BOUND + , 0 + , IndexRangeList("") + , IndexRangeList("") + , "" + , false, false, true + ); + + SeedHandlerMfe sHM(energy, sC); + + ////////////////////////////////////////////////////////////////////////////////////////////////////////// + ///////////////////////////////////////////// FILLSEED //////////////////////////////////////////////// + ////////////////////////////////////////////////////////////////////////////////////////////////////////// + REQUIRE(sHM.fillSeed(0,energy.size1()-1, 0,energy.size2()-1) == 0); + + + ////////////////////////////////////////////////////////////////////////////////////////////////////////// + ///////////////////////////////////////////// TRACEBACK /////////////////////////////////////////////// + ////////////////////////////////////////////////////////////////////////////////////////////////////////// + size_t si1=9999999, si2=9999999; + while( sHM.updateToNextSeed(si1,si2,0,r1.size()-1,0,r2.size()) ) { + // prepare trace + Interaction i(r1,r2); + i.basePairs.push_back(energy.getBasePair(si1,si2)); + i.basePairs.push_back(energy.getBasePair(si1+sHM.getSeedLength1(si1,si2)-1,si2+sHM.getSeedLength2(si1,si2)-1)); + // trace base pairs + sHM.traceBackSeed(i, si1, si2); + // check number of traced base pairs + REQUIRE(i.basePairs.size() == bp); + } - SeedHandlerMfe sHM(energy, sC); - - ////////////////////////////////////////////////////////////////////////////////////////////////////////// - ///////////////////////////////////////////// FILLSEED //////////////////////////////////////////////// - ////////////////////////////////////////////////////////////////////////////////////////////////////////// - REQUIRE(sHM.fillSeed(0,energy.size1()-1, 0,energy.size2()-1) == 0); } } @@ -122,8 +155,9 @@ TEST_CASE( "SeedHandlerMfe", "[SeedHandlerMfe]") { ReverseAccessibility racc(acc2); InteractionEnergyBasePair energy(acc1, racc); + size_t seedBP=4; // seedBP / seedMaxUP / seedTMaxUP / seedQMaxUP / seedMaxE / seedMaxED / seedTRange / seedQRange / seedTQ - SeedConstraint sC(4,3,3,0,0 + SeedConstraint sC(seedBP,3,3,0,0 , AccessibilityDisabled::ED_UPPER_BOUND , 0 , IndexRangeList("") @@ -138,5 +172,23 @@ TEST_CASE( "SeedHandlerMfe", "[SeedHandlerMfe]") { ///////////////////////////////////////////// FILLSEED //////////////////////////////////////////////// ////////////////////////////////////////////////////////////////////////////////////////////////////////// REQUIRE(sHM.fillSeed(0,energy.size1()-1, 0,energy.size2()-1) == 2); + + + ////////////////////////////////////////////////////////////////////////////////////////////////////////// + ///////////////////////////////////////////// TRACEBACK /////////////////////////////////////////////// + ////////////////////////////////////////////////////////////////////////////////////////////////////////// + size_t si1=9999999, si2=9999999; + while( sHM.updateToNextSeed(si1,si2,0,r1.size()-1,0,r2.size()) ) { + // prepare trace + Interaction i(r1,r2); + i.basePairs.push_back(energy.getBasePair(si1,si2)); + i.basePairs.push_back(energy.getBasePair(si1+sHM.getSeedLength1(si1,si2)-1,si2+sHM.getSeedLength2(si1,si2)-1)); + // trace base pairs + sHM.traceBackSeed(i, si1, si2); + // check number of traced base pairs + REQUIRE(i.basePairs.size() == seedBP); + } } + + } From 44f16f2590d6308cd9bf8f2ee462df3025e097b3 Mon Sep 17 00:00:00 2001 From: martin-mann Date: Fri, 31 Jan 2020 16:20:49 +0100 Subject: [PATCH 7/8] * 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) --- ChangeLog | 25 +++++++++++ src/IntaRNA/SeedHandlerMfe.cpp | 81 ++++++++++++++++++---------------- src/IntaRNA/SeedHandlerMfe.h | 32 ++++++++++---- 3 files changed, 91 insertions(+), 47 deletions(-) diff --git a/ChangeLog b/ChangeLog index c037bc21..bddefa1d 100644 --- a/ChangeLog +++ b/ChangeLog @@ -10,10 +10,35 @@ # changes in development version since last release ################################################################################ +# 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/src/IntaRNA/SeedHandlerMfe.cpp b/src/IntaRNA/SeedHandlerMfe.cpp index 5c276ecd..a3bfaca8 100644 --- 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 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 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,7 +131,7 @@ 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 } @@ -135,7 +139,7 @@ fillSeed( const size_t i1min, const size_t i1max, const size_t i2min, const size 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,11 +293,11 @@ 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; } @@ -300,24 +305,24 @@ traceBackSeed( Interaction & interaction 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; diff --git a/src/IntaRNA/SeedHandlerMfe.h b/src/IntaRNA/SeedHandlerMfe.h index 0f770176..4938bd82 100644 --- a/src/IntaRNA/SeedHandlerMfe.h +++ b/src/IntaRNA/SeedHandlerMfe.h @@ -148,8 +148,8 @@ class SeedHandlerMfe : public SeedHandler /** * Provides the seed energy during recursion. * - * @param i1 the seed left end in seq 1 (index including offset) - * @param i2 the seed left end in seq 2 (index including offset) + * @param i1 the seed left end in seq 1 + * @param i2 the seed left end in seq 2 * @param bpInbetween the number of seed base pairs enclosed by the left- * and right-most base pair, ie. bpSeed-2 * @param u1 the number of unpaired bases within seq 1 @@ -168,8 +168,8 @@ class SeedHandlerMfe : public SeedHandler * you have to call the method in appropriate order depending on your seed * recursion. * - * @param i1 the seed left end in seq 1 (index including offset) - * @param i2 the seed left end in seq 2 (index including offset) + * @param i1 the seed left end in seq 1 + * @param i2 the seed left end in seq 2 * @param bpInbetween the number of seed base pairs enclosed by the left- * and right-most base pair, ie. bpSeed-2 * @param u1 the number of unpaired bases within seq 1 @@ -281,7 +281,7 @@ traceBackSeed( Interaction & interaction const size_t seedBps = getConstraint().getBasePairs(); // trace back the according seed - traceBackSeed( interaction, i1-offset1, i2-offset2 + traceBackSeed( interaction, i1, i2 , seedBps-2 , getSeedLength1(i1,i2)-seedBps , getSeedLength2(i1,i2)-seedBps ); @@ -334,9 +334,16 @@ E_type SeedHandlerMfe:: getSeedE( const size_t i1, const size_t i2, const size_t bpInbetween, const size_t u1, const size_t u2 ) const { +#if INTARNA_IN_DEBUG_MODE + if ( i1 < offset1 ) throw std::runtime_error("SeedHandlerMfe::getSeedE(i1="+toString(i1)+") is out of range (>"+toString(offset1)+")"); + if ( i1-offset1 >= seed.size1() ) throw std::runtime_error("SeedHandlerMfe::getSeedE(i1="+toString(i1)+") is out of range (<"+toString(seed.size1()+offset1)+")"); + if ( i2 < offset2 ) throw std::runtime_error("SeedHandlerMfe::getSeedE(i2="+toString(i2)+") is out of range (>"+toString(offset2)+")"); + if ( i2-offset2 >= seed.size2() ) throw std::runtime_error("SeedHandlerMfe::getSeedE(i2="+toString(i2)+") is out of range (<"+toString(seed.size2()+offset2)+")"); +#endif + return seedE_rec( SeedIndex({{ - (SeedRecMatrix::index) i1 - , (SeedRecMatrix::index) i2 + (SeedRecMatrix::index) (i1-offset1) + , (SeedRecMatrix::index) (i2-offset2) , (SeedRecMatrix::index) bpInbetween , (SeedRecMatrix::index) u1 , (SeedRecMatrix::index) u2 }}) ); @@ -349,9 +356,16 @@ void SeedHandlerMfe:: setSeedE( const size_t i1, const size_t i2, const size_t bpInbetween, const size_t u1, const size_t u2, const E_type E ) { +#if INTARNA_IN_DEBUG_MODE + if ( i1 < offset1 ) throw std::runtime_error("SeedHandlerMfe::setSeedE(i1="+toString(i1)+") is out of range (>"+toString(offset1)+")"); + if ( i1-offset1 >= seed.size1() ) throw std::runtime_error("SeedHandlerMfe::setSeedE(i1="+toString(i1)+") is out of range (<"+toString(seed.size1()+offset1)+")"); + if ( i2 < offset2 ) throw std::runtime_error("SeedHandlerMfe::setSeedE(i2="+toString(i2)+") is out of range (>"+toString(offset2)+")"); + if ( i2-offset2 >= seed.size2() ) throw std::runtime_error("SeedHandlerMfe::setSeedE(i2="+toString(i2)+") is out of range (<"+toString(seed.size2()+offset2)+")"); +#endif + seedE_rec( SeedIndex({{ - (SeedRecMatrix::index) i1 - , (SeedRecMatrix::index) i2 + (SeedRecMatrix::index) (i1-offset1) + , (SeedRecMatrix::index) (i2-offset2) , (SeedRecMatrix::index) bpInbetween , (SeedRecMatrix::index) u1 , (SeedRecMatrix::index) u2 }}) ) = E; From af97c327e3972f0a4a4172e1b6f52235de0b7e38 Mon Sep 17 00:00:00 2001 From: martin-mann Date: Fri, 31 Jan 2020 20:57:00 +0100 Subject: [PATCH 8/8] v3.1.5 --- ChangeLog | 8 +++++++- 1 file changed, 7 insertions(+), 1 deletion(-) diff --git a/ChangeLog b/ChangeLog index bddefa1d..7784a54d 100644 --- a/ChangeLog +++ b/ChangeLog @@ -10,11 +10,17 @@ # 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