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 2 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
83 changes: 37 additions & 46 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,67 +727,59 @@ 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
// residue number of the monomer we're binding to
int res_num_offset = 1;
Comment thread
ethan-schrodinger marked this conversation as resolved.
Outdated
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)) {
// 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
} 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_BASE &&
new_monomer_ap_name == ap_model_name_for(NA_BASE_AP_N1_9)) {
res_num_offset = -1;
} else if (monomer_type == MonomerType::NA_PHOSPHATE &&
new_monomer_ap_name ==
ap_model_name_for(NAPhosphateAP::TO_PREV_SUGAR)) {
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;
}
}
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 the ideal residue number is not available, just use one more than the
// highest current residue number, since we know that's available
if (new_res_num < 0 || existing_res_nums.contains(new_res_num)) {
Comment thread
ethan-schrodinger marked this conversation as resolved.
Outdated
new_res_num = *std::ranges::max_element(existing_res_nums) + 1;
}
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:
* - non-negative (since get_residue_number returns unsigned ints)
* - 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
121 changes: 121 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 @@ -4800,6 +4800,127 @@ 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 peptide chains.
*
* The function picks a residue number based on the residue number of the
* monomer being bound to, an offset that depends on the chain type and the
* attachment point used, and the set of residue numbers already in the chain.
* For peptides:
* - connecting via the new monomer's N attachment point (R1) uses an offset
* of -1
* - any other attachment point uses the default offset of +1
* If the resulting "ideal" number is negative or already taken, the function
* instead falls back to one more than the highest existing residue number.
*
* Note that the function only reads the residue number and polymer of
* bound_to_monomer; it does not require the attachment points to form a valid
* connection, so we can drive each branch independently.
*/
BOOST_AUTO_TEST_CASE(test_get_residue_number_for_new_monomer_peptide)
{
using rdkit_extensions::get_residue_number;
const auto N_AP = ap_model_name_for(PeptideAP::N);
const auto C_AP = ap_model_name_for(PeptideAP::C);

QUndoStack undo_stack;
TestMolModel model(&undo_stack);
// residue numbers {1, 2, 3} for the three alanines
add_text_to_mol_model(model, "PEPTIDE1{A.A.A}$$$$V2.0");
const auto* res1 = model.getMol()->getAtomWithIdx(0);
const auto* res3 = model.getMol()->getAtomWithIdx(2);
BOOST_REQUIRE(get_residue_number(res1) == 1u);
BOOST_REQUIRE(get_residue_number(res3) == 3u);

// Offset of -1 (new monomer's N attachment point) applied to residue 1. The
// ideal number (1 - 1 = 0) is available, so it is used as-is. This also
// covers the documented guarantee that the residue number is non-negative
// at its lower boundary.
BOOST_TEST(get_residue_number_for_new_monomer("G", ChainType::PEPTIDE, N_AP,
res1) == 0);

// Default offset of +1 (any non-N attachment point) applied to residue 3,
// the highest in the chain. The ideal number (3 + 1 = 4) is available.
BOOST_TEST(get_residue_number_for_new_monomer("G", ChainType::PEPTIDE, C_AP,
res3) == 4);

// Fallback when the ideal number is already taken: the default +1 offset
// applied to residue 1 yields 2, which is already used, so the function
// falls back to one more than the highest existing number (3 + 1 = 4).
BOOST_TEST(get_residue_number_for_new_monomer("G", ChainType::PEPTIDE, C_AP,
res1) == 4);
}

/**
* Test get_residue_number_for_new_monomer() for RNA chains, where the offset
* depends on both the monomer type (sugar/base/phosphate, derived from the new
* monomer's residue name) and the attachment point used:
* - new sugar via its 3' attachment point (R2): offset -2 (a number is
* skipped to leave room for the sugar's base)
* - new sugar via its 1' attachment point (R3): offset -1
* - new base via its N1/N9 attachment point (R1): offset -1
* - new phosphate via its attachment point to the previous sugar (R1):
* offset -1
* - anything else: default offset of +1
*/
BOOST_AUTO_TEST_CASE(test_get_residue_number_for_new_monomer_rna)
{
using rdkit_extensions::get_residue_number;

QUndoStack undo_stack;
TestMolModel model(&undo_stack);
// RNA1{R(A)P} gives a sugar (residue 1), a base (residue 2), and a
// phosphate (residue 3)
add_text_to_mol_model(model, "RNA1{R(A)P}$$$$V2.0");
const auto* sugar = model.getMol()->getAtomWithIdx(0);
const auto* base = model.getMol()->getAtomWithIdx(1);
const auto* phosphate = model.getMol()->getAtomWithIdx(2);
BOOST_REQUIRE(get_residue_number(sugar) == 1u);
BOOST_REQUIRE(get_residue_number(base) == 2u);
BOOST_REQUIRE(get_residue_number(phosphate) == 3u);

// New base via its N1/N9 attachment point (R1): offset -1 applied to the
// sugar (residue 1). Ideal number 1 - 1 = 0 is available.
BOOST_TEST(get_residue_number_for_new_monomer(
"A", ChainType::RNA, ap_model_name_for(NA_BASE_AP_N1_9),
sugar) == 0);

// New phosphate via its attachment point to the previous sugar (R1): offset
// -1 applied to the sugar (residue 1). Ideal number 1 - 1 = 0 is available.
BOOST_TEST(get_residue_number_for_new_monomer(
"P", ChainType::RNA,
ap_model_name_for(NAPhosphateAP::TO_PREV_SUGAR),
sugar) == 0);

// New sugar via its 1' attachment point (R3): offset -1 applied to the
// sugar (residue 1). Ideal number 1 - 1 = 0 is available.
BOOST_TEST(get_residue_number_for_new_monomer(
"R", ChainType::RNA, ap_model_name_for(NASugarAP::ONE_PRIME),
sugar) == 0);

// New sugar via its 3' attachment point (R2): offset -2 applied to the base
// (residue 2). Ideal number 2 - 2 = 0 is available.
BOOST_TEST(get_residue_number_for_new_monomer(
"R", ChainType::RNA,
ap_model_name_for(NASugarAP::THREE_PRIME), base) == 0);

// New phosphate via its attachment point to the next sugar (R2): not a
// special case, so the default +1 offset applies. Applied to the phosphate
// (residue 3, the highest), ideal number 3 + 1 = 4 is available.
BOOST_TEST(get_residue_number_for_new_monomer(
"P", ChainType::RNA,
ap_model_name_for(NAPhosphateAP::TO_NEXT_SUGAR),
phosphate) == 4);

// Fallback to the negative-number branch: the new sugar's 3' offset of -2
// applied to the sugar (residue 1) gives an ideal number of 1 - 2 = -1.
// Since that is negative, the function falls back to one more than the
// highest existing number (3 + 1 = 4).
BOOST_TEST(get_residue_number_for_new_monomer(
"R", ChainType::RNA,
ap_model_name_for(NASugarAP::THREE_PRIME), sugar) == 4);
}

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