I have used the minirosetta application with the membrane_highres_t1.wts scoring function to create 1817 models. Among them only 48 had negative score, which ranged between -948.737 and -633.723. In contrast the rest of the models had score between 37700.552 and 197471.510. How can this happen? Does it mean that the scoring function is not suitable for my transmembrane protein?
Below is my parameter file:
-in:file:fasta FurA_20-500.fasta # target sequence
-in:file:template_pdb PDB_FILES/2X79_A.pdb # template structure
-in:file:fullatom # input will be fullatom
-in:file:psipred_ss2 SAM_files/FurA_20-500_renumbered.dssp-ehl2.ss2 # optional: SAM predicted ss converted to 'psipred ss2' format
-in:file:alignment FurA_20-500-truncated_Mhp1_2X79_A_publication_alignment.aln # input alignment in Clustal format
## Loop specific options
-loops:frag_sizes 9 3 1
-loops:frag_files fragments_SAM_rdb/FurA_.200.9mers fragments_SAM_rdb/FurA_.200.3mers none
-loops:remodel quick_ccd # use ccd to remodel loops
-loops:idealize_after_loop_close # idealize structure after closing loops (structure will have ideal rosetta bond lengths and angles)
-loops:extended true # force extended on loops (phi-psi angles set to 180 degrees) independent of loop input file, for rebuilding loops entirely.
-loops:build_initial true # precede loop modeling with initial round of removing missing densities and building simple loops
-loops:refine refine_ccd # small movements to remodel loop
-silent_decoytime # Add time since the last structure was written to score line
-select_best_loop_from 100 # Keep building loops until N and choose best (by score)
-loops:relax fastrelax # fastrelax loops (5 cycles at default)
-cm:min_loop_size 4 # minimum size of loops to consider for rebuild
## Membrane specific input information
## Output format
-out:file:fullatom # output file will be fullatom
-out:file:scorefile FurA_20-500-2X79_A_threading_models.fasc # output a table of Rosetta scores
-nstruct 1000 # number of structures generated
-run:constant_seed # Use a constant seed (1111111 unless specified)
-run:jran 39596021 # this is good for testing since you should always get the same result
-overwrite # overwrite any already-existing results having the same name
Scores that large are almost certainly due to clashes (atom overlaps). Most of the large score should be in the fa_rep term. Look at the score function summaries and tell us which term it's in. What do the awful ones look like?
Presumably Rosetta is making poor early-stage (centroid) models that have an irresolvable clash in fullatom mode. (You are using quick_ccd instead of a longer protocol, try perturb_ccd). The fullatom scorefunction is very rugged and may not be able to find its way out of a bad clash.
If you scores are really distributed with some negative and the rest hugely positive, but nothing in the middle (no weak positives)...that probably indicates something is wrong but I'm not quite sure what.
Here are the scores of the top scored model:
SCORE: score fa_atr fa_rep fa_intra_rep fa_mbenv fa_mbsolv pro_close fa_pair hbond_sr_bb hbond_lr_bb hbond_bb_sc hbond_sc dslf_ss_dst dslf_cs_ang dslf_ss_dih dslf_ca_dih rama omega fa_dun p_aa_pp ref aln_ident aln_len aln_perc brlx_irms cen_irms final_looprelax_score irms nres silent_score time aln_id template description
SCORE: -948.737 -1671.695 284.025 4.667 0.130 585.420 4.367 -13.844 -355.482 -15.681 -30.315 -28.284 0.000 0.000 0.000 0.000 22.159 49.422 228.958 -5.984 -6.600 13.987 0.222 12.843 13.236 409.000 12.843 2.782 449.000 -948.736 -948.736 1.098 2JLO_A.pdb 2JLO_A empty_19
And the best positively scored model:
SCORE: 37700.552 -1727.293 9243.872 36.113 310.972 415.777 28433.781 -9.460 -247.533 -3.073 -9.287 -6.299 0.000 0.000 0.000 0.000 216.842 40.914 975.099 36.727 -6.600 0.222 449.000 1.098 12.057 12.057 37700.555 12.057 409.000 37700.555 6174.000 2JLO_A 2JLO_A.pdb S_2JLO_A_0545_1
I just had a quick look at these two models as well as 2 other negatively scored, but couldn't spot any suspicious difference. Half of the differences are in the transmembrane part of the protein. Do you think these transmembrane loops are disfavored by the scoring function and hence increase the score that much?
Interesting! It's not all in fa_rep, it's mostly in pro_close. pro_close is a term that ensures proline's ring stays closed. Rosetta uses a directed acyclic graph (http://en.wikipedia.org/wiki/Directed_acyclic_graph) in the AtomTree - the lack of cycles means that ring systems have a pseudobond. pro_close is the term that ensures proline residues chi angles leave the "predicted" virtual nitrogen NV atom atom the real N atom - ensuring the pseudobond between the last carbon and nitrogen stays in place.
I've never seen pro_close this horrible...can you attach to the forum or send me one of the bad PDBs? (There's a way to email me directly over on the left). If not that, can you use the per-residue scores section to ID a bad proline and clip out it and a residue to either side?
It may be that you're getting a cis-trans isomerization at a proline (should be ok...), or something else is fluky near a proline.
I have emailed you the 5 most negatively scored and the 5 less positively scored models, along with the score.sc file that score_jd2 produces. What flag prints the per-residue scores? In another thread you said it is -out:pdb but I don't get the scores in the pdbs. These are the flags I used:
-in:file:tags empty_19 empty_40 empty_16 empty_37 empty_13 empty_34 empty_12 empty_32 empty_46 empty_43 empty_29 empty_8 empty_47 empty_45 empty_11 empty_33 empty_10 empty_31 empty_14 empty_35 empty_42 empty_25 empty_4 empty_24 empty_3 empty_1 empty_22 empty_30 empty_9 empty_20 empty_41 empty_15 empty_36 empty_27 empty_6 empty_2 empty_23 empty_26 empty_5 empty_18 empty_39 empty_17 empty_38 empty_28 empty_7 empty_44 empty empty_21
Actually, now that I had a closer look at the models I've sent you, proline 271 is broken in all the positively scored models but intact in the negatively scored. Also the bad models have unclosed loops (203-204, 302-306), which means the alignment can be improved in that region.
They're in the PDBs you sent me. Open it up in a text editor (or Excel, using "text import", space-delimited, merge delimiters), and you can see all the per-residues scores in a table at the bottom. We print them space-delimited instead of table formatted so you need a spreadsheet to get them aligned visibly with the headers.
The PDB problems:
Wow, this isn't just prolines - in the "failed" models (the ones not called empty), there are a ton of broken sidechains - in particular, the CB-CG (or OG) bond is being shattered. For example, in model #535, look at W61 - the CB-CG "bond" is 20 angstroms!
The scorefunction is only picking this up through the occasional clash (fa_rep) and proline closure (which is destroyed by ludicrous bond lengths) - normally, Rosetta uses only ideal bond lengths, and so this sort of problem isn't "scored for" because it's not supposed to ever occur.
Unfortunately I have no clue what is causing this - it's a very unlikely thing for Rosetta to do, because that bond length isn't a degree of freedom.
How long are your loops? I wonder if there might be a bug in -select_best_decoy_from that is inappropriately caching sidechains or something.
If you want to repair the "bad" models, you can easily do so with:
fixbb -database $DATABASE -ex1 -ex2 -packing:repack_only -ndruns 10 -nstruct 1 -s *pdb
depending on where the models went wrong, this might repair them, or they might still be awful. It will take 30 s/structure at most.
I tried to fix one of the bad models, the sidechains were repaired but placed over other neighboring atoms. How can I fix that without relaxing the structure? Also the broken loops are not repaired.
I have also recompiled minirosetta and started a job. I will know tomorrow if the terminal loop is remodeled.
"I tried to fix one of the bad models, the sidechains were repaired but placed over other neighboring atoms. How can I fix that without relaxing the structure?"
You can't - I was hoping there'd be space to put the sidechains back in, but apparently there isn't. The mainchains are probably too closely packed.
I have tried the -loops:remodel perturb_ccd option instead of -loops:remodel quick_ccd to see if that improves the models, but Rosetta complains about the terminal loop:
ERROR: !loop.is_terminal( pose ) || pose.fold_tree().is_cutpoint( loop_cut )
ERROR:: Exit from: src/protocols/loops/loop_mover/perturb/LoopMover_CCD.cc line: 212
Is there any way to circumvent that problem?
The assert statement on that line is wrong. Remove the exclamation point (!), recompile, and see if it works. I filed a bug report and am testing the fix in the developer codebase, but I don't have a terminal loop handy to test on. https://carbon.structbio.vanderbilt.edu/mantisbt/view.php?id=230
I tested that removing the exclamation point should fix the runtime_assert failure. Unfortunately in my test case Rosetta does not actually modify the terminal loop any - no crash, but no modeling - let me know how it turns out for you.
OK, I have recompiled minirosetta and created 75 models with "-loops:remodel perturb_ccd". 74 models had scores ranging between -888.885 and -212.126, and just one had score 25.175 but was completely unfolded. Are these scores normal for a transmembrane protein of 409 amino acids?
Also, no loop was broken and all the stidechains were intact. However, as you can see in the models I have emailed you, the long loops (50-70, 131-158, 202-229) are unstructured compared with the ones created with "-loops:remodel quick_ccd" option. Is this normal?
"Are these scores normal for a transmembrane protein of 409 amino acids?"
The rule of thumb for score12 is that -2 per residue is okay. Your scorefunction appears closely related to score12, so -888 seems good. I don't do membrane work, obviously...
These awful-looking loops in the newer structure are very consistent with your loop lengths. I can tell you that in my experience (and that of others) that neither CCD or KIC loop closure work well on loops beyond 12 residues. I am not satisfied with the performance of any loop modeling I've seen for loops beyond 12. Mike Tyka's LoopHash stuff *might* work but it's not released at the moment anyway.
I agree with you that the "loops" in the newer structure look terrible.
I haven't run the membrane_homology_modeling demo that I assume you are using as a base, but if you look at the loops file, you see 6 viable loops, one N-terminus, and one huge insertion - I'm willing to guess Rosetta does a moderately poor job on those latter two. (Have you tried the demo to see what sort of results it makes?)
I noticed that if I reduce -select_best_loop_from to 10 I get more negatively scored models. With -select_best_loop_from 1 I get even more. What is a reasonable value for this flag?
I think 1 is perfectly reasonable. The flag lets Rosetta do "extra sampling" - basically extra outer-outermost Monte Carlo loops - around some types of loop sampling. If loop sampling cheaply occurs in the middle of an expensive protocol, then the flag becomes a slider to control how many sheap loop samples to do within the context of the broader expensive model. If it's causing problems, we should definitely jettison it in favor of a value of 1.