From f544ff46908a76a1d03269acc255603183e89b60 Mon Sep 17 00:00:00 2001 From: Paul Eisenhuth Date: Tue, 16 Jul 2024 14:38:54 +0200 Subject: [PATCH 1/2] Changed rmsd calculation to automorphic and added informative tracer output --- source/src/protocols/ligand_docking/Transform.cc | 14 +++++++++----- 1 file changed, 9 insertions(+), 5 deletions(-) diff --git a/source/src/protocols/ligand_docking/Transform.cc b/source/src/protocols/ligand_docking/Transform.cc index f5cbc8e078..842bc7d0ad 100644 --- a/source/src/protocols/ligand_docking/Transform.cc +++ b/source/src/protocols/ligand_docking/Transform.cc @@ -27,6 +27,7 @@ #include #include #include +#include #include #include #include @@ -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); @@ -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++; @@ -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); @@ -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 ); - 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; From c7621f0281958d0f122422384038c0f958369bd7 Mon Sep 17 00:00:00 2001 From: Paul Eisenhuth Date: Wed, 17 Jul 2024 16:10:21 +0200 Subject: [PATCH 2/2] Added rmsd filter criteria to HighResDocker --- .../protocols/ligand_docking/HighResDocker.cc | 56 +++++++++++++++++-- .../protocols/ligand_docking/HighResDocker.hh | 7 +++ 2 files changed, 59 insertions(+), 4 deletions(-) diff --git a/source/src/protocols/ligand_docking/HighResDocker.cc b/source/src/protocols/ligand_docking/HighResDocker.cc index 28ff4a8229..d29d0e57d4 100644 --- a/source/src/protocols/ligand_docking/HighResDocker.cc +++ b/source/src/protocols/ligand_docking/HighResDocker.cc @@ -44,6 +44,7 @@ #include #include +#include // Package Headers #include @@ -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(); } @@ -185,6 +188,12 @@ HighResDocker::parse_my_tag( if ( tag->hasOption("resfile") ) { resfile_= tag->getOption("resfile"); } + + /// RMSD /// + if ( tag->hasOption("rmsd") ) { + check_rmsd_ = true; + rmsd_max_ = tag->getOption("rmsd"); + } } @@ -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 rigid_body_movers= create_rigid_body_movers(pose); + utility::vector1 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; + } + } + core::conformation::Residue original_ligand(pose.residue(ligand_index)); + 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_); @@ -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); @@ -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 ); } @@ -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 \ No newline at end of file diff --git a/source/src/protocols/ligand_docking/HighResDocker.hh b/source/src/protocols/ligand_docking/HighResDocker.hh index af109bcebe..ef10f4a6d9 100644 --- a/source/src/protocols/ligand_docking/HighResDocker.hh +++ b/source/src/protocols/ligand_docking/HighResDocker.hh @@ -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. @@ -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);