Skip to content

SKETCH-2782: get_residue_number_for_new_monomer now uses one residue number per monomer, even for nucleic acids #372

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

Merged
KevKeating merged 4 commits into from
Jun 29, 2026
Merged

SKETCH-2782: get_residue_number_for_new_monomer now uses one residue number per monomer, even for nucleic acids #372

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
4 commits
Select commit Hold shift + click to select a range
fad48f7
SKETCH-2782: get_residue_number_for_new_monomer now uses one residue …
KevKeating Jun 11, 2026
9e7f9f9
clang-format
KevKeating Jun 11, 2026
a5a02da
fixing residue numbering and rewriting the tests so the numbering mis…
KevKeating Jun 25, 2026
6364721
adding comments explaining test numbering
KevKeating Jun 25, 2026
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
99 changes: 50 additions & 49 deletions src/schrodinger/sketcher/model/mol_model.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -712,12 +712,11 @@ get_monomer_type_from_chain_type(const rdkit_extensions::ChainType chain_type,
}

/**
* @return a set of all residue numbers for monomers of the specified type that
* are in the same polymer as the given atom.
* @return a set of all residue numbers for monomers that are in the same
* polymer as the given atom.
*/
static std::unordered_set<int>
get_all_residue_numbers_of_monomer_type_in_polymer(
const MonomerType monomer_type, const RDKit::Atom* const atom_in_polymer)
get_all_residue_numbers_in_polymer(const RDKit::Atom* const atom_in_polymer)
{
// we make a copy of the molecule so we can flag it as monomeric
auto mol = atom_in_polymer->getOwningMol();
Expand All @@ -728,68 +727,70 @@ get_all_residue_numbers_of_monomer_type_in_polymer(
std::unordered_set<int> residue_numbers;
for (auto atom_idx : atom_idxs) {
auto* atom = mol.getAtomWithIdx(atom_idx);
if (get_monomer_type(atom) == monomer_type) {
auto res_num = rdkit_extensions::get_residue_number(atom);
residue_numbers.insert(res_num);
}
auto res_num = rdkit_extensions::get_residue_number(atom);
residue_numbers.insert(res_num);
}
return residue_numbers;
}

/**
* Determine the appropriate residue number to use for a new monomer that will
* be connected to an existing monomer.
* @param res_name The residue name of the monomer to be added
* @param chain_type The type of the monomer to be added
* @param new_monomer_ap_name The new monomer's attachment point that will be
* used to connect it to bound_to_monomer
* @param bound_to_monomer The existing monomer that the new monomer will be
* bound to
* @return the residue number, which is guaranteed to be:
* - non-negative (since get_residue_number returns unsigned ints)
* - unique amongst all residues with the same monomer type (This allows,
* e.g., an RNA base and sugar to have the same residue number, but ensures
* that peptide monomers are uniquely numbered.)
*/
static int
get_residue_number_for_new_monomer(const std::string_view res_name,
const rdkit_extensions::ChainType chain_type,
const std::string_view new_monomer_ap_name,
const RDKit::Atom* const bound_to_monomer)
int get_residue_number_for_new_monomer(
const std::string_view res_name,
const rdkit_extensions::ChainType chain_type,
const std::string_view new_monomer_ap_name,
const RDKit::Atom* const bound_to_monomer)
{
using rdkit_extensions::ChainType;

auto monomer_type = get_monomer_type_from_chain_type(chain_type, res_name);

@rachelnwalker rachelnwalker Jun 22, 2026

Copy link
Copy Markdown
Collaborator

Choose a reason for hiding this comment

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

General comment but not urgent - I'm not sure whether you will always be able to determine whether something is a base/sugar/phosphate based off residue name. I'm not sure how common it will be but there are customers with modified phosphate linkers and sugars that won't adhere to a certain naming convention - you maybe need to deduce by a substructure match, or just rely completely on the attachment points used

@KevKeating KevKeating Jun 24, 2026

Copy link
Copy Markdown
Collaborator Author

Choose a reason for hiding this comment

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

Yeah, agreed. We might want to record information about the monomer type when the monomer is added, since Sketcher and/or the HELM converter would presumably be in the best position to figure that out (and the Sketcher would likely have additional information about the monomer based on the tool being used). I've filed SKETCH-2796 for this.

auto existing_res_nums = get_all_residue_numbers_of_monomer_type_in_polymer(
monomer_type, bound_to_monomer);
auto existing_res_nums =
get_all_residue_numbers_in_polymer(bound_to_monomer);

// First, determine what the "ideal" residue number would be based on the
// offset from the residue number of the monomer we're binding to. For
// example, the ideal residue number of a new C-terminal peptide is one
// higher than the bound peptide, and the ideal residue number of a new
// N-terminal peptide is one lower than the bound peptide. We refer to it as
// an "ideal" residue number since we don't know whether we'll actually be
// able to use that number; it may be non-positive, which isn't allowed, or
// it may already be taken by an existing monomer. In those cases, we'll
// just use one more than the highest current residue number, since we know
// that's available and valid.

int res_num_offset = 1;
auto bound_to_monomer_chain_type =
rdkit_extensions::getChainType(*bound_to_monomer);
if (chain_type == ChainType::PEPTIDE &&
bound_to_monomer_chain_type == ChainType::PEPTIDE &&
new_monomer_ap_name == ap_model_name_for(PeptideAP::N)) {
new_monomer_ap_name == ap_model_name_for(PeptideAP::C)) {
// a new C terminal peptide residue
res_num_offset = -1;
} else if (chain_type == ChainType::RNA &&
bound_to_monomer_chain_type == ChainType::RNA) {
if (monomer_type == MonomerType::NA_SUGAR &&
new_monomer_ap_name == ap_model_name_for(NASugarAP::THREE_PRIME)) {
// sugars can link to the previous residue
res_num_offset = -1;
} else if (!((monomer_type == MonomerType::NA_PHOSPHATE &&

@rachelnwalker rachelnwalker Jun 22, 2026

Copy link
Copy Markdown
Collaborator

Choose a reason for hiding this comment

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

Wouldn't these offsets be +1 and +2? Totally possible I am not understanding the MonomerType enum correctly - is this a base connecting to a sugar then a phosphate connecting to the same sugar? if the sugar is 0, I would think that the base is 1 then the phosphate is 2

@KevKeating KevKeating Jun 24, 2026 •
edited
Loading

Copy link
Copy Markdown
Collaborator Author

Choose a reason for hiding this comment

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

Yeah, it looks like a bunch of the numbers here are off. Good catch! I'll fix that and update the tests, both to check for the correct results and hopefully to make those types of issues more obvious. (I also temporarily tweaked monomer drawing locally so it includes the residue numbers to verify that the current numbering is wrong, since I kept getting myself massively confused trying to reason through the numbering here.)

new_monomer_ap_name ==
ap_model_name_for(NAPhosphateAP::TO_NEXT_SUGAR)) ||
(monomer_type == MonomerType::NA_BASE &&
new_monomer_ap_name == NA_BASE_AP_PAIR))) {
// phosphates and bases can link to the next residue, but all other
// linkages are probably within the same residue
res_num_offset = 0;
} else if (chain_type == ChainType::RNA) {
if (monomer_type == MonomerType::NA_SUGAR) {
if (new_monomer_ap_name ==
ap_model_name_for(NASugarAP::THREE_PRIME)) {
// a new sugar bound to the next phosphate. We skip a number
// since the base is typically given the number immediately
// after the sugar
res_num_offset = -2;
Comment on lines +766 to +771

@ethan-schrodinger ethan-schrodinger Jun 18, 2026

Copy link
Copy Markdown
Collaborator

Choose a reason for hiding this comment

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

Should we be worried about this "typically"?

@KevKeating KevKeating Jun 24, 2026

Copy link
Copy Markdown
Collaborator Author

Choose a reason for hiding this comment

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

Not really. This would only be an issue if the user is intentionally leaving an abasic site (i.e. a missing base), and even then, the only result of this would be that the residues aren't numbered sequentially, which shouldn't have any real consequences. If we want a guarantee that we have the ideal residue numbering, i.e. starting at one and continuing sequentially in the correct order, we'll need to renumber the entire strand after it's been built. That's probably not worth worrying about unless we know of something downstream that's requires it.

} else if (new_monomer_ap_name ==
ap_model_name_for(NASugarAP::ONE_PRIME)) {
// a new sugar bound to its base
res_num_offset = -1;
}
} else if (monomer_type == MonomerType::NA_PHOSPHATE) {
if (new_monomer_ap_name ==
ap_model_name_for(NAPhosphateAP::TO_PREV_SUGAR)) {
res_num_offset = 2;
} else if (new_monomer_ap_name ==
ap_model_name_for(NAPhosphateAP::TO_NEXT_SUGAR)) {
res_num_offset = -1;
}
}
}
auto bound_to_res_num =
rdkit_extensions::get_residue_number(bound_to_monomer);
int new_res_num = bound_to_res_num + res_num_offset;
if (new_res_num < 0 || existing_res_nums.contains(new_res_num)) {

if (new_res_num <= 0 || existing_res_nums.contains(new_res_num)) {
// the ideal residue number isn't valid or isn't available. Note that
// residue numbers of 0 confuse HELM generation, so we avoid them.
new_res_num = *std::ranges::max_element(existing_res_nums) + 1;
}
return new_res_num;
Expand Down
22 changes: 22 additions & 0 deletions src/schrodinger/sketcher/model/mol_model.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -1714,5 +1714,27 @@ SKETCHER_API void add_text_to_mol_model(
const std::optional<RDGeom::Point3D> position = std::nullopt,
const bool recenter_view = true);

/**
* Determine the appropriate residue number to use for a new monomer that will
* be connected to an existing monomer via a standard connection
* @param res_name The residue name of the monomer to be added
* @param chain_type The type of the monomer to be added, which must be the same
* as the chain type of bound_to_monomer (since otherwise the new residue should
* be a separate chain)
* @param new_monomer_ap_name The new monomer's attachment point that will be
* used to connect it to bound_to_monomer
* @param bound_to_monomer The existing monomer that the new monomer will be
* bound to
* @return the residue number, which is guaranteed to be:
* - positive
* - unique amongst all residues in the same chain (The HELM reader and writer
* expect one monomer per residue number, even for nucleic acids.)
*/
SKETCHER_API int
get_residue_number_for_new_monomer(const std::string_view res_name,
const rdkit_extensions::ChainType chain_type,
const std::string_view new_monomer_ap_name,
const RDKit::Atom* const bound_to_monomer);

} // namespace sketcher
} // namespace schrodinger
125 changes: 125 additions & 0 deletions test/schrodinger/sketcher/model/test_mol_model.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -26,6 +26,7 @@
#include <rdkit/GraphMol/FileParsers/FileParsers.h>
#include <rdkit/GraphMol/FileParsers/MolSupplier.h>
#include <rdkit/GraphMol/MolTransforms/MolTransforms.h>
#include <rdkit/GraphMol/MonomerInfo.h>
#include <rdkit/GraphMol/QueryAtom.h>
#include <rdkit/GraphMol/QueryBond.h>
#include <rdkit/GraphMol/ROMol.h>
Expand Down Expand Up @@ -4800,6 +4801,130 @@ BOOST_AUTO_TEST_CASE(test_addBoundMonomer_RNA)
BOOST_TEST(helm == "RNA1{R(C)P.R(A)P.R(G)P}$$$$V2.0");
}

/**
* Test get_residue_number_for_new_monomer() for a peptide
*/
BOOST_AUTO_TEST_CASE(test_get_residue_number_for_new_monomer_peptide)
{
QUndoStack undo_stack;
TestMolModel model(&undo_stack);

// create an amino acid with a residue number of 10
auto mol = rdkit_extensions::to_rdkit("PEPTIDE1{A}$$$$V2.0", Format::HELM);
auto* res_info = static_cast<RDKit::AtomPDBResidueInfo*>(
mol->getAtomWithIdx(0)->getMonomerInfo());
res_info->setResidueNumber(10);
model.addMol(*mol);
auto peptide = model.getMol()->getAtomWithIdx(0);
BOOST_REQUIRE(rdkit_extensions::get_residue_number(peptide) == 10u);

// check residue numbers for the previous and next peptide
auto prev_peptide_num = get_residue_number_for_new_monomer(
"C", ChainType::PEPTIDE, ap_model_name_for(PeptideAP::C), peptide);
auto next_peptide_num = get_residue_number_for_new_monomer(
"C", ChainType::PEPTIDE, ap_model_name_for(PeptideAP::N), peptide);
BOOST_TEST(prev_peptide_num == 9u);
BOOST_TEST(next_peptide_num == 11u);

// add a new amino acid so we can check what happens when the ideal residue
// numbers is already assigned to an existing monomer
model.addBoundMonomer("C", ChainType::PEPTIDE, {1.0, 0.0, 0.0},
ap_model_name_for(PeptideAP::N), peptide,
ap_model_name_for(PeptideAP::C));
peptide = model.getMol()->getAtomWithIdx(0);
next_peptide_num = get_residue_number_for_new_monomer(
"C", ChainType::PEPTIDE, ap_model_name_for(PeptideAP::N), peptide);
BOOST_TEST(next_peptide_num == 12u);
}

/**
* Test get_residue_number_for_new_monomer() for a nucleic acid sugar
*/
BOOST_AUTO_TEST_CASE(test_get_residue_number_for_new_monomer_na_sugar)
{
QUndoStack undo_stack;
TestMolModel model(&undo_stack);

// create a sugar with a residue number of 10
auto mol = rdkit_extensions::to_rdkit("RNA1{R}$$$$V2.0", Format::HELM);
auto* res_info = static_cast<RDKit::AtomPDBResidueInfo*>(
mol->getAtomWithIdx(0)->getMonomerInfo());
res_info->setResidueNumber(10);
model.addMol(*mol);
auto sugar = model.getMol()->getAtomWithIdx(0);
BOOST_REQUIRE(rdkit_extensions::get_residue_number(sugar) == 10u);

// check residue numbers for a base and phosphates bound to the sugar
auto base_num = get_residue_number_for_new_monomer(
"A", ChainType::RNA, ap_model_name_for(NA_BASE_AP_N1_9), sugar);
auto next_phos_num = get_residue_number_for_new_monomer(
"P", ChainType::RNA, ap_model_name_for(NAPhosphateAP::TO_PREV_SUGAR),
sugar);
auto prev_phos_num = get_residue_number_for_new_monomer(
"P", ChainType::RNA, ap_model_name_for(NAPhosphateAP::TO_NEXT_SUGAR),
sugar);
BOOST_TEST(base_num == 11u);
BOOST_TEST(next_phos_num == 12u);
BOOST_TEST(prev_phos_num == 9u);
}

/**
* Test get_residue_number_for_new_monomer() for a nucleic acid phosphate
*/
BOOST_AUTO_TEST_CASE(test_get_residue_number_for_new_monomer_na_phosphate)
{
QUndoStack undo_stack;
TestMolModel model(&undo_stack);

// create a phosphate with a residue number of 10
auto mol = rdkit_extensions::to_rdkit("RNA1{P}$$$$V2.0", Format::HELM);
auto* res_info = static_cast<RDKit::AtomPDBResidueInfo*>(
mol->getAtomWithIdx(0)->getMonomerInfo());
res_info->setResidueNumber(10);
model.addMol(*mol);
auto phosphate = model.getMol()->getAtomWithIdx(0);
BOOST_REQUIRE(rdkit_extensions::get_residue_number(phosphate) == 10u);

// check residue numbers for sugars bound to the phosphate
auto prev_sugar_num = get_residue_number_for_new_monomer(
"R", ChainType::RNA, ap_model_name_for(NASugarAP::THREE_PRIME),
phosphate);
auto next_sugar_num = get_residue_number_for_new_monomer(
"R", ChainType::RNA, ap_model_name_for(NASugarAP::FIVE_PRIME),
phosphate);
// the base is expected to be numbered in between the previous sugar and
// this phosphate, so we skip 9 and instead assign 8 to the previous sugar
BOOST_TEST(prev_sugar_num == 8u);
// the next base should get 12 (i.e. be numbered after the next sugar), so
// we don't need to skip a number between this phosphate and the next sugar
BOOST_TEST(next_sugar_num == 11u);
Comment on lines +4888 to +4900

@ethan-schrodinger ethan-schrodinger Jun 25, 2026

Copy link
Copy Markdown
Collaborator

Choose a reason for hiding this comment

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

I think it would be worth adding comments here (and on the test below) that briefly explain why we should expect these specific residue numbers.

@KevKeating KevKeating Jun 25, 2026

Copy link
Copy Markdown
Collaborator Author

Choose a reason for hiding this comment

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

Yeah, that seems helpful. Added.

}

/**
* Test get_residue_number_for_new_monomer() for a nucleic acid base
*/
BOOST_AUTO_TEST_CASE(test_get_residue_number_for_new_monomer_na_base)
{
QUndoStack undo_stack;
TestMolModel model(&undo_stack);

// create a base with a residue number of 10
auto mol = rdkit_extensions::to_rdkit("RNA1{A}$$$$V2.0", Format::HELM);
auto* res_info = static_cast<RDKit::AtomPDBResidueInfo*>(
mol->getAtomWithIdx(0)->getMonomerInfo());
res_info->setResidueNumber(10);
model.addMol(*mol);
auto base = model.getMol()->getAtomWithIdx(0);
BOOST_REQUIRE(rdkit_extensions::get_residue_number(base) == 10u);

// check residue number for a sugar bound to the base
auto sugar_num = get_residue_number_for_new_monomer(
"R", ChainType::RNA, ap_model_name_for(NASugarAP::ONE_PRIME), base);
// the base should be numbered after the sugar, so the sugar should get 9
// (instead of 11)
BOOST_TEST(sugar_num == 9);
}

/**
* Confirm that combining two peptide chains via a standard backbone connection
* produces HELM output with only a single chain
Expand Down