I started working with my non-canonical residue REL and found out that it doesn't bond to other residues.
Shall I write some other atoms in the end of mol file?
I am afraid only visualization can help here so I attach the mol file.
Also I attach the created with this mol file params file and one molecule with REL residue.
Thank you in advance,
If you want a bonded residue, the params file has to have outgoing/incoming bonds (usually CONN# atoms, look at one of the regular residue types for comparison). The documentation I beleive you were looking at earlier (https://www.rosettacommons.org/demos/latest/public/relax_around_chemically_bound_ligand/README ) explains how to add connections. I think the molfile_to_params_polymer that you already tried (that wasn't working) is supposed to do it automatically. In this recent thread, I don't think the user ever got it working, but I explain what needs to happen to make a modified serine, which may be relevant advice here. https://www.rosettacommons.org/node/9672
We are having an email argument about this that I'll summarize for you once it's finished. In the meantime: you're trying to make a retinal-lysine conjugate? That would have 12 rotatable torsions? Are you going to be trying to pack that? 12 torsions is really too many for rotamer-based packing.
Yes, I'm trying to make a retinal-lysine conjugate.
Ideally, I want rosetta to know that it is a very conjugate and rigid molecule. But I didn't think about rotamer libraries.
By the way, you've mentioned that molfile_to_params_polymer can do the connections. Maybe I can find the working code somewhere?
I assume the lys-retinal is tightly packed anyway, in which case packing it isn't likely to be a problem.
If you want to make your lys-retinal polymeric manually, instructions are in that demo. Briefly, just figure out which of those atom names maps to which real atoms in lysine, and then add the UPPER and LOWER bits copied from the lysine parameters accordingly.
I've asked where the working molfile_to_params that can handle polymers is, but I haven't gotten a response.
So you mean I shall just find the lysine params file and copy the UPPER and LOWER parts of the file and the paste it in the REL.params?
But what is the name of the lysine params file here?
Rosetta's main parameters are at Rosetta/main/database/chemical/residue_type_sets/fa_standard/. Lysine specifically is at residue_types/l-caa/LYS.params. You can't quite copy it, because your REL has numbered atom names, but LYS uses the traditional CA, CB, CG, etc. You need to figure out which of your atoms is CA, etc, before pasting UPPER and LOWER in, because the UPPER and LOWER lines refer to atoms by those names. For example, for UPPER:
ICOOR_INTERNAL UPPER 149.999985 63.799984 1.328685 C CA N
ICOOR_INTERNAL O 179.999969 59.200024 1.231015 C CA UPPER
The UPPER connection (the upstream residues)'s placement depends on C, CA, and N. The backbone O atom depensd on C, CA, and UPPER. Your REL has C1, C2, C3 not CA, C, etc, so you'll need to figure out which are which.
The easiest way to figure out which is why is overlay them in PyMOL with atom name labeling turned on (I beleive you did have an output PDB with REL in it?)
Well, renaming is not a problem, I've already done it.
So I need just to copy paste these lines (2 for UPPER and 2 for LOWER) in my REL.params?
Read all of LYS.params and copy anything using UPPER and LOWER (there's one other place setting the UPPER CONNECT and LOWER CONNECT or similar, from memory). Then we see if it still fails :)
Well, it still doesn't connect with other residues. I guess the problem is occuring because in REL.params I have the full Retinal + Lysine amino acid.
It is not a radical, there is an -OH group on the C-terminal atom and extra -H on the N-terminal atom. I think that's the reason why rosetta cannot create bonds.
But on what stage shall I delete them? Because when I tryed to equilibrate it in Gaussian without these atom groups, everything bacame not realistic.
"there is an -OH group on the C-terminal atom and extra -H on the N-terminal atom"
Delete those atoms from your parameters file, then.
Another possible problem if Rosetta isn't drawing the bonds is that your input PDB says the REL should be its own chain. Is it in-chain and in-order with the rest of the peptide backbone? For example, if REL is attached to a LYS at position 65, does the REL-LYS residue in the PDB occur between 64 and 66, with the same chain ID, with the same name as the params file says it will have?
Well, I've tried to delete all the lines containing these 3 extra atoms from the params file. But then I get the error message:
"Error setting proton chi: Chi to set as proton chi does not exist." What does it mean -- proton chi? Maybe I can define another Chi as proton one?
What about the occurance of REL in the chain -- yes, everything is allright with it, the chain ID, the place and the atom numeration is right.
Try commenting the proton CHI lines out.
Proton CHIs are rotatable protons whose positions are important, like serine's alcohol proton. Methyl rotation isn't so important, but getting alcohol protons right is very important.
Well, yes, I commented the line.
And then I changed in the params file type from LIGAND to POLYMER and everything works now!
And one more question.
Previous time I had REL coord right matching its future position. Is there any way to guess this position using Rosetta, so that I won't need to find it manually
using VMD? Because I tryed to do all the same for a situation when REL molecule wasnt't near and nothing worked.
I'm not sure I understand the question. I think you are saying "the REL residue in my PDB had garbage coordinates. How do I get Rosetta to fix the coordinates when building in the residue?"
It's certainly something Rosetta can do in the broader sense, but not something that's straightforward from command line. Assuming the INTERNAL coordinates for the residue were correct, you'd use an internal C++ function that's called something like "apply_ideal_geometry_at_polymer_bond" or similar (I forget the exact name). If the internal coordinates weren't already correct you'd also need to make it flexible, treating it either as a ligand or building the already-discussed somewhat prohibitive rotamer library.
This is the sort of thing that's not too hard to describe and probably simple with PyRosetta but unlikely to be exposed through an existing command line. There might be a mode of the Remodel application that will do it. The loop modeling application's "build_initial" flag is a reliable way to fix your "bad starting coordinates" problem, but not a great tool in the middle of a tightly packed protein.
Well, yes, I'm sure about my internal coords, the only thing I want is to put it in my residue without intersections. Because even if it won't be that perfect,
I'll do lot's of relaxation and refinement anyway. The only thing I want is not to put it manually, because I want to create an automized protocol for building my structures.
So as I understand, there is no easy way to just put mol or pdb coords of my ligand, define the upper and lower connectivities, corresponding residues?
Well, there's no easy way that I know of off the top of my head. As I said, the -build_initial mode of loop modeling can do this in some fashion, and Remodel can probably do it too. A PyRosetta script wouldn't be hard either. There's a documentation event next week, I'll put this on their docket to see if one of them knows of an easy solution.
(Another way to look at this is, if nobody writing the code has needed this exact behavior before, there probably isn't a way to do it immediately from command line).
Here are a few more things that are probably ultimately automatable:
A) use the loop modeling method with build_initial and an alanine in place, then use sed to replace the alanine atom names into the REL atom names for the backbone (delete the sidechain) - this will cause rosetta to build in the sidechain in the right position once the backbone is right.
B) You can get the backbone atoms aligned with pair_align in pymol (probably in a scriptable fashion)? Again, if the backbone is fine, Rosetta can build in the sidechain from your params file.
Well, I think option A suits me better.
So, if I understand right, I shall do smth like:
main/souce/bin/loopmodel.default.linuxgccrelease -in:file:s 1e12_save.pdb -in:file:fullatom -loops:build_initila -loops: loopfile X.loop
With my loop file (if REL is residue number 219) containing just one this line:
LOOP 218 220
And with alanine instead of lysine of rel in 1e12_save.pdb
And then I shall delete all the atoms except O=C-CA-N2- and rename the resname in REL.
And then run for example relax application with REL.params and it will build the whole sidechain for me.
Is everything right?
There is some workflow similar to this that will work. Obviously I haven't tested this recently.
I guess if there's a lysine there to start, that will work fine instead of alanine - it doesn't really matter since you're getting rid of the sidechain anyway.
Look for the "loop length changing" subdemo inside the AnchoredDesign demo for an example of doing this with a normal loop.
You are likely to need a fragments file to get this to work, but the fragments don't need to be valid / relevant (in other words, just literally any fragments file that can be read in will work - it isn't going to use fragments for your problem anyway).
What was there in your input structure to begin with?
Well, I have just full lysine in my input structure.
I tryed the following now: just renamed LYS in REL (with all the atom names), and ran loopmodel.default.linuxgccrelease with -extra_res_fa REL.params and options:
-loops::frag_sizes 9 3
-loops::frag_files 1E12_9.frags 1E12_3.frags
However, I got an ERROR message:
[ERROR] Exception caught by JobDistributor for job 1E12_ops1E12_ros_0001
[ERROR] EXCN_utility_exit has been thrown from: src/core/util/SwitchResidueTypeSet.cc line: 324
ERROR: core::util::switch_to_residue_type_set fails
You need the centroid residue type for your REL as well as your existing fullatom params if you are trying to loop model with actual REL (that's why I said lysine or alanine or whatever. I know that's a question we've answered in these forums before - if you can't find it, post again and I'll look.
I don't quite understand what you're doing at this point - if there's already a lysine in your structure, why do you need loop modeling to place a REL there? You wanted to do relax, right? Just go straight into relax with REL and you don't need to mess with centroids.
Loop modeling goes through "centroid mode" where each residue is represented by a reduced-sidechain representation. Your params file is only for fullatom, so it doesn't work for centroid mode. The quick and dirty fix is to just copy the lysine centroid parameter file from the database (https://www.rosettacommons.org/docs/latest/scripting_documentation/RosettaScripts/RosettaScripts) and edit it to replace the LYS therein with your residue name. It's not completely good, as the behavior of your residue type is different from lysine, but it should keep Rosetta from crashing.
That said, I'm not completely sure what you're trying to do with the loop modeling.
If it's just changing the sequence of the position to your residue, then loop modeling is probably overkill. Rosetta will fill in any missing atoms (and will discard any extra atoms). Often, all you need to do is edit the three letter code in the input PDB to be your code, and as long as the backbone atoms are present and match the atom names of the backbone, Rosetta should fill in the rest of the residue.
Loop modeling should only be useful in this case if you don't actually have any backbone atoms for that positions. (e.g. if you have a deleted loop, or missing density in the input structure.) In that case, I might suggest just doing the loop modeling with that position being an alanine or a regular lysine, and then doing the text editing substitution to your new amino acid later. (I think this is what Steven was suggesting.)
What you might be concerned about is the conformation which Rosetta builds when it builds the missing sidechain. The default conformation for this sidechain is encoded in the ICOOR lines in the params file. The ones for the default amino acids are deliberately bad, so that you can easily recognize them being the default builds. If you want another default, just change the torsion setting (the first number) in the ICOOR lines to the value that you want. Note that that will really only affect the sidechain and the newly-built atoms - all the existing atoms (e.g. the backbone) will be kept in their original locations when Rosetta fills in missing atoms.
At any rate, if there are heavy atoms missing in the sidechain, Rosetta will repack that sidechain. That is, it will pick the lowest energy rotamer from the rotamer library you provided for that residue. If there's no rotamer library (e.g. you didn't specify one in the params file), only the default conformation (the ICOOR one) will be availible, so that's the one which will be used.
Well, now I have my REL residue in rhosopsin.
1) The problem is that when I run relax application, where Rosetta builds the left part of my REL residue, even 5 cycles of FastRelax are enough to absolutely
distort the area near retinlal in my molecule in a very bad way.
Yes, if I used it then only in Rosetta apps, I would just run another app without using fastrelax.
But the problem is that sometimes I want to use Rosetta to just attach retinal, as I really don't know another easy and free way to attach it.
And if there might be the changes in the structure, they should be very very small.
In what app I can build my REL without extra problems?
2) Another problem is: when I have LYS pre-built in structure, and Rosetta adds only retinal part, there randomly occurs problem in the junction point.
I want NZ hydrogen from LYS and H40 hydrogen from RET to be absolutely trans, absolutely 180 degrees from one another. However, it doesn't happen always.
I know, it is already a very long discussion, and I want to thank you for your advices and help.
"I know, it is already a very long discussion, and I want to thank you for your advices and help. "
You're welcome. It's long because you're doing complex stuff. Getting noncanonicals right is "easy" for maybe 5 people in the world.
Problem 1): Do you mean you want to JUST build retinal with NO other changes to the protein? If that's what you meant, use fixbb with the flag -packing:prevent_repacking, or use the score or score_jd2 applications with whatever flags they need to force PDB output. That will run the PDB loading process (where retinal gets built in) without doing much or any later modification.
problem 2): This is probably an issue with the params file (in other words, the ICOOR_INTERNAL lines for the hydrogens in question aren't specified quite right). I don't know how to actually fix that, though. The first column in H40's line is a pseduodihedral defined on the atoms at the end of the line - play with that value to see if you can get it to position it the way you want? I'm especially uncertain because you say "sometimes" - this type of placement ought to have been deterministic.