RNA-FRETools / MASH-FRET

MATLAB package to analyze single-molecule FRET data
https://rna-fretools.github.io/MASH-FRET/
GNU General Public License v3.0
8 stars 2 forks source link

simulation pre-set file not read properly #104

Closed LeonieUniFR closed 1 year ago

LeonieUniFR commented 1 year ago

Hello, I encountered some problems when running simulations with MASH-FRET. I created presetFiles.mat with createSimPrm.m.

  1. When I entered values for trans_prob and left trans_rates empty nothing was imported when using 'import presets'. Is it not sufficient to give only the probabilities?
  2. When I entered values for trans_rates and left trans_prob empty rates were imported when using 'import presets' (numbers appeared in the Thermodynamic model and turned grey). BUT the resulting traces don’t fit my numbers). When I add the same rates by hand into the table on the GUI completely different traces (ones which fit my rates) are generated. I set the frame rate to 1 to prevent some frames^-1 to seconds^-1 miscalculations.
  3. I created a transition matrix with 4 states where the molecule can only go in one direction. The state sequence should always be 13241324… But when I generate traces the state sequence in the exported ASCII does not match. Are some states not in the traces because their dwell time is too short and they are missed?#

Thank you for your help, Leonie Email: leonie.vollmar@pc.uni-freiburg.de

Rates in GUI -> traces with long dwell times ratesinGUI Rates imported via presets -> very short dwell times ratesinpresets

mca-sh commented 1 year ago

Hello Leonie and thank you for reporting this issue.

I have to admit that the way presets are imported in MASH is not very clear. And thanks to your report, I see that it is also not correct and I must fix it.

  1. trans_prob is not supposed to contain the transition probabilities but the transition partition probabilities/factors. They are transition probabilities but normalized by the row sum and with 0 diagonal probabilities. It must be imported together with state lifetimes that are stored in field lifetimes of the preset file. This is one incorrect part: it was not defined in the createSimPrm.mat. But apparently, you have transition probabilities ready in your hands, so you can use the field trans_rates only. Keep in mind that value in trans_rates have no units as they stand for the number of transitions per time interval (and not per second).
  2. This is absolutely true, apparently MASH is not handling properly preset rate matrix with non-null diagonal terms. this is now corrected in the new version (see below).
  3. Yes indeed, when multiple transition occur within one time interval, the exported state sequence retains the state that populate majorly this time interval. To verify if it is the case for you, your can export the dwell time files and verify that the state sequence in them is correct and that the dwell times are shorter than 1s in your case. If it still does not correspond to the parameters you've defined, I will have to look deeper into it.

You can download the corrected version: https://github.com/RNA-FRETools/MASH-FRET/archive/refs/heads/review-sim-presets.zip, or fork the branch review-sim-presets. Please tell me if it is working for you so I close the issue and merge the corrections in the main branch.

Mélodie

LeonieUniFR commented 1 year ago

Hello Mélodie, Thank you for you quick reply.

  1. Okay, I will only use trans_rates from now on. Also thank you for pointing out that the units of trans_rates is per time interval and not per second. This led to some confusion on my side.
  2. I tested your new version. Unfortunately it crashed and gave the following error in buildModel (see screenshot). I then stopped the program right bevor the error. When I understand it correctly the size of ip seems to be the problem. Shouldn't it be ip =ip_all(n,:)' in line 86 in buildModel instead of ip = ip_all(:,n)'? error_buildModel
  3. Thanks. It totally makes sense when looking at the exported dwell times.

Best, Leonie

mca-sh commented 1 year ago

Ouch.. I made a mess here Almost, actually it should be ip_all(n,:), Thank you for your feedback, modifications have been made to the same branch, You can use the same download link I gave you previously. Mélodie

LeonieUniFR commented 1 year ago

Hi Mélodie, unfortunately the dimentions still produced an error in buildModel and a new error turned up:

Matrix dimensions must agree. Error in buildModel (line 72) ip_all = tau_all./repmat(sum(tau_all,1),[J,1,1])'; % [J-by-N]

I then had to make a change in line 72 ip_all = tau_all./repmat(sum(tau_all,1),[J,1,1]); % removed ' and add an additional line 73 ip_all = ip_all'; % now transposed

now ip_allhas the same dimentions when being calculated or when coming from presets. Now it works :)

Best, Leonie

mca-sh commented 1 year ago

Oh my.. Thank you Leonie for the feedback AND the bug fix. I will add this correction to the source code and close this issue then. Do not hesitate to post a new issue anytime you encounter such malfunctioning, this is very helpful for us. Best, Mélodie