Hello,
I have downloaded a PDB with a chemically bound ligand (small molecule)
I have followed (in every way that I can see) all the details presented in the tutorial/demo here:
https://www.rosettacommons.org/demos/latest/public/relax_around_chemically_bound_ligand/README
I have been able to get Rosetta3 to complete a basic relax using my 1) modified PDB file along with 2) a params file and 3) a chemical constrainst file.
The resulting, relaxed pdb shows generally the correct geometry (i.e. the angles and distances are close to the crystal structure from which I measured the specified parameters) between the two atoms that should be covalently bonded (a cystine from the protein and a carbon on the ligand). However, in the relaxed structure the cystine just ends up getting protenated (not bonded).
EDIT: I ran the tutorial scripts and the result is the same--the relaxed structure leaved the lignad un-bonded--is actually the correct behavior here?
You can see by the warnings when I run relax, that the connection I'm asking for is not getting made (not sure it the other warnings are relevant)
core.import_pose.import_pose: File '4cxt_modded.pdb' automatically determined to be of type PDB
core.conformation.Conformation: [ WARNING ] missing heavyatom: CG on residue LYS 83
core.conformation.Conformation: [ WARNING ] missing heavyatom: CD on residue LYS 83
core.conformation.Conformation: [ WARNING ] missing heavyatom: CE on residue LYS 83
core.conformation.Conformation: [ WARNING ] missing heavyatom: NZ on residue LYS 83
core.conformation.Conformation: [ WARNING ] missing heavyatom: OXT on residue ASP:CtermProteinFull 132
core.conformation.Conformation: [ WARNING ] Failed to find a residue connection for residue 133 with connection point 1
core.pack.pack_missing_sidechains: packing residue number 83 because of missing atom number 6 atom name CG
Residue 133 should be my ligand (in Rosetta numbering) and my constraint file it's trying to pair it with the SG of my protein (residue 103):
AtomPair SG 103 C1 133 HARMONIC 1.828 0.01
Angle SG 103 C1 133 C2 133 HARMONIC 2.0534322 0.034906585
Angle CB 103 SG 103 C1 133 HARMONIC 1.9713319 0.034906585
I feel that I am missing some small detail that may not be spelled out in the turotiral.
I'm happy to provide whatever pther files would be helpfui but I'm what is expected here (i.e. I just don't want to dump my whole working directory ).
If someone coud please help me out, either by suggesting a detail that is not covered in the tutorial, or by letting me know what files (or other details) I should provide here to trouble shoot this.
Thanks
UPDATE: Running the script
$ROSETTA3_TOOLS/protein_tools/scripts/pdb_renumber.py
on my input pbd reveals that the ligand I cut-and-pasted into the original PDB is getting renumbered as "1" (not "133" as the tutorial would suggest) . To confirm that this was expected behavior, I ran the same renumbering sctipt on the input pdb from the tutorial, and thesame thing happened to the ligand there too (numbered as "1" rather than the"147" specified in the constrainst file).
So is this why rosetta cant find residue "133" in my original relax? Nope--I re-ran the tutorial script and it resulted in no warnings about not being able to connet the ligand residues.
So still stuck
Regarding pdb_renumber.py, by default it restarts numbering for each chain. If you want to get true Rosetta numbering, also pass the --norestart option. (Alternatively, the -renumber_pdb option passed to any Rosetta binary executable will output the renumbered version. And if you just want to check numbering, the energies table at the end of the PDB from Rosetta output will have the Rosetta numbering.)
That demo's a little bit old. The code for reading in covalently linked residues has changed a bit since then, and I'm wondering if things are cross purposes now.
Most notably, Rosetta now pays more attention to the LINK records. The example file has the residual LINK record from the original PDB, which doesn't quite work for setting up covalent connections, as the atom name and residue number for the FMN is not properly set. Additionally, Rosetta should automatically add connection information for PDBs containing LINK records, so the whole CONN editing should be unnecessary. That's also a benefit to you, as it should also automatically remove hydrogens from the connection, whereas the manual CONN editing approach doesn't change the hydrogens listed. (It will leave the hydrogens as they were in the free-ligand state. This didn't matter much for the FMN example, as the connection was to an atom without a hydrogen to begin with.)
If it's working better for you, I believe that the approach in the demo should still work. You may just need to remove the LINK records. And if the hydrogen is bothering you, you would have to manually remove it from the params file (or from the input SDF prior to generation) when you do the CONN changes.
Thank you rmoretti.
The --norestart option on renumber.py convinced me that my renumbering is in sync with Rosetta's (so the C1 of my ligand should at least be "found")
The LINK option was actually something I was testing yesterday (became suspicious that Rosetta renumbering wouldn't be able to correct the number in that LINK line). However, changing that line to match the re-numbered residue didn't change anything, so I gave up on that idea.
So I'm still stuck as the tutorial's approach doesn't produce the desired connection (with or without the LINK line). But if I'm getting what you are saying, I should go back to messing with the LINK line? Maybe even remove the CONN and ICOOR_INTERNAL lines from the params file...?
EDIT: I've now removed both the CONNECT and ICOOR_INTERNAL lines (corresponding to CONN1) from the params file and the relax runs without any errors. However I now get this warning
That's the rosetta numbered position of my ligand as specified in the LINK line (my ligand residue name is SXJ)
Still, the relax runs and I get a structure, but the output-ed relaxed structure still shows no covalent bond. So I'm assuming it's still not being interpreted correctly but I can't see what needs to be adjusted.
If I could just get a line on a working example of a pdb/params/cst file set, I could probably figure it out.
I finally got relax to run without it complaining about missing any connections. Tas was after I rebuilt the input PDB file using cleanpdb and then renumber_pdb to make sure there was nothing to surprise Rosetta (or to confuse me).
Unfortunatly, my output pdb still does not show the ligand as being bound and the inclusion of the LINK line doesn't seem to make a difference in my output socres. For example, to test out whether the LINK line was making any difference, I ran 3 relaxes (on the same starting PDB) with and then without the LINK line. The average scores for the 3 relaxes on structures with the LINK line were ~-246 REU, while the average scores for the structures without the LINK line were ~-247 REU. So not very differnet.
Then, as a final sanity check I ran relax without the CST file and (unsurprisingly) ended up with scores around ~-415 REUs.
So I'm tempted to move forward on this, rationalizing the decision becasue the CST file (which is at least reflective of the covalently bound ligand) is exerting a substantial effect on the resulting model.
Got everything working.
The LINK line does iin fact work, but the formatting/spacing of the line is critical. This seems obvious in retrospect, but no warnings or errors evert got thrown by Rosetta about an un-parsable LINK line--since we all know that Rosetta likes to warn us of everything, I was assuming that something like this would trigger a warning as well.