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

MFE / pfunc macrostate computation wrong?! #26

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

MFE / pfunc macrostate computation wrong?! #26

sjanssen2 opened this issue Apr 21, 2021 · 3 comments
Labels
documentation Improvements or additions to documentation

Comments

@sjanssen2
Copy link
Member

sjanssen2 commented Apr 21, 2021

The macrostate MFE and pfunc computation is incorrect (when not using rna_turner1999 parameters).

Assume the following structure (lonely base pairs are allowed; last line gives names for structural components)

(...).(...).(...).(...)
CaaaGaCaaaGaCaaaGaCaaaG
aaaaaXbbbbbYcccccZddddd

The enum candidate for this structure in the macrostate grammar is (<- are manual annotations to relate sub-structures):

trafo(
  ambd_Pr(
    drem( <0,0> hl( <0,1> <1,4> <4,5> ) <5,5> )   <- (...) = aaaaa
    <5,6>                                         <- X
    ambd_Pr(
      drem( <6,6> hl( <6,7> <7,10> <10,11> ) <11,11> )   <- (...) = bbbbb
      <11,12>  						 <- Y
      ambd_Pr( 
        drem( <12,12> hl( <12,13> <13,16> <16,17> ) <17,17> )   <- (...) = ccccc
        <17,18> 					        <- Z
        cadd_Pr_Pr_Pr(
          drem( <18,18> hl( <18,19> <19,22> <22,23> ) <23,23> )    <- (...) = ddddd
	  nil( <23,23> ) 
        )
      )
    )
  )
)

The MFE algebra needs to internally choose from multiple options (instead of passing these options as individual candidates to the general choice function, as microstate does, which comes with an inflation of candidates in the search space rendering probability computation impossible): consider the leftmost two stems aaaaaXbbbbb which correspond the topmost ambd_Pr call. The unpaired base X can either dangle from right onto stem a or from left onto stem b or not dangle at all. However, if it dangles from left onto b we also need to consider if this is a more favorable situation than dangling X and Y onto b and forming an "external mismatch".
With Turner1999 parameters, the external mismatch was identical to left + right dangling, thus it could be computed stepwise. With the newer Turner2004 parameters, this is no longer the case and an external mismatch (ViennaPackage) can yield better energy than left and right dangle (ViennaPackage) together.
How can macrostate correctly decide here, without ending up with exponential many situations?

The problem recurses: while making the decision for X, we need to consider Y. However, the dangling of Z is dependent on Y :-/

@sjanssen2 sjanssen2 added the documentation Improvements or additions to documentation label Apr 21, 2021
@sjanssen2 sjanssen2 changed the title document macrostate issue MFE / pfunc macrostate computation wrong?! Apr 21, 2021
@sjanssen2
Copy link
Member Author

a more realistic, small example is the RNA sequence CCcCCaaaGGCCaaaGGuuGG which can form the structure ((.((...))((...))..)). Evaluation via RNAeval (version 2.4.17) and -d1 yields a free energy of 6.1, which is the same that gets reported by RNAshapes_eval_macrostate CCcCCaaaGGCCaaaGGuuGG "((.((...))((...))..))".

However, if we change the temperature setting from the default 37 to 57, RNAeval reports 7.86 kcal/mol, while RNAshapes_eval_macrostate gives 7.95 :-/

Here, the problem is within the algebra function mladldr and is caused by the asynchronity of mismatch_multi and dangling bases from left and right onto the same stem, i.e. there are multiple situations / algebra functions with the described issue.

@sjanssen2
Copy link
Member Author

Björn Voß confirms the above problem with his macrostate grammar (private communication).

@sjanssen2
Copy link
Member Author

referenced here: https://github.com/jlab/gapc/blob/80bc6668c59923572411d4333c98d836f5cf8862/librna/rnalib.c#L80 please don't break this link!

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
documentation Improvements or additions to documentation
Projects
None yet
Development

No branches or pull requests

1 participant