This question is regarding the implementation of the Dunbrack rotamer energy term (fa_dun).
Based on my understanding of the implementation in Rosetta 2.3, I noticed that the fa_dun energy is comprised of two components:
1) A backbone-dependent probability of the rotamer : E = -log (p)
2) An "un-normalized" gaussian penalty parametrized using the average chi angles and the standard deviations for the rotamer: E_penalty = -log (gaussian)
My concern is regarding the use of an "un-normalized" gaussian penalty term. My understanding of the Dunbrack rotamer libraries is that the probabilities are for the "entire" rotamer(i.e. over the entire distribution of the chi-angles for a rotamer within the cutoffs that classify an angle as g+, g- or t), not just the rotamer conformation *at* the mean values listed in the rotamer library file. Therefore, I believe that there is a correction term that is missing for the penalty terms that would be comprised of the normalization constants of the gaussian:
E_correction = -log (normalization of the gaussian).
As long as one stays within the same rotamer, the correction term would cancel out between the different side-chain perturbations. However, as one samples different rotamers or samples different amino-acids(for a design calculation), the correction terms would be different since the chi-angle standard deviations are different.
Please let me know if this is a known issue or if I am misunderstanding either the rotamer libraries or their implementation in Rosetta.
-Tushar.
The guy who knows best (Andrew) answers:
"
I've run across this a number of times; the behavior described here is present both in Rosetta2 and Rosetta3.
Yes: the normalization of the gaussian is omitted. I don't know the history behind this, but have experimented several times by adding the gaussian back in. Strangely, it hurts rotamer recovery. For this reason, I haven't tried to persuade the community to fix this behavior.
But you're right: it's not a true log-likelihood score since the integral of the probability over all chi angles doesn't come to 1."
Thanks for the reply. I would just like some clarifications on the response.
1) When you say that it hurts rotamer recovery, is that after the energy coefficients have been refitted with the normalization? I ask this because changing the term would require a new set of coefficients. If you do have a new set of coefficients, would you mind sharing the wts_patch file for us to try in some of the problems that we are looking at?
2) Do rotamer recoveries decline across the board for most amino-acids? Or is it more of a mixed-bag?
3) How about the influence on design problems where I think that this will have the biggest impact?
1) Andrew is referring to rotamer recovery used as a test while fitting those coefficients, so, yes. We don't have a weights file for it around anymore (and other code will have shifted, so it will have bitrotted anyway).
We don't have answers to 2 and 3 - if you find interesting things let us know!
0#said:
This question is regarding the implementation of the Dunbrack rotamer energy term (fa_dun).
Based on my understanding of the implementation in Rosetta 2.3, I noticed that the fa_dun energy is comprised of two components:
1) A backbone-dependent probability of the rotamer : E = -log (p)
2) An "un-normalized" gaussian penalty parametrized using the average chi angles and the standard deviations for the rotamer: E_penalty = -log (gaussian)
My problem:
My concern is the use of an "un-normalized" gaussian penalty term.
1、Calculating fa_dun int the default case , whether to use penalty term. if without penalty term, how can I add it to calculat the fa_dun score;
2、I want to know the calculation formula of fa_dun with penalty.In the code of rosettta, I didn't find out about this calculation, whether you can give some hints for my learning(I want to know something about Scoring functions ).
Thank you very much.
I can't tell you much about how the Dunbrack term was implemented in 2.3. It wasn't much changed for the 3.x series, at least.
3.x has supported 3 different Dunbrack libraries through its history. The 02 library (default) is presumably the one you are asking about. The 08 library was found to be "over-smoothed" and is deprecated in favor of the 10 library.
The code that describes how the Dunbrack Energy is calculated is, as of 3.4, in core/pack/dunbrack/RotamericSingleResidueDunbrackLibrary.tmpl.hh. The function is named eval_rotameric_energy_deriv, and the part you care about starts at line 569 or so.
The penalty function cannot be turned on or off without changing the code (there's no command line option for it) - it's just always on. The code is this:
for ( Size ii = 1; ii <= T; ++ii ) {
/// chidev penalty: define a gaussian with a height of 1 and a standard deviation of chisd[ ii ];
/// exp( -1 * (chi_i - chi_mean_i )**2 / ( 2 * chisd_ii**2 ) )
/// Pretend this gaussian defines a probability for having an angle at a certain deviation from the mean.
/// Convert that probability into an energy:
/// -ln p = -ln exp( -1* (chi_i - chi_mean_i)**2/(2 chisd_i**2) ) = (chi_i - chi_mean_i)**2/(2 chisd_i**2)
chidevpen[ ii ] = chidev[ ii ]*chidev[ ii ] / ( 2 * chisd[ ii ] * chisd[ ii ] );
}
Real chidevpensum( 0.0 );
for ( Size ii = 1; ii <= T; ++ii ) {
chidevpensum += chidevpen[ ii ];
}
scratch.fa_dun_tot() = -std::log( scratch.rotprob() ) + chidevpensum;
scratch.fa_dun_rot() = -std::log( scratch.rotprob() );
scratch.fa_dun_dev() = chidevpensum;
Real const score( -std::log( scratch.rotprob() ) + chidevpensum );
if ( ! eval_deriv ) return score;
Thank you for your patient guidance, these are very useful, very gratefull!