Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

window mode +/-1 offset? #24

Open
sjanssen2 opened this issue Apr 21, 2021 · 8 comments
Open

window mode +/-1 offset? #24

sjanssen2 opened this issue Apr 21, 2021 · 8 comments
Labels
bug Something isn't working

Comments

@sjanssen2
Copy link
Member

foldgrammars: eval gives correct overdangle free enery of -13.30 kcal/mol:

(gapc) t490s x86_64 /Daten/Git/jlab/fold-grammars/Misc/Applications/RNAshapes>./x86_64-linux-gnu/RNAshapes_eval_overdangle -P /home/sjanssen/share/gapc/librna/rna_turner2004.par AAGGGCGUCGUCGCCCCGAGUCGUAGCAGUUGACUACUGU "..(((((....))))).........(((((.....)))))"
Answer: 
( ( ..(((((....))))).........(((((.....))))) , -1330 ) , [][] )

using a longer sequence (with the very same first 40 nucleotides) in window mode returns wrong free energy: -13.70

(gapc) t490s x86_64 /Daten/Git/jlab/fold-grammars/Misc/Applications/RNAshapes>./x86_64-linux-gnu/RNAshapes_mfe_overdangle_window -w 40 -i 20 -P /home/sjanssen/share/gapc/librna/rna_turner2004.par AAGGGCGUCGUCGCCCCGAGUCGUAGCAGUUGACUACUGUUuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuu
Answer (0, 40) :
( -1370 , ( ( ..(((((....))))).........(((((.....))))) , [][] ) , 28778.4 ) )
Answer (20, 60) :
( -470 , ( ( ....((((((.....))))))................... , [] ) , 0.0130966 ) )
Answer (40, 80) :
( 0 , ( ( ........................................ , _ ) , 6.38663e-06 ) )
Answer (60, 89) :
( 0 , ( ( ............................. , _ ) , 0.000171326 ) )

also, RNAeval from Vienna package (2.4.17) confirms energy of 13.30:

(base) t490s x86_64 /Daten/Git/jlab/fold-grammars/Misc/Test-Suite/StefanStyle/x86_64-linux-gnu>~/src/ViennaRNA-2.4.17/src/bin/RNAeval -v -d2 -T 37 -P /home/sjanssen//share/gapc/librna/rna_turner2004.par

Use '&' to connect 2 sequences that shall form a complex.
Input sequence (upper or lower case) followed by structure; @ to quit
....,....1....,....2....,....3....,....4....,....5....,....6....,....7....,....8
AAGGGCGUCGUCGCCCCGAGUCGUAGCAGUUGACUACUGU
..(((((....))))).........(((((.....)))))

length = 40
External loop                           :   -80
Interior loop (  3, 16) GC; (  4, 15) GC:  -330
Interior loop (  4, 15) GC; (  5, 14) GC:  -330
Interior loop (  5, 14) GC; (  6, 13) CG:  -340
Interior loop (  6, 13) CG; (  7, 12) GC:  -240
Hairpin  loop (  7, 12) GC              :   400
Interior loop ( 26, 40) GU; ( 27, 39) CG:  -250
Interior loop ( 27, 39) CG; ( 28, 38) AU:  -210
Interior loop ( 28, 38) AU; ( 29, 37) GC:  -210
Interior loop ( 29, 37) GC; ( 30, 36) UA:  -220
Hairpin  loop ( 30, 36) UA              :   480
AAGGGCGUCGUCGCCCCGAGUCGUAGCAGUUGACUACUGU
..(((((....))))).........(((((.....)))))
 energy = -13.30 kcal/mol

Is the window mode broken in gapc?

@sjanssen2 sjanssen2 added the bug Something isn't working label Apr 21, 2021
@sjanssen2
Copy link
Member Author

I think the issue is with dall which calls ext_mismatch_energy from rnalib.c. The later probably runs out of the window sub-sequence for P->mismatchExt[closingBP][lbase][rbase].

@sjanssen2
Copy link
Member Author

I am addressing this issue via PR in gapc repo

@sjanssen2
Copy link
Member Author

last window isn't printed by gapc binary :-/

@sjanssen2
Copy link
Member Author

sjanssen2 commented Apr 22, 2021

looks like it makes a difference if product is compiled with --kbacktrace or not!
Issue transfers from pAliKiss to RNAalishapes.
we loose the best solution for window (4, 19)

RNAalishapes>./x86_64-linux-gnu/RNAalishapes_subopt_overdangle_window -i 1  -P /home/sjanssen/share/gapc/librna/rna_turner1999.par -u 1        "G_____CGGGAGAC____CG#G_____CGGGAGACG___CG#G_____CGGGAGACG___CG#" -C 0 -w 15 -e 0.1
Answer (0, 15) :
( ( -93.3333 = energy: -93.3333 + covar.: 0 ) , ( ( ........(....). , [] ) , 0.0512456 ) )
Answer (1, 16) :
( ( -93.3333 = energy: -93.3333 + covar.: 0 ) , ( ( .......(....).. , [] ) , 0.0512456 ) )
Answer (2, 17) :
( ( -93.3333 = energy: -93.3333 + covar.: 0 ) , ( ( ......(....)... , [] ) , 0.0512456 ) )
Answer (3, 18) :
( ( -93.3333 = energy: -93.3333 + covar.: 0 ) , ( ( .....(....).... , [] ) , 0.0512456 ) )
Answer (4, 19) :
( ( -86.6667 = energy: -86.6667 + covar.: 0 ) , ( ( ...((....)....) , [] ) , 0.0459917 ) )
Answer (5, 20) :
( ( -296.667 = energy: -296.667 + covar.: 0 ) , ( ( .(((....)....)) , [] ) , 1.38822 ) )
(base) t490s x86_64 /Daten/Git/jlab/fold-grammars/Misc/Applications/RNAalishapes>./x86_64-linux-gnu/RNAalishapes_subopt_overdangle_window -i 1  -P /home/sjanssen/share/gapc/librna/rna_turner1999.par -u 1        `cat ~/eclipse-workspace/librna_change/Debug/mininput.msa | tr "\n" "#"` -C 0 -w 15 -e 0
Answer (0, 15) :
( ( -93.3333 = energy: -93.3333 + covar.: 0 ) , ( ( ........(....). , [] ) , 0.0512456 ) )
Answer (1, 16) :
( ( -93.3333 = energy: -93.3333 + covar.: 0 ) , ( ( .......(....).. , [] ) , 0.0512456 ) )
Answer (2, 17) :
( ( -93.3333 = energy: -93.3333 + covar.: 0 ) , ( ( ......(....)... , [] ) , 0.0512456 ) )
Answer (3, 18) :
( ( -93.3333 = energy: -93.3333 + covar.: 0 ) , ( ( .....(....).... , [] ) , 0.0512456 ) )
Answer (4, 19) :
Answer (5, 20) :
( ( -296.667 = energy: -296.667 + covar.: 0 ) , ( ( .(((....)....)) , [] ) , 1.38822 ) )

If compiled without --kbacktrack we find the best solution:

Answer (4, 19) :
( ( -93.3333 = energy: -93.3333 + covar.: 0 ) , ( ( ....(....)..... , [] ) , 0.0512456 ) )

@sjanssen2
Copy link
Member Author

if we aim to look at the full search space by increasing -e 100000, we see that every window lacks candidates that are present a) if compiled without --kbacktrace or if we print the forward res in out_main.cc.

@sjanssen2
Copy link
Member Author

A non alignment, gapc self contained example:

(base) t490s x86_64 /Daten/Git/jlab/gapc_libtest/testdata/regresstest/temp>gapc -p "mfe * pretty" ../../grammar/adpf.gap    --window-mode -t --kbacktrack 
Warning: 
  scoring choice [double] h([double] l)
                          ^
../../grammar/adpf.gap:590.27: Declared role of algebra choice function h (in algebra p_func_sample) does not match autodetected role none(n).
(base) t490s x86_64 /Daten/Git/jlab/gapc_libtest/testdata/regresstest/temp>make -f out.mf CXXFLAGS_EXTRA="-I ../../gapc_filter/"
Including global makefile /home/sjanssen/share/gapc/config_linux-gnu.mf
echo '#include "out.hh"' > out_main.cc
cat /home/sjanssen/include/rtlib/generic_main.cc >> out_main.cc
g++  -I/usr/include  -I/usr/include -I/home/sjanssen/include/ -I/home/sjanssen/include/rtlib -I/home/sjanssen/include/librna  -O3  -std=c++17 -D_XOPEN_SOURCE=500 -MMD -MP -Wall -Wnon-virtual-dtor -Wno-unused-variable -Wno-parentheses -DNDEBUG -DBISONNEW  -I ../../gapc_filter/ out_main.cc -c -o out_main.o 
g++  -I/usr/include  -I/usr/include -I/home/sjanssen/include/ -I/home/sjanssen/include/rtlib -I/home/sjanssen/include/librna  -O3  -std=c++17 -D_XOPEN_SOURCE=500 -MMD -MP -Wall -Wnon-virtual-dtor -Wno-unused-variable -Wno-parentheses -DNDEBUG -DBISONNEW  -I ../../gapc_filter/ out.cc -c -o out.o 
g++  -o out out_main.o out.o string.o  -L/usr/lib/x86_64-linux-gnu  -L/home/sjanssen/lib/ -L/home/sjanssen/lib/rtlib -L/home/sjanssen/lib/librna   -Xlinker -rpath -Xlinker /usr/lib/x86_64-linux-gnu      -Xlinker -rpath -Xlinker /home/sjanssen/lib/ -Xlinker -rpath -Xlinker /home/sjanssen/lib/rtlib -Xlinker -rpath -Xlinker /home/sjanssen/lib/librna    -L/usr/lib/x86_64-linux-gnu -lgsl -lgslcblas -lm -lrna 
(base) t490s x86_64 /Daten/Git/jlab/gapc_libtest/testdata/regresstest/temp>./out ucgaagacuugguccuuuggug -P ../../../librna/vienna/rna_turner1999.par -w 20 -i 1 
Answer (0, 20) :
( -300 , .(((((........))))). )
Answer (1, 21) :
Answer (2, 22) :

@sjanssen2
Copy link
Member Author

iff ext_mismatch_energy always returns the same value, the problem vanishes?!

@sjanssen2
Copy link
Member Author

here comes a really small example containing this error:

Save the following code as adpf.gap, compile via

(base) t490s x86_64 /Daten/Git/jlab/gapc_libtest/testdata/regresstest/temp>gapc -p "mfe * (enum)" ../../grammar/adpf.gap --window-mode -t --kbacktrack && make -f out.mf CXXFLAGS_EXTRA="-I ../../gapc_filter/"

and execute binary via

(base) t490s x86_64 /Daten/Git/jlab/gapc_libtest/testdata/regresstest/temp>./out ucgaag -w 4 -i 1 which will result in the output

Answer (0, 4) :
( -160 , cadd( dlr( <0,0> hl( <0,2> ) <2,2> ) cadd( dlr( <2,2> hl( <2,4> ) <4,4> ) nil( 1 ) ) ) )
Answer (1, 5) :
( -160 , cadd( dlr( <1,1> hl( <1,3> ) <3,3> ) cadd( dlr( <3,3> hl( <3,5> ) <5,5> ) nil( 1 ) ) ) )
Answer (2, 6) :

adpf.gap:

import rna
import adpf_filter

input rna

signature FS (alphabet, comp) {
  comp sadd(Subsequence, comp);
  comp cadd(comp, comp);
  comp dlr(Subsequence, comp, Subsequence);
  comp hl(Subsequence);
  comp nil(void);
  choice [comp] h([comp]);
}

algebra mfe implements FS(alphabet = char, comp = int) {
  int sadd(Subsequence lb, int e) {
    return e;
  }
  int cadd(int x, int y) {
    return x + y;
  }
  int dlr(Subsequence lb, int x, Subsequence rb) {
    return x;
  }
  int hl(Subsequence lb) {
    return fake_energy(lb, lb);
  }
  int nil(void) {
    return 0;
  }
  choice [int] h([int] i) {
    return list(minimum(i));
  }
}

algebra enum auto enum;
algebra count auto count;

grammar fold uses FS(axiom = struct) {
   struct = sadd(REGION with minsize(2) with maxsize(2), struct)
          | cadd(dangle, struct)
          | nil(EMPTY) 
	  # h ;

   dangle = dlr(LOC, hl( REGION with minsize(2) with maxsize(2)), LOC) 
	  # h;
}

instance ppmfe = fold (  mfe ) ;

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
bug Something isn't working
Projects
None yet
Development

No branches or pull requests

1 participant