Dear Everyone,
I am going through the PyRosetta Tutorial kindly written by the Gray Lab. I am going through it one chapter at a time sequencially and i have reached the folding tutorial: http://graylab.jhu.edu/pyrosetta/downloads/documentation/pyrosetta4_online_format/PyRosetta4_Workshop4_Folding.pdf
My question:
I am stuck in the point where the tutorial asks to set up my own Monte Carlo method (the first step of the tutorial), I do understand Monte Carlo, but i don't understand how to set it up for this particular tutorial, I have looking through the D060_Folding.py script: http://graylab.jhu.edu/pyrosetta/downloads/scripts/demo/D060_Folding.py but i could not isolate the monte carlo script (is it built into PyRosetta (the tutorial does not discuss it)? Do i have to write a new one? i don't understand).
Can anyone help? how should i start the tutorial? is there a written script i can see how the monte carlo is setup? the way the tutorial asks?
Excerpt from the tutorial that confuses me:
Now, write a program that implements a Monte Carlo algorithm to optimize the protein conformation. In the main program, create a loop with 100 iterations. Each iteration should call a subroutine to make a random trial move, score the protein, and then accept or reject the new conformation based on the Metropolis criterion
Is this supposed to be something very easy as to not discuss in detail? where can i read about it in detial?
Thank you
AC
Rosetta does have a MonteCarlo class for the purpose, but the language you have pasted implies to me that implementing Monte Carlo is meant to be an exercise to the student. It seems to me to be a good choice: it's a mathematical algorithm that is pretty easy to describe in words, so translating that into code should be pretty easy.
You do understand Monte Carlo - so let me make this suggestion. First, think of how it works in human terms. Then, write that down in pseudocode. Next, think about how the pseudocode units translate into Rosetta concepts - the "random trial move" is a mover apply() function; the "score the protein" is using the scorefunction, etc. Then, translate your pseudocode into proper Python and try it out. I guess I can give more detail here if you need it?
I don't think this is supposed to be "very easy" per se, but I think it's supposed to be composed of parts you'll already understand. (I'm guessing - I've never looked at the tutorial). So I assume the challenge is meant to be stringing the parts together into a coherent whole.
Thank you very much for your replay.
So i took your advice and been persistantly working trying to write the script myself. The only description available is in the first page (not very detailed). Link to this particular tutorial: http://graylab.jhu.edu/pyrosetta/downloads/documentation/pyrosetta4_online_format/PyRosetta4_Workshop4_Folding.pdf
I have attached my script (it does not have the PyMol mover since it is disracting me from the actual math).
I started with a loop (as instructed) with 100 itirations (monte carlo sample size, correct?).
inside it i attempt to write the three subroutines in the tutorial (1 - Random Trial Move: randomly choose a residue, then choose either phi or psi to change. 2 - score the structure before and after the angle change. 3 - Metropolis Criterion to choose the best structure).
I think (but i am not confident) that i got everything ok except last step (Metropolis Criterion). The Random Trial Move updates the pose perminantly correct? or is there an additional step to revert back to the original starting structure?
The last step (Metropolis Criterion) is my biggest issue right now. i cannot continue the script without knowing how to use the Metropolis Criterion correctly to choose the best structure. To be honest, i do not understand the Metropolis Criterion itself so i am having difficulty implementing it.
Can you be so kind and check my work and see if i am on the right path?
AC
Reverting the pose: You are correct that the pose is updated correctly. You've already got the idea of a "starting energy" - why not a "starting pose"? Then reverting the move if it rejected is as simple as copying the old pose over the rejected new one.
Metropolis monte carlo: Monte Carlo is the random bit. The Metropolis criterion is the mathy bit. I'm not sure you've written it correctly (it looks like you've written it so that negE is always negative, which can't be right) but it's close. The point of the Metropolis decision is this: GOOD changes (new energy < old energy) are always accepted. BAD changes (new > old) are accepted SOMETIMES, with that SOMETIMES being less frequent as the energy increase gets larger. This means we always go "downhill" in the energy landscape when possible, and uphill sometimes, to escape local minima.
You're clearly on the right path.
Ok, I think I got it.
This Metropolis Criterion broke my brain, but i think i nailed it.
So i divided the script into two sections (one for phi angel and the other for psi angel), that way i can revert back to (pre angle change state) if the structure is rejected (this is to over come the issue that the pose is updated globally), or continue to next itiration if structure is accepted.
i wrote the Decision(probability) function which takes a number between 0 and 1 (i.e: P) where closer to 1 give randomly more True and closer to 0 gives randomly more False and use that to make the Accept/Reject decision.
If everything is ok so far i think i am only left with the Metropolis Criterion's function: I am not sure i am using the E < 0 equation correctly, i tried to write what was written in the tutorial: it seems the equation is written where the change in E is always negative, which is why i make the variable negE always negative (answer to your previous question). Is this correct? or am i miss reading the equation in the tutorial?
also kt = 1 correct? i am a bit confused because the tutorial says 1 Rosetta Energy Unit (is this a different quanta? or just simply 1?)
AC
Why is the energy always negative? It's because we've already passed through the "is the change good or bad" check. If the change is good, we just always accept it. If the change is bad, we use the Metropolis criterion. So we only get to the Metropolis critierion when the energy is increasing. If pay attention to the math - the probability will ALWAYS be on (0,1) for a change where the energy increases. If the energy DECREASED, the sign flips, and you end up with a "probability" above 1, on (1, infinity]. I think that is the source of your confusion about how the equation works. I still think you should just multiply by negative one - NOT take the absolute value and THEN multiply by negative one.
kt is arbitrary, because REU are arbitrary. 1 is reasonable. I've been told that 0.8 matches physiologic temperatures for REU calculations, but it's a tenuous connection. I think with REF2015 there is a definite value for kt, but I don't remember what it is.
If you want to peek at code, let me suggest finding the MonteCarlo class in your Rosetta install (it's probably MonteCarlo.cc somewhere in the protocols library). It implements this function in C++.
Your code would have a lot less duplication if you reset your poses via pose copy instead of resetting phi/psi manually. Think of extending this to other types of changes - clearly you can't write an undoer for any possible action, but copying a stored pose over the current pose can undo any change.
Thank you very much for all your explanations, I finished the tutorial and i think i got the script right (attached).
As you advised i found a way to revert back to previous pose: (if move is accepted update main pose - if rejected revert back to previous one).
I also removed the absolute statement from the Metropolis Criterion, but left kt=1 just to follow the tutorial.
finally i reached the point of the tutorial where i insert the 3-mer fragment.
I have been using the C++ Abinitio binary for a year now and this tutorial (and your explanations) have actually explained A LOT of gaps in my knowlege, so thank you very much.
Final question: Does everyone write their own PyRosetta Abinitio script? or is there a standarnd class/function/method that preforms an Abinitio fold that comes with PyRosetta? or does everyone just use the C++ Abinitio binary by calling it from a PyRosetta script?
PyRosetta exposes most of C++ rosetta. Rosetta is structured as "libraries". The C++ abinitio executable calls C++ code in those libraries; PyRosetta can also call C++ code in those libraries. So I imagine most users that want access to standard ab initio via Python will do so by calling the ab initio protocol/mover (in code), but not the executable (with a shell or exectuable call). (Not that I know what it's named - in a sane world it would be AbInitioMover, or similar, but that's not the world we live in.)