Skip to content

Commit 30b2d92

Browse files
KevKeatingKevKeating
authored
SKETCH-2782: get_residue_number_for_new_monomer now uses one residue number per monomer, even for nucleic acids (#372)
- Linked Case: SKETCH-2782 ### Description Based on [this discussion](#367 (comment)), I've modified `get_residue_number_for_new_monomer` to use one residue number per monomer, even for nucleic acids. ### Testing Done Added unit tests coverage for the new logic.
1 parent 6cd608c commit 30b2d92

3 files changed

Lines changed: 197 additions & 49 deletions

File tree

‎src/schrodinger/sketcher/model/mol_model.cpp‎

Lines changed: 50 additions & 49 deletions
Original file line numberDiff line numberDiff line change
@@ -712,12 +712,11 @@ get_monomer_type_from_chain_type(const rdkit_extensions::ChainType chain_type,
712712
}
713713

714714
/**
715-
* @return a set of all residue numbers for monomers of the specified type that
716-
* are in the same polymer as the given atom.
715+
* @return a set of all residue numbers for monomers that are in the same
716+
* polymer as the given atom.
717717
*/
718718
static std::unordered_set<int>
719-
get_all_residue_numbers_of_monomer_type_in_polymer(
720-
const MonomerType monomer_type, const RDKit::Atom* const atom_in_polymer)
719+
get_all_residue_numbers_in_polymer(const RDKit::Atom* const atom_in_polymer)
721720
{
722721
// we make a copy of the molecule so we can flag it as monomeric
723722
auto mol = atom_in_polymer->getOwningMol();
@@ -728,68 +727,70 @@ get_all_residue_numbers_of_monomer_type_in_polymer(
728727
std::unordered_set<int> residue_numbers;
729728
for (auto atom_idx : atom_idxs) {
730729
auto* atom = mol.getAtomWithIdx(atom_idx);
731-
if (get_monomer_type(atom) == monomer_type) {
732-
auto res_num = rdkit_extensions::get_residue_number(atom);
733-
residue_numbers.insert(res_num);
734-
}
730+
auto res_num = rdkit_extensions::get_residue_number(atom);
731+
residue_numbers.insert(res_num);
735732
}
736733
return residue_numbers;
737734
}
738735

739-
/**
740-
* Determine the appropriate residue number to use for a new monomer that will
741-
* be connected to an existing monomer.
742-
* @param res_name The residue name of the monomer to be added
743-
* @param chain_type The type of the monomer to be added
744-
* @param new_monomer_ap_name The new monomer's attachment point that will be
745-
* used to connect it to bound_to_monomer
746-
* @param bound_to_monomer The existing monomer that the new monomer will be
747-
* bound to
748-
* @return the residue number, which is guaranteed to be:
749-
* - non-negative (since get_residue_number returns unsigned ints)
750-
* - unique amongst all residues with the same monomer type (This allows,
751-
* e.g., an RNA base and sugar to have the same residue number, but ensures
752-
* that peptide monomers are uniquely numbered.)
753-
*/
754-
static int
755-
get_residue_number_for_new_monomer(const std::string_view res_name,
756-
const rdkit_extensions::ChainType chain_type,
757-
const std::string_view new_monomer_ap_name,
758-
const RDKit::Atom* const bound_to_monomer)
736+
int get_residue_number_for_new_monomer(
737+
const std::string_view res_name,
738+
const rdkit_extensions::ChainType chain_type,
739+
const std::string_view new_monomer_ap_name,
740+
const RDKit::Atom* const bound_to_monomer)
759741
{
760742
using rdkit_extensions::ChainType;
761743

762744
auto monomer_type = get_monomer_type_from_chain_type(chain_type, res_name);
763-
auto existing_res_nums = get_all_residue_numbers_of_monomer_type_in_polymer(
764-
monomer_type, bound_to_monomer);
745+
auto existing_res_nums =
746+
get_all_residue_numbers_in_polymer(bound_to_monomer);
747+
748+
// First, determine what the "ideal" residue number would be based on the
749+
// offset from the residue number of the monomer we're binding to. For
750+
// example, the ideal residue number of a new C-terminal peptide is one
751+
// higher than the bound peptide, and the ideal residue number of a new
752+
// N-terminal peptide is one lower than the bound peptide. We refer to it as
753+
// an "ideal" residue number since we don't know whether we'll actually be
754+
// able to use that number; it may be non-positive, which isn't allowed, or
755+
// it may already be taken by an existing monomer. In those cases, we'll
756+
// just use one more than the highest current residue number, since we know
757+
// that's available and valid.
765758

766759
int res_num_offset = 1;
767-
auto bound_to_monomer_chain_type =
768-
rdkit_extensions::getChainType(*bound_to_monomer);
769760
if (chain_type == ChainType::PEPTIDE &&
770-
bound_to_monomer_chain_type == ChainType::PEPTIDE &&
771-
new_monomer_ap_name == ap_model_name_for(PeptideAP::N)) {
761+
new_monomer_ap_name == ap_model_name_for(PeptideAP::C)) {
762+
// a new C terminal peptide residue
772763
res_num_offset = -1;
773-
} else if (chain_type == ChainType::RNA &&
774-
bound_to_monomer_chain_type == ChainType::RNA) {
775-
if (monomer_type == MonomerType::NA_SUGAR &&
776-
new_monomer_ap_name == ap_model_name_for(NASugarAP::THREE_PRIME)) {
777-
// sugars can link to the previous residue
778-
res_num_offset = -1;
779-
} else if (!((monomer_type == MonomerType::NA_PHOSPHATE &&
780-
new_monomer_ap_name ==
781-
ap_model_name_for(NAPhosphateAP::TO_NEXT_SUGAR)) ||
782-
(monomer_type == MonomerType::NA_BASE &&
783-
new_monomer_ap_name == NA_BASE_AP_PAIR))) {
784-
// phosphates and bases can link to the next residue, but all other
785-
// linkages are probably within the same residue
786-
res_num_offset = 0;
764+
} else if (chain_type == ChainType::RNA) {
765+
if (monomer_type == MonomerType::NA_SUGAR) {
766+
if (new_monomer_ap_name ==
767+
ap_model_name_for(NASugarAP::THREE_PRIME)) {
768+
// a new sugar bound to the next phosphate. We skip a number
769+
// since the base is typically given the number immediately
770+
// after the sugar
771+
res_num_offset = -2;
772+
} else if (new_monomer_ap_name ==
773+
ap_model_name_for(NASugarAP::ONE_PRIME)) {
774+
// a new sugar bound to its base
775+
res_num_offset = -1;
776+
}
777+
} else if (monomer_type == MonomerType::NA_PHOSPHATE) {
778+
if (new_monomer_ap_name ==
779+
ap_model_name_for(NAPhosphateAP::TO_PREV_SUGAR)) {
780+
res_num_offset = 2;
781+
} else if (new_monomer_ap_name ==
782+
ap_model_name_for(NAPhosphateAP::TO_NEXT_SUGAR)) {
783+
res_num_offset = -1;
784+
}
787785
}
788786
}
789787
auto bound_to_res_num =
790788
rdkit_extensions::get_residue_number(bound_to_monomer);
791789
int new_res_num = bound_to_res_num + res_num_offset;
792-
if (new_res_num < 0 || existing_res_nums.contains(new_res_num)) {
790+
791+
if (new_res_num <= 0 || existing_res_nums.contains(new_res_num)) {
792+
// the ideal residue number isn't valid or isn't available. Note that
793+
// residue numbers of 0 confuse HELM generation, so we avoid them.
793794
new_res_num = *std::ranges::max_element(existing_res_nums) + 1;
794795
}
795796
return new_res_num;

‎src/schrodinger/sketcher/model/mol_model.h‎

Lines changed: 22 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -1714,5 +1714,27 @@ SKETCHER_API void add_text_to_mol_model(
17141714
const std::optional<RDGeom::Point3D> position = std::nullopt,
17151715
const bool recenter_view = true);
17161716

1717+
/**
1718+
* Determine the appropriate residue number to use for a new monomer that will
1719+
* be connected to an existing monomer via a standard connection
1720+
* @param res_name The residue name of the monomer to be added
1721+
* @param chain_type The type of the monomer to be added, which must be the same
1722+
* as the chain type of bound_to_monomer (since otherwise the new residue should
1723+
* be a separate chain)
1724+
* @param new_monomer_ap_name The new monomer's attachment point that will be
1725+
* used to connect it to bound_to_monomer
1726+
* @param bound_to_monomer The existing monomer that the new monomer will be
1727+
* bound to
1728+
* @return the residue number, which is guaranteed to be:
1729+
* - positive
1730+
* - unique amongst all residues in the same chain (The HELM reader and writer
1731+
* expect one monomer per residue number, even for nucleic acids.)
1732+
*/
1733+
SKETCHER_API int
1734+
get_residue_number_for_new_monomer(const std::string_view res_name,
1735+
const rdkit_extensions::ChainType chain_type,
1736+
const std::string_view new_monomer_ap_name,
1737+
const RDKit::Atom* const bound_to_monomer);
1738+
17171739
} // namespace sketcher
17181740
} // namespace schrodinger

‎test/schrodinger/sketcher/model/test_mol_model.cpp‎

Lines changed: 125 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -26,6 +26,7 @@
2626
#include <rdkit/GraphMol/FileParsers/FileParsers.h>
2727
#include <rdkit/GraphMol/FileParsers/MolSupplier.h>
2828
#include <rdkit/GraphMol/MolTransforms/MolTransforms.h>
29+
#include <rdkit/GraphMol/MonomerInfo.h>
2930
#include <rdkit/GraphMol/QueryAtom.h>
3031
#include <rdkit/GraphMol/QueryBond.h>
3132
#include <rdkit/GraphMol/ROMol.h>
@@ -4800,6 +4801,130 @@ BOOST_AUTO_TEST_CASE(test_addBoundMonomer_RNA)
48004801
BOOST_TEST(helm == "RNA1{R(C)P.R(A)P.R(G)P}$$$$V2.0");
48014802
}
48024803

4804+
/**
4805+
* Test get_residue_number_for_new_monomer() for a peptide
4806+
*/
4807+
BOOST_AUTO_TEST_CASE(test_get_residue_number_for_new_monomer_peptide)
4808+
{
4809+
QUndoStack undo_stack;
4810+
TestMolModel model(&undo_stack);
4811+
4812+
// create an amino acid with a residue number of 10
4813+
auto mol = rdkit_extensions::to_rdkit("PEPTIDE1{A}$$$$V2.0", Format::HELM);
4814+
auto* res_info = static_cast<RDKit::AtomPDBResidueInfo*>(
4815+
mol->getAtomWithIdx(0)->getMonomerInfo());
4816+
res_info->setResidueNumber(10);
4817+
model.addMol(*mol);
4818+
auto peptide = model.getMol()->getAtomWithIdx(0);
4819+
BOOST_REQUIRE(rdkit_extensions::get_residue_number(peptide) == 10u);
4820+
4821+
// check residue numbers for the previous and next peptide
4822+
auto prev_peptide_num = get_residue_number_for_new_monomer(
4823+
"C", ChainType::PEPTIDE, ap_model_name_for(PeptideAP::C), peptide);
4824+
auto next_peptide_num = get_residue_number_for_new_monomer(
4825+
"C", ChainType::PEPTIDE, ap_model_name_for(PeptideAP::N), peptide);
4826+
BOOST_TEST(prev_peptide_num == 9u);
4827+
BOOST_TEST(next_peptide_num == 11u);
4828+
4829+
// add a new amino acid so we can check what happens when the ideal residue
4830+
// numbers is already assigned to an existing monomer
4831+
model.addBoundMonomer("C", ChainType::PEPTIDE, {1.0, 0.0, 0.0},
4832+
ap_model_name_for(PeptideAP::N), peptide,
4833+
ap_model_name_for(PeptideAP::C));
4834+
peptide = model.getMol()->getAtomWithIdx(0);
4835+
next_peptide_num = get_residue_number_for_new_monomer(
4836+
"C", ChainType::PEPTIDE, ap_model_name_for(PeptideAP::N), peptide);
4837+
BOOST_TEST(next_peptide_num == 12u);
4838+
}
4839+
4840+
/**
4841+
* Test get_residue_number_for_new_monomer() for a nucleic acid sugar
4842+
*/
4843+
BOOST_AUTO_TEST_CASE(test_get_residue_number_for_new_monomer_na_sugar)
4844+
{
4845+
QUndoStack undo_stack;
4846+
TestMolModel model(&undo_stack);
4847+
4848+
// create a sugar with a residue number of 10
4849+
auto mol = rdkit_extensions::to_rdkit("RNA1{R}$$$$V2.0", Format::HELM);
4850+
auto* res_info = static_cast<RDKit::AtomPDBResidueInfo*>(
4851+
mol->getAtomWithIdx(0)->getMonomerInfo());
4852+
res_info->setResidueNumber(10);
4853+
model.addMol(*mol);
4854+
auto sugar = model.getMol()->getAtomWithIdx(0);
4855+
BOOST_REQUIRE(rdkit_extensions::get_residue_number(sugar) == 10u);
4856+
4857+
// check residue numbers for a base and phosphates bound to the sugar
4858+
auto base_num = get_residue_number_for_new_monomer(
4859+
"A", ChainType::RNA, ap_model_name_for(NA_BASE_AP_N1_9), sugar);
4860+
auto next_phos_num = get_residue_number_for_new_monomer(
4861+
"P", ChainType::RNA, ap_model_name_for(NAPhosphateAP::TO_PREV_SUGAR),
4862+
sugar);
4863+
auto prev_phos_num = get_residue_number_for_new_monomer(
4864+
"P", ChainType::RNA, ap_model_name_for(NAPhosphateAP::TO_NEXT_SUGAR),
4865+
sugar);
4866+
BOOST_TEST(base_num == 11u);
4867+
BOOST_TEST(next_phos_num == 12u);
4868+
BOOST_TEST(prev_phos_num == 9u);
4869+
}
4870+
4871+
/**
4872+
* Test get_residue_number_for_new_monomer() for a nucleic acid phosphate
4873+
*/
4874+
BOOST_AUTO_TEST_CASE(test_get_residue_number_for_new_monomer_na_phosphate)
4875+
{
4876+
QUndoStack undo_stack;
4877+
TestMolModel model(&undo_stack);
4878+
4879+
// create a phosphate with a residue number of 10
4880+
auto mol = rdkit_extensions::to_rdkit("RNA1{P}$$$$V2.0", Format::HELM);
4881+
auto* res_info = static_cast<RDKit::AtomPDBResidueInfo*>(
4882+
mol->getAtomWithIdx(0)->getMonomerInfo());
4883+
res_info->setResidueNumber(10);
4884+
model.addMol(*mol);
4885+
auto phosphate = model.getMol()->getAtomWithIdx(0);
4886+
BOOST_REQUIRE(rdkit_extensions::get_residue_number(phosphate) == 10u);
4887+
4888+
// check residue numbers for sugars bound to the phosphate
4889+
auto prev_sugar_num = get_residue_number_for_new_monomer(
4890+
"R", ChainType::RNA, ap_model_name_for(NASugarAP::THREE_PRIME),
4891+
phosphate);
4892+
auto next_sugar_num = get_residue_number_for_new_monomer(
4893+
"R", ChainType::RNA, ap_model_name_for(NASugarAP::FIVE_PRIME),
4894+
phosphate);
4895+
// the base is expected to be numbered in between the previous sugar and
4896+
// this phosphate, so we skip 9 and instead assign 8 to the previous sugar
4897+
BOOST_TEST(prev_sugar_num == 8u);
4898+
// the next base should get 12 (i.e. be numbered after the next sugar), so
4899+
// we don't need to skip a number between this phosphate and the next sugar
4900+
BOOST_TEST(next_sugar_num == 11u);
4901+
}
4902+
4903+
/**
4904+
* Test get_residue_number_for_new_monomer() for a nucleic acid base
4905+
*/
4906+
BOOST_AUTO_TEST_CASE(test_get_residue_number_for_new_monomer_na_base)
4907+
{
4908+
QUndoStack undo_stack;
4909+
TestMolModel model(&undo_stack);
4910+
4911+
// create a base with a residue number of 10
4912+
auto mol = rdkit_extensions::to_rdkit("RNA1{A}$$$$V2.0", Format::HELM);
4913+
auto* res_info = static_cast<RDKit::AtomPDBResidueInfo*>(
4914+
mol->getAtomWithIdx(0)->getMonomerInfo());
4915+
res_info->setResidueNumber(10);
4916+
model.addMol(*mol);
4917+
auto base = model.getMol()->getAtomWithIdx(0);
4918+
BOOST_REQUIRE(rdkit_extensions::get_residue_number(base) == 10u);
4919+
4920+
// check residue number for a sugar bound to the base
4921+
auto sugar_num = get_residue_number_for_new_monomer(
4922+
"R", ChainType::RNA, ap_model_name_for(NASugarAP::ONE_PRIME), base);
4923+
// the base should be numbered after the sugar, so the sugar should get 9
4924+
// (instead of 11)
4925+
BOOST_TEST(sugar_num == 9);
4926+
}
4927+
48034928
/**
48044929
* Confirm that combining two peptide chains via a standard backbone connection
48054930
* produces HELM output with only a single chain

0 commit comments

Comments
 (0)