You are here

Problems when running Rosetta VIP with cofactors

15 posts / 0 new
Last post
Problems when running Rosetta VIP with cofactors
#1

Dear Rosetta users,

I'm trying to run Rosetta VIP on a heme protein. I started by relaxing the structure according to the protocol in the "relax_around_chemically_bound_ligand" demo. My chemically bound ligand would be the heme group, to a histidine. For this, I created *.params files for the heme group and the covalently bound histidine, and created a contraints file according to the protocol. The relaxation seemed to work fine (log.txt available), as Rosetta was able to make the connection between heme and his and apply the constraints.
The problem is when I try to run Rosetta VIP on the relaxed protein (log file available). I give Rosetta both *.paramas files, the contranints file, a resfile (resfile says "don't touch heme nor his") and the starting structure. My flag file for VIP reads:

-database /.../rosetta_database/
#-ignore_unrecognized_res
-s startingStructure.pdb
-cst_fa_file contraints.cst
-resfile resfile.txt
-extra_res_fa HEM.params # *.parmas file for HEME
-extra_res_fa HIS.params # *.parmas file for modified HIS
-cp:cutoff 6.0
-sasa_calculator_probe_radius 1.0
-run:silent
-cp:print_intermediate_pdbs
-cp:local_relax
-cp:relax_mover "relax"
-cp:print_reports
-mute protocols.simple_moves protocols::checkpoint protocols.relax core.scoring.etable basic.io.database core.io.pdb core.chemical protocols.packstat #core.pack

The first problem appears when core.conformation.Conformation attempts to connect the Fe atom from the heme to the N atom from the histidine:

core.conformation.Conformation: Failed to find a residue connection for residue 154 with connection point 1

(residue 154 is the heme)

What surprised me is that core.conformation.Conformation is the same thinf that makes the bond properly in the relax application (this can be seen in the relax log file).

Then, VIP keeps running and I have a second error, which is the following:

core.pack.task.ResfileReader: APOLAR mode read for residue 154 which has been instructed to use non-canonical amino acids.
ERROR:: Exit from: src/core/pack/task/ResfileReader.cc line: 1485

VIP stops running there. The error seems (IMO) to be related to the params file, because I can comment the resfile line in the flag file and nothing changes.

Let me know if you want me to post the complete log file from VIP.

Thanks,
Benjamin

Post Situation: 
Fri, 2013-04-26 13:31
Basantab

"The first problem appears when core.conformation.Conformation attempts to connect the Fe atom from the heme to the N atom from the histidine:"

The PDB file format says nothing about where chemical bonds are and are not present, and what is polymer versus ligand. (Technically you can use CONECT records, but nobody actually does). So, Rosetta is forced to guess, and this particular instance is guessing wrong. Conformation is a class partly responsible for this sort of data. VIP does not use the standard PDB reader, for whatever reason. If relax worked, and VIP didn't, you can try replacing this line from vip.cc:

core::io::pdb::build_pose_from_pdb_as_is( in_pose, option[ OptionKeys::in::file::s ]().vector().front() );

with something like this:

core::import_pose::pose_from_pdb( in_pose, option[ OptionKeys::in::file::s ]().vector().front() );

You may need to add a #include for core/import_pose/import_pose.hh.

"because I can comment the resfile line in the flag file and nothing changes."
This doesn't make any sense - you say that removing the resfile flag from your options file still leaves a resfile-related crash? Can you confirm that? Maybe also post your resfile?

Mon, 2013-04-29 07:58
smlewis

Hi smlewis,

Regarding the first error: the relax protocol interpreted the bond between Fe and N correcly with the bonding information contained only in the *.params files, I did not use CONECT in the starting pdb for relax. That's why I gessed it would work in VIP too. I will try your solution.

What I meant by "because I can comment the resfile line in the flag file and nothing changes." is that I can make a new flag file leaving the "-resfile resfile.txt" line out (by adding "#" at the begining) and get the same error, as if rosetta was not taking it into account in the first place. On the contrary, if I leave out the *.params file, rosetta doesn't recognize the heme and stops running (which was rather spected). I was not aware that the error from ResfileReader.cc was strictly related to the resfile option, still, the error is there with or whithout including the resfile flag. Does it make sense now?

Another piece of information: The "resfile" error occurs when rosetta is making some rotamers, you can see this in the attached log file.

This is my resfile.txt:

#Leave out heme and histidine 93

start

93 A NATRO
154 A NATRO

The intention of this is to not allow VIP to move or mutate the heme group or histidine 93, which binds it.

Benjamin

Mon, 2013-04-29 09:04
Basantab

First error: if that thing I suggested doesn't work (and I think it will), let's try making a single his+heme params type instead (I'm surprised that's not what the demo already did).

Second error:

UGH. VIP is totally ignoring the input resfile...and confusing the issue by needlessly using resfile-labeled machinery (the APOLAR setting) to get the behavior it wants. You'll have to modify the code to fix this.

rosetta3.4/rosetta_source/src/protocols/vip/VIP_Mover.cc, function try_point_mutants:

This if statement:

if( j != void_mutatables[aa] ){
task->nonconst_residue_task(j).prevent_repacking();
}

Needs to be modified to exclude positions 93 and 154. If 93 and 154 are resids (continuous numbering from 1) as well as PDB numbers (the error suggests this!), then just:

if( j != void_mutatables[aa] || j=93 || j=154 ){

will work as a quick hack. If you are uncomfortable with C++, or need a more robust hack, we can work it out.

Mon, 2013-04-29 09:42
smlewis

I see that in the solution that you proposed for the first error you refer to the VIP.cc file, but I only found VIP_mover.cc and VIP_report.cc in the protocols folder, in addition to that, I found the "core::io::pdb::build_pose_from_pdb_as_is..." line in the VIP_mover.cc file. I guess that making the change in VIP_mover is what you meant, right?

I can try doing the quick hack. Correct me if I'm wrong, but both changes would require to recompile VIP. If so, I would like to avoid recompiling all the applications. Can you give me some instructions to recompile VIP alone in a new folder in the /bin directory? something like /bin/vipAlt. Or should I go to the building forum for this?

Thanks,
Benjamin

**EDIT**
I would like to put the modified *.cc files in a separate folder and tell scons to build it as a new app, if that's not too complicated.

Mon, 2013-04-29 12:36
Basantab

Scons should be smart enough to only build those items which it needs. Recompiling after compiling should take less time than trying to recompile from scratch in a separate folder.

You can create a separate application by copying the relevant application code to a new name, and then adding the appropriate line to rosetta_source/src/apps.src.settings. Recompiling should create the new application in the bin directory.

That's only going to work if the changes you're making are limited to the rosetta_source/src/apps directory. If you're making changes in the rosetta_source/src/protocols directory, then all application are going to share the changes. If that isn't going to work for you, I probably would recommend making a copy of the rosetta_source directory tree, making the changes there, and recompiling in the copied rosetta_source directory.

Mon, 2013-04-29 18:11
rmoretti

Hi,

I applied the changes you proposed and compiled in a separate folder:

-Issue with the HIS-HEME bond: Changing the original line in vip.cc for "core::import_pose::pose_from_pdb( in_pose, option[ OptionKeys::in::file::s ]().vector().front() );" didn't fix the problem, I just got a new output line under the previous error message saying something similar:

core.conformation.Conformation: Failed to find a residue connection for residue 154 with connection point 1 //same line as before
core.import_pose.import_pose: Can't find a chemical connection for residue 154 HEM //new line after changing vip.cc

-Avoiding HIS and HEME repack: The quick fix worked and I have some output structures to look at. I can give you a hand on finding a general solution for this problem, mainly in the testing end.

Going back to the HIS-HEME bond issue, shouldn't I be able to keep the geometry of the HIS-HEME bond by using just the constraints file and avoiding repack of HIS and HEME? I tried this and Rosetta doesn't seem to be taking the constraints.cst file into account, as I can see in the distorted angles and distances between HIS and HEME in the output structures I was able to get. Also, I don't see any output regarding constraints as I did during the relax protocol (output muted during VIP run: protocols.simple_moves protocols::checkpoint protocols.relax core.scoring.etable basic.io.database core.io.pdb core.chemical protocols.packstat core.scoring.rms_util).

New issue 1: I realized that relaxing the structure before running VIP on it is unnecessary and might affect the outcome of the algorithm, so I added hydrogens via an on-line server leaving the HEME out and then added hydrogens to the separate HEME with pymol. I ran the mol_to_params.py script on the *.mol version of the HEME and got the *.params file and HEM_0001.pdb. To get the input structure I just pasted the content of HEM_0001.pdb before the "END" tag of my hydrogenated protein. With this input structure I ran Rosetta VIP with the following flag:

-database /.../rosetta_database/
#-ignore_unrecognized_res
-s nonRelaxedInput.pdb
-resfile resfile.txt
-extra_res_fa HEM.params
-cst_fa_file constraints.cst
-cp:cutoff 6.0
-sasa_calculator_probe_radius 1.0
-run:silent
-cp:print_intermediate_pdbs
-cp:local_relax
-cp:relax_mover "relax"
-cp:print_reports
-mute protocols.simple_moves protocols::checkpoint protocols.relax core.scoring.etable basic.io.database core.io.pdb core.chemical protocols.packstat core.pack

With this flag and the repack fix VIP works and I can get some structures, but I also get a lot of output from core.scoring.rms_util, always the same line, many times:

core.scoring.rms_util: WARNING: In CA_rmsd, residue range 1 to 154 requested but only 153 protein CA atoms found.

I couldn't find the source of this warning, but I didn't get it when using the pre-relaxed structure. As I said, VIP runs anyway, but should I worry about it?

New Issue 2: I'm getting some mutations in my final_mutant.pdb that are not listed among the accepted mutations in the program output. They do appear in the reports.txt file and, apparently, they are the last mutations of each run that are supposedly rejected but for some reason appear in the final pdb. This also happens when using the provided flag and input structure in the demo. A bug maybe?

Thanks,
Benjamin

Wed, 2013-05-01 21:57
Basantab

Enforcing a bond with constraints should be okay ... though you may have problems with atom overlap. If the hydrogens that would have been removed on bond formation still exist, (and even the bonded atoms that come in closer proximity) they'll have a high repulsive energy because of the bonded distances. You may need to have very high weights on the bonded interaction to overcome this repulsion.

Regarding the warning, it's just telling you it's ignoring the HEME for the purposes of rmsd, as it doesn't have a Calpha atom. It's entirely expected when using ligands, and isn't a problem as long as you're fine with having the reported rmsd be across the protein only.

Thu, 2013-05-02 11:25
rmoretti

"A bug maybe":

Sounds like it...I've brought it to the attention of the authors.

Fri, 2013-05-03 14:24
smlewis

Hi rmoretti,

I would like to make sure Rosetta is finding and using the constraints file, what part of the standard output should I take out of the "-mute" tag to see it?
Also, would you please point me in the right direction on how to manipulate the weight of the bonded interaction?

Thanks,
Benjamin

Thu, 2013-05-02 13:58
Basantab

The best way to ensure Rosetta is using constraints is to write a deliberately garbage constraint and ensure it crashes - add "TEST_FAIL_LINE" or some other garbage and make sure it errors. None of those muted tracers look likely to be muting constraints.

Brief inspection of the VIP code suggests it never loads constraints anyway (as you suggested). Adding constraints is nontrivial, because the code uses a lot of different scorefunctions in different places.

If you want to try constraints, try adding "core::scoring::constraints::add_constraints_from_cmdline_to_pose( in_pose );" after the pose creation line you already modified in apps/vip.cc, (this may require a #included header core/scoring/constraints/util.hh)

and adding

core::scoring::constraints::add_constraints_from_cmdline_to_scorefunction( some_scorefunction_variable_name ); (with the RIGHT VARIABLE NAME)

immediately after all instances of create_score_function in apps/vip.cc and protocols/VIP_Mover.cc (and possibly other places I didn't catch).

The use of global minimization in the protocol suggests to me that you need kinematic connection between heme and his, not just a constraint. (At minimum a constraint would need to constrain distance, angle, and dihedral separately!) Are you willing to try a unified residue type that has the heme attached to the histidine? You make it the same way as you made your heme.params, except using a polymer version of molfile_to_params, which is here: ./rosetta_demos/design_with_ncaa/scripts/python/apps/public/molfile_to_params_polymer.py

This post: https://www.rosettacommons.org/content/error-recognizing-rotamer-library... has attached a similar polymer ligand - in this case it's a cysteine with a fluorophore attached. I had rotamers here too, which can be generated by the noncanonicals machinery but are a bit of a hassle.

Fri, 2013-05-03 14:40
smlewis

I went for the unified residue type option and tried to make the *.params file with molfile_to_params_polymer.py

First I had to move it to the folder where the rest of the public python apps are: /rosetta_source/src/apps/python/public so it would be able to import everything it needs. molfile_to_params_polymer.py is originally in the folder you suggested and running it from there gave this error:

File "molfile_to_params_polymer.py", line 19, in
from rosetta_py.io.mdl_molfile import *
ImportError: No module named rosetta_py.io.mdl_molfile

After moving molfile_to_params_polymer.py, I ran it with my modified .mol file (1) as input, getting the following error:

File "molfile_to_params_polymer.py", line 1995, in
sys.exit(main(sys.argv[1:]))
File "molfile_to_params_polymer.py", line 1919, in main
m = molfiles[0]
TypeError: 'generator' object has no attribute '__getitem__'

I repeated this with a .mol file that had previously worked in with molfile_to_params.py and got the same error. I wasn't able to figure out what's wrong, but it looks like there is a problem when reading the input. I wonder if moving molfile_to_params_polymer.py to a different folder is the cause of this.

(1) modified *.mol file: I followed the instructions in the readme.txt file that comes with the ncaa demo. At one point it requires to add some lines starting with "M " at the end of the *.mol file. I tried to add these lines the best I could, but it's hard to extrapolate the instructions from the example to my ncaa, i.e.: what are exactly the atoms in the POLY_IGNORE line? part of the molecule that will be erased from the .params file?
I saw at the post you cited that additional lines can be used, like "POLY_UPPER" and "POLY_LOWER", what are these lines for? Is there an alternative tutorial to use molfile_to_params_polymer.py different from the one in the demo folder? maybe one that is more extensive and less focused in the example?

Attached is my *.mol file.

Thanks,
Benjamin

Wed, 2013-05-08 07:09
Basantab

Solved:

I copied all the files from the /rosetta_demo/.../scripts folder to the corresponding subfolders in the /rosetta_source/src/python/rosetta_py folder. Then, added the "UPPER_CONNECT" and "LOWER_CONNECT" lines to the .params file. Now rosetta keeps the structure of the his-heme complex as if it was a residue.

Thu, 2013-05-16 05:59
Basantab

I was about to use Rosetta Holes to generate pdbs and score my input and output structures, but there seems to be no holes executable, as I would expect from the pdf available with the holes demo. I could not find a holes.cc file either. Should I add something to rosetta_source/src/apps.src.settings and recompile?

Thanks,
Benjamin

Fri, 2013-05-03 06:54
Basantab

holes has been in-and-out of the release a few times for reasons that escape me. I think it was in earlier releases (maybe 3.3?).

Fri, 2013-05-03 14:49
smlewis