Hi, I am trying to learn protein-ligand docking using complex 3dau (protein+NADPH+MTX)
My first task is to calculate the binding energy of cofactor binding to a protein.
I did the following:
1.I relaxed the protein+cofactor+ligand (whole complex).
2.Then i generated the params file for the ligand and the cofactor. (Using the coordinates from the output of relaxed pdb)
3.Then i appended the generated cofactor pdb(from2) to protein pdb(protein coordinates got from otput of relax).
4.Then i used the dock application:
./ligand_dock.linuxgccrelease -database ~/rosetta-3.4/rosetta_database/ -s ~/input/3dau_plus_nap.pdb -in:file:extra_res_fa ~/input/NAP.params
I got a silent file as a output: silent.out; score= -507.698.
I calculated the score of the (only)protein pdb using ./score.linuxgccrelease : -330.170.
As i said, I have to calculate the binding energy of cofactor binding to the protein.
Will the difference in score of the protein and docked protein+cofactor give the binding energy?
If not, how do the binding energy calculation?
Also what are the ouputs of a typical dock run?
What are the information we can get from the silent.out?
How to convert the silent.out to .pdb.
First, I would switch steps 1 + 2. That way when you relax the whole complex, you have will have something there. From my experience, if Rosetta doesn't understand a cofactor or ligand, it will ignore it and it won't even be in the output PDB.
For the relax, you may want to edit the scorefunction weights file (or copy it and rename it) (rosetta_database/scoring/weights/standard.wts and replace the statistical term fa_pair with the coulombic term hack_elec at about the same weights. This will allow electrostatic scoring of your cofactor + ligand during relax. When you relax the complex, make sure to use the option -constrain_relax_to_start_coords. For info on how to use the your modified scorefunction: https://www.rosettacommons.org/manuals/archive/rosetta3.4_user_guide/d5/....
To get the binding energy, if your cofactor + ligand are different chains from the protein chains, which I think they should be anyway, you can use Steven's Interface Analyzer application. See https://www.rosettacommons.org/manuals/archive/rosetta3.4_user_guide/df/...
You may also want to heed Rocco's advice and just relax the cofactor+protein:
'A final note - if you're trying to emulate a "blind" docking, you don't want to include the docking partner when you pre-relax the protein. There's a good chance if you include the partner that you'll over-optimize to the docked state, resulting in results that are artificially better, and wouldn't be applicable to a blind docking case.'
Though, I didn't realize that constrain_relax_to_start_coords will renumber the PDB in Rosetta 3.4. So, just so you know -
In the ligand docking output, there should be annotations of the binding energy (along with other parameters) in the scoring lines. This will be the energy of the docked ligand to the cofactor+protein complex, though, instead of the energy of the ligand+cofactor to the protein.
Rescoring in the apo state is a reasonable way of finding the binding energy. You have to be a little careful of energy terms which involve the internal energy of the ligand/cofactor. The standard ligand docking scorefunction doesn't really have them (there's often a negligible contribution from fa_intra_rep), but if you're using other terms like the MM terms, the ligand may have appreciable internal energy. The best bet is to just count the edges in the energy graph between the ligand and the protein (how ligand docking normally does it), or compare the bound complex energy to that of a complex where the ligand and cofactor have been translated far away. (Generally, at least a 10 Ang space between atoms means no interaction energy.)
That said, a -170 REU binding energy for ligand+cofactor is not reasonable. I'm guessing the difference is mainly in the energy functions use. Ligand docking uses it's own scoring function, and the score application typically uses score12. These won't result in comparable weights. I would recommend using the same scorefunction for each - e.g. rescore both the holo and apo complexes with score12, with the ligand docking scorefunction ("-score:weights ligand"), or with the modified score function Jared suggests. Note that depending on flags, ligand docking does some energy function manipulation, so output from the ligand docking application won't necessarily be the same as that from the score application with the ligand weights.
To extract the information from the output, you can use the extract_atomtree_diffs application. See the ligand docking manual (https://www.rosettacommons.org/manuals/archive/rosetta3.4_user_guide/d4/...) for details.