Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

Rmsd in smallmol docking #95

Open
wants to merge 2 commits into
base: main
Choose a base branch
from
Open

Rmsd in smallmol docking #95

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
f544ff4
Changed rmsd calculation to automorphic and added informative tracer …
Jul 16, 2024
c7621f0
Added rmsd filter criteria to HighResDocker
Jul 17, 2024
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
56 changes: 52 additions & 4 deletions source/src/protocols/ligand_docking/HighResDocker.cc
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -44,6 +44,7 @@
#include <core/pack/rotamer_set/UnboundRotamersOperation.hh>

#include <core/scoring/ScoreFunction.hh>
#include <core/scoring/rms_util.hh>

// Package Headers
#include <core/pose/Pose.hh>
Expand Down Expand Up @@ -127,7 +128,9 @@ HighResDocker::HighResDocker(HighResDocker const & that):
num_cycles_(that.num_cycles_),
repack_every_Nth_(that.repack_every_Nth_),
score_fxn_(that.score_fxn_),
movemap_builder_(that.movemap_builder_)
movemap_builder_(that.movemap_builder_),
rmsd_max_(that.rmsd_max_),
check_rmsd_(that.check_rmsd_)
{
initialize_from_options();
}
Expand Down Expand Up @@ -185,6 +188,12 @@ HighResDocker::parse_my_tag(
if ( tag->hasOption("resfile") ) {
resfile_= tag->getOption<std::string>("resfile");
}

/// RMSD ///
if ( tag->hasOption("rmsd") ) {
check_rmsd_ = true;
rmsd_max_ = tag->getOption<core::Real>("rmsd");
}
}


Expand Down Expand Up @@ -282,11 +291,26 @@ HighResDocker::apply(core::pose::Pose & pose) {
score_fxn_->score( pose ); // without this neither of the movers below were working
// I believe that this may have been related to adding constraints incorrectly at other places in my code.
// Rigid body exploration
utility::vector1<protocols::moves::MoverOP> rigid_body_movers= create_rigid_body_movers(pose);
utility::vector1<protocols::moves::MoverOP> rigid_body_movers = create_rigid_body_movers(pose);

core::Size ligand_index = 0;
for ( core::Size i(1), imax(pose.total_residue()); i<=imax; ++i ) {
core::conformation::Residue const & residue( pose.residue(i) );
if ( residue.is_ligand() ) {
ligand_index = i;
Comment on lines +299 to +300
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

There can theoretically be multiple ligands in the pose, only one of which (not necessarily the last) is the one which is being modified. Unfortunately, things are indirect enough that I don't think there's a good way of getting just the primary ligand from the info provided, so there likely isn't a need to change things. Just make sure to document the issue thoroughly. -- Though I might suggest breaking this out into a private function, which might make updating things later easier.

}
}
core::conformation::Residue original_ligand(pose.residue(ligand_index));
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Prefer using a ResidueOP and .clone() to make a copy. I don't think it makes a difference here, but it does in other contexts where there's multiple derived types, it plays better with subtypes (the approach here would lose subtype info.)

TR.Debug << "Save residue number " << ligand_index << " as ligand." << std::endl;

core::Size rejected_rmsd_moves = 0;

for ( core::Size cycle = 1; cycle <= num_cycles_; ++cycle ) {

TR.Debug << "Start cycle " << cycle << ".";
if ( repack_every_Nth_ != 0 && cycle % repack_every_Nth_ == 1 ) TR.Debug << " Repacking!";
TR.Debug << std::endl;

core::pack::task::PackerTaskOP packer_task;

packer_task = make_packer_task(pose, all_residues_);
Expand Down Expand Up @@ -320,10 +344,22 @@ HighResDocker::apply(core::pose::Pose & pose) {
min_mover->apply(pose);
}

if ( check_rmsd_ ) {
core::Real rmsd = core::scoring::automorphic_rmsd( original_ligand, pose.residue(ligand_index), /*superimpose=*/false );
TR.Debug << "Current rmsd is " << rmsd << ". Max rmsd: " << rmsd_max_ << std::endl;
if (rmsd > rmsd_max_){
pose = monteCarlo->last_accepted_pose();
rejected_rmsd_moves++;
}
}
monteCarlo->boltzmann( pose );

}

if (check_rmsd_) {
TR << "Rejected " << rejected_rmsd_moves << " of " << num_cycles_ << " pose updates due to rmsd limitations." << std::endl;
}

// keep the best structure we found, not the current one
monteCarlo->show_scores();
monteCarlo->recover_low(pose);
Expand Down Expand Up @@ -628,7 +664,8 @@ void HighResDocker::provide_xml_schema( utility::tag::XMLSchemaDefinition & xsd
+ XMLSchemaAttribute::required_attribute("movemap_builder", xs_string, "Name of a previously defined MoveMaoBuilder.")
+ XMLSchemaAttribute("resfile", xs_string, "Name (path to) the resfile.")
+ XMLSchemaAttribute("cycles", xsct_non_negative_integer, "Number of cycles to run.")
+ XMLSchemaAttribute("repack_every_Nth", xsct_non_negative_integer, "Perform side chain repacking every Nth cycle.");
+ XMLSchemaAttribute("repack_every_Nth", xsct_non_negative_integer, "Perform side chain repacking every Nth cycle.")
+ XMLSchemaAttribute("rmsd", xsct_real, "Maximum RMSD to be sampled away from the starting position.");
protocols::moves::xsd_type_definition_w_attributes( xsd, mover_name(), "Randomly connects a fragment from the library to the growing ligand.", attlist );
}

Expand Down Expand Up @@ -670,5 +707,16 @@ void HighResDocker::set_all_residues(bool input)
all_residues_ = input;
}

void HighResDocker::set_rmsd_limit(core::Real rmsd_limit)
{
check_rmsd_ = true;
rmsd_max_ = rmsd_limit;
}

void HighResDocker::disable_rmsd_limit()
{
check_rmsd_ = false;
}

} //namespace ligand_docking
} //namespace protocols
} //namespace protocols
7 changes: 7 additions & 0 deletions source/src/protocols/ligand_docking/HighResDocker.hh
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -95,6 +95,10 @@ public:

void set_all_residues(bool input);

void set_rmsd_limit(core::Real rmsd_limit);

void disable_rmsd_limit();

public: //For CitationManager:

/// @brief Provide the citation.
Expand Down Expand Up @@ -153,6 +157,9 @@ private:
bool allow_repacking_ = true;
bool all_residues_ = false;

core::Real rmsd_max_ = 99999.999;
bool check_rmsd_ = false;

};

//void setup_native_residue_favoring(core::pose::Pose & pose, core::pack::task::PackerTaskOP task);
Expand Down
14 changes: 9 additions & 5 deletions source/src/protocols/ligand_docking/Transform.cc
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -27,6 +27,7 @@
#include <core/scoring/ScoreType.hh>
#include <core/scoring/constraints/ConstraintSet.hh>
#include <core/scoring/constraints/ConstraintIO.hh>
#include <core/scoring/rms_util.hh>
#include <core/conformation/UltraLightResidue.hh>
#include <core/pose/Pose.hh>
#include <core/pose/PDBInfo.hh>
Expand Down Expand Up @@ -236,6 +237,7 @@ void Transform::apply(core::pose::Pose & pose)
core::Size accepted_moves = 0;
core::Size rejected_moves = 0;
core::Size outside_grid_moves = 0;
core::Size rmsd_rejects = 0;

// core::pose::Pose best_pose(pose);

Expand Down Expand Up @@ -385,6 +387,7 @@ void Transform::apply(core::pose::Pose & pose)
if ( check_rmsd_ && !check_rmsd(original_ligand, ligand_residue) ) { //reject the new pose
ligand_residue = last_accepted_ligand_residue;
rejected_moves++;
rmsd_rejects++;
} else if ( probability < 1 && numeric::random::rg().uniform() >= probability ) { //reject the new pose
ligand_residue = last_accepted_ligand_residue;
rejected_moves++;
Expand Down Expand Up @@ -431,6 +434,10 @@ void Transform::apply(core::pose::Pose & pose)
TR.Warning << " and a box size of at least " << utility::Real2string(recommended_box_size( accept_ratio ),1) << " are recommended." << std::endl;
}
}
if ( rmsd_rejects > 0 ) {
core::Real rmsd_reject_ratio = (core::Real)rmsd_rejects / ( (core::Real)accepted_moves + (core::Real)rejected_moves );
TR << "Moves rejected for going beyond max RMSD (" << transform_info_.rmsd << "): " << rmsd_rejects << " " << rmsd_reject_ratio << std::endl;
}

protocols::jd2::add_string_real_pair_to_current_job("Transform_accept_ratio", accept_ratio);

Expand Down Expand Up @@ -612,12 +619,9 @@ bool Transform::check_rmsd(core::conformation::UltraLightResidue const & start,
{
debug_assert(start.natoms() == current.natoms());

core::Real total_distance =0.0;
for ( core::Size atomno = 1; atomno <= start.natoms(); ++atomno ) {
total_distance += start[atomno].distance(current[atomno]);
}
core::Real rmsd = core::scoring::automorphic_rmsd( *start.residue(), *current.residue(), /*superimpose=*/false );
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Part of the purpose of UltraLightResidue is to not have to update the coordinates of the underlying conformation::Residue and take the overhead of that. As such, the residue() coordinates are set at initialization and not updated -- your change here isn't doing what you want it to (it's using stale coordinates).


core::Real rmsd = sqrt(total_distance/start.natoms());
TR.Debug << "Current rmsd " << rmsd << ". Limit: " << transform_info_.rmsd << std::endl;

if ( rmsd <= transform_info_.rmsd ) {
return true;
Expand Down