Generate and read scaffold networks

Cargo.toml

@@ -1,6 +1,6 @@
 [package]
 name = "rdkit"
-version = "0.4.12"
+version = "0.4.13"
 edition = "2021"
 authors = ["Xavier Lange <xrlange@gmail.com>", "Javier Pineda <javierjpineda13@gmail.com"]
 license = "MIT"
@@ -16,7 +16,7 @@ byteorder = "1.4.3"
 cxx = "1"
 flate2 = "1"
 log = "0.4"
-rdkit-sys = { path = "rdkit-sys", version = "0.4.9" }
+rdkit-sys = { path = "rdkit-sys", version = "0.4.13" }
 thiserror = "1"
 
 [dev-dependencies]

rdkit-sys/CHANGELOG.md

@@ -1,3 +1,10 @@
+0.4.13
+---
+
+ * Scaffold network generation: `create_scaffold_network`/`update_scaffold_network`, plus
+   accessors that return the nodes, counts, molCounts and edges as plain data
+ * `rdkit` gains a safe `scaffold_network_for_mol`/`scaffold_network_for_mols` wrapper
+
 0.4.0
 ---
 

rdkit-sys/Cargo.toml

@@ -1,7 +1,7 @@
 [package]
 name = "rdkit-sys"
 authors = ["Xavier Lange (xrlange@gmail.com)", "chrissly31415"]
-version = "0.4.12"
+version = "0.4.13"
 edition = "2021"
 license = "MIT"
 description = "RDKit CFFI library builder and bindings"

rdkit-sys/src/bridge/scaffold_network.rs

@@ -1,8 +1,10 @@
 #[cxx::bridge(namespace = "RDKit")]
 pub mod ffi {
     unsafe extern "C++" {
+        include!("wrapper/include/ro_mol.h");
         include!("wrapper/include/scaffold_network.h");
 
+        pub type ROMol = crate::ro_mol_ffi::ROMol;
         pub type ScaffoldNetworkParams;
         //
         // pub fn default_scaffold_network_params() -> SharedPtr<ScaffoldNetworkParams>;
@@ -39,5 +41,28 @@ pub mod ffi {
         pub type ScaffoldNetworkClass;
         pub fn default_scaffold_network() -> SharedPtr<ScaffoldNetworkClass>;
 
+        pub fn create_scaffold_network(
+            mol: &SharedPtr<ROMol>,
+            params: &SharedPtr<ScaffoldNetworkParams>,
+        ) -> Result<SharedPtr<ScaffoldNetworkClass>>;
+
+        pub fn update_scaffold_network(
+            network: &mut SharedPtr<ScaffoldNetworkClass>,
+            mol: &SharedPtr<ROMol>,
+            params: &SharedPtr<ScaffoldNetworkParams>,
+        ) -> Result<()>;
+
+        pub fn scaffold_network_nodes(network: &SharedPtr<ScaffoldNetworkClass>) -> Vec<String>;
+        pub fn scaffold_network_counts(network: &SharedPtr<ScaffoldNetworkClass>) -> Vec<u32>;
+        pub fn scaffold_network_mol_counts(network: &SharedPtr<ScaffoldNetworkClass>) -> Vec<u32>;
+        pub fn scaffold_network_edge_begin_indices(
+            network: &SharedPtr<ScaffoldNetworkClass>,
+        ) -> Vec<usize>;
+        pub fn scaffold_network_edge_end_indices(
+            network: &SharedPtr<ScaffoldNetworkClass>,
+        ) -> Vec<usize>;
+        pub fn scaffold_network_edge_types(
+            network: &SharedPtr<ScaffoldNetworkClass>,
+        ) -> Vec<String>;
     }
 }

rdkit-sys/tests/test_scaffold_network.rs

@@ -1,5 +1,23 @@
 use rdkit_sys::scaffold_network_ffi::*;
 
+/// 2-phenylquinoline: one fused bicycle plus a pendant ring, so the network has
+/// something to actually fragment.
+const PHENYLQUINOLINE: &str = "c1ccc(-c2ccc3ncccc3c2)cc1";
+
+fn mol(smiles: &str) -> cxx::SharedPtr<rdkit_sys::ro_mol_ffi::ROMol> {
+    cxx::let_cxx_string!(smiles_cxx_string = smiles);
+    rdkit_sys::ro_mol_ffi::smiles_to_mol(&smiles_cxx_string).unwrap()
+}
+
+/// The stock parameters minus the dummy-atom scaffolds, which keeps the
+/// assertions below readable.
+fn params_without_generics() -> cxx::SharedPtr<ScaffoldNetworkParams> {
+    let mut params = default_scaffold_network_params();
+    set_include_generic_scaffolds(&mut params, false);
+    include_generic_bond_scaffolds(&mut params, false);
+    params
+}
+
 #[test]
 fn test_scaffold_network() {
     default_scaffold_network_params();
@@ -18,3 +36,110 @@ fn test_scaffold_network() {
 
     let _scaffold_network = default_scaffold_network();
 }
+
+#[test]
+fn test_create_scaffold_network() {
+    let params = params_without_generics();
+    let network = create_scaffold_network(&mol(PHENYLQUINOLINE), &params).unwrap();
+
+    let nodes = scaffold_network_nodes(&network);
+    assert_eq!(
+        nodes,
+        [
+            PHENYLQUINOLINE,
+            "*c1ccccc1",
+            "c1ccccc1",
+            "*c1ccc2ncccc2c1",
+            "c1ccc2ncccc2c1",
+        ]
+    );
+
+    // one entry per node, so the downstream can zip them against `nodes`
+    assert_eq!(scaffold_network_counts(&network).len(), nodes.len());
+    assert_eq!(scaffold_network_mol_counts(&network).len(), nodes.len());
+    assert_eq!(scaffold_network_mol_counts(&network), [1, 1, 1, 1, 1]);
+
+    // edges run from the more specific scaffold to the more general one
+    let begin_indices = scaffold_network_edge_begin_indices(&network);
+    let end_indices = scaffold_network_edge_end_indices(&network);
+    let types = scaffold_network_edge_types(&network);
+    assert_eq!(begin_indices, [0, 1, 0, 3]);
+    assert_eq!(end_indices, [1, 2, 3, 4]);
+    assert_eq!(
+        types,
+        [
+            "Fragment",
+            "RemoveAttachment",
+            "Fragment",
+            "RemoveAttachment"
+        ]
+    );
+}
+
+#[test]
+fn test_create_scaffold_network_single_ring() {
+    let params = params_without_generics();
+    let network = create_scaffold_network(&mol("c1ccccc1"), &params).unwrap();
+
+    // nothing to strip: the molecule is already its own scaffold
+    assert_eq!(scaffold_network_nodes(&network), ["c1ccccc1"]);
+    assert!(scaffold_network_edge_begin_indices(&network).is_empty());
+}
+
+#[test]
+fn test_create_scaffold_network_no_rings() {
+    let params = params_without_generics();
+    let network = create_scaffold_network(&mol("CCO"), &params).unwrap();
+
+    // RDKit still emits a node for a ring-free molecule, it is just the empty
+    // SMILES
+    assert_eq!(scaffold_network_nodes(&network), ["CCO", ""]);
+    assert_eq!(scaffold_network_edge_types(&network), ["Initialize"]);
+}
+
+#[test]
+fn test_update_scaffold_network() {
+    let params = params_without_generics();
+    let mut network = create_scaffold_network(&mol("c1ccccc1"), &params).unwrap();
+    update_scaffold_network(&mut network, &mol("c1ccc2ncccc2c1"), &params).unwrap();
+
+    assert_eq!(
+        scaffold_network_nodes(&network),
+        ["c1ccccc1", "c1ccc2ncccc2c1"]
+    );
+    assert_eq!(scaffold_network_mol_counts(&network), [1, 1]);
+}
+
+/// `new_scaffold_network_params` hands its argument straight to the
+/// `ScaffoldNetworkParams(std::vector<std::string>)` constructor, so an empty
+/// vec leaves RDKit with no bond breaker reactions and nothing ever fragments.
+/// Callers who want the stock fragmentation want
+/// `default_scaffold_network_params`.
+#[test]
+fn test_empty_bond_breakers_do_not_fragment() {
+    let params = new_scaffold_network_params(&vec![]);
+    let network = create_scaffold_network(&mol(PHENYLQUINOLINE), &params).unwrap();
+    assert_eq!(scaffold_network_nodes(&network), [PHENYLQUINOLINE]);
+
+    let params =
+        new_scaffold_network_params(&vec!["[!#0;R:1]-!@[!#0:2]>>[*:1]-[#0].[#0]-[*:2]".into()]);
+    let network = create_scaffold_network(&mol(PHENYLQUINOLINE), &params).unwrap();
+    assert!(scaffold_network_nodes(&network).len() > 1);
+}
+
+/// RDKit throws when asked for neither flavor of scaffold; the bridge turns
+/// that into an `Err` rather than letting it unwind through Rust.
+#[test]
+fn test_create_scaffold_network_bad_params() {
+    let mut params = default_scaffold_network_params();
+    include_scaffolds_with_attachments(&mut params, false);
+    include_scaffolds_without_attachments(&mut params, false);
+
+    match create_scaffold_network(&mol(PHENYLQUINOLINE), &params) {
+        Err(e) => assert_eq!(
+            e.what(),
+            "must include at least one of scaffolds with attachments or scaffolds without attachments"
+        ),
+        Ok(_) => panic!("expected err variant"),
+    }
+}

rdkit-sys/wrapper/include/scaffold_network.h

@@ -24,4 +24,24 @@ void collect_mol_counts(std::shared_ptr<ScaffoldNetworkParams> &params, bool inp
 using ScaffoldNetworkClass = ScaffoldNetwork::ScaffoldNetwork;
 
 std::shared_ptr<ScaffoldNetworkClass> default_scaffold_network();
+
+// RDKit::ScaffoldNetwork::createScaffoldNetwork is a template over a container of molecules, but only
+// std::vector<ROMOL_SPTR> (boost::shared_ptr) and std::vector<std::shared_ptr<ROMol>> are explicitly
+// instantiated. cxx hands us one std::shared_ptr<ROMol> at a time, so we wrap it in a one element
+// vector here rather than fighting cxx over Vec<SharedPtr<ROMol>> (see BACKLOG.md). Callers who want
+// a network spanning several molecules build it up with update_scaffold_network.
+std::shared_ptr<ScaffoldNetworkClass> create_scaffold_network(const std::shared_ptr<ROMol> &mol,
+                                                              const std::shared_ptr<ScaffoldNetworkParams> &params);
+void update_scaffold_network(std::shared_ptr<ScaffoldNetworkClass> &network, const std::shared_ptr<ROMol> &mol,
+                             const std::shared_ptr<ScaffoldNetworkParams> &params);
+
+// The ScaffoldNetwork members come back as plain data, one parallel vector per member, so that no
+// C++ type crosses the bridge. edges[i] runs from edge_begin_indices[i] to edge_end_indices[i], both
+// of which index into the nodes vector.
+rust::Vec<rust::String> scaffold_network_nodes(const std::shared_ptr<ScaffoldNetworkClass> &network);
+rust::Vec<std::uint32_t> scaffold_network_counts(const std::shared_ptr<ScaffoldNetworkClass> &network);
+rust::Vec<std::uint32_t> scaffold_network_mol_counts(const std::shared_ptr<ScaffoldNetworkClass> &network);
+rust::Vec<std::size_t> scaffold_network_edge_begin_indices(const std::shared_ptr<ScaffoldNetworkClass> &network);
+rust::Vec<std::size_t> scaffold_network_edge_end_indices(const std::shared_ptr<ScaffoldNetworkClass> &network);
+rust::Vec<rust::String> scaffold_network_edge_types(const std::shared_ptr<ScaffoldNetworkClass> &network);
 } // namespace RDKit

rdkit-sys/wrapper/src/scaffold_network.cc

@@ -62,4 +62,79 @@ std::shared_ptr<ScaffoldNetworkClass> default_scaffold_network() {
 	return std::shared_ptr<ScaffoldNetworkClass>(scaffold_network);
 }
 
+std::shared_ptr<ScaffoldNetworkClass> create_scaffold_network(const std::shared_ptr<ROMol> &mol,
+                                                              const std::shared_ptr<ScaffoldNetworkParams> &params) {
+	std::vector<std::shared_ptr<ROMol>> mols{mol};
+	ScaffoldNetworkClass *scaffold_network =
+	    new ScaffoldNetworkClass(ScaffoldNetwork::createScaffoldNetwork(mols, *params));
+	return std::shared_ptr<ScaffoldNetworkClass>(scaffold_network);
+}
+
+void update_scaffold_network(std::shared_ptr<ScaffoldNetworkClass> &network, const std::shared_ptr<ROMol> &mol,
+                             const std::shared_ptr<ScaffoldNetworkParams> &params) {
+	std::vector<std::shared_ptr<ROMol>> mols{mol};
+	ScaffoldNetwork::updateScaffoldNetwork(mols, *network, *params);
+}
+
+rust::Vec<rust::String> scaffold_network_nodes(const std::shared_ptr<ScaffoldNetworkClass> &network) {
+	rust::Vec<rust::String> nodes;
+	nodes.reserve(network->nodes.size());
+
+	for (const auto &node : network->nodes) { nodes.push_back(rust::String(node)); }
+
+	return nodes;
+}
+
+rust::Vec<std::uint32_t> scaffold_network_counts(const std::shared_ptr<ScaffoldNetworkClass> &network) {
+	rust::Vec<std::uint32_t> counts;
+	counts.reserve(network->counts.size());
+
+	for (const auto count : network->counts) { counts.push_back(count); }
+
+	return counts;
+}
+
+rust::Vec<std::uint32_t> scaffold_network_mol_counts(const std::shared_ptr<ScaffoldNetworkClass> &network) {
+	rust::Vec<std::uint32_t> mol_counts;
+	mol_counts.reserve(network->molCounts.size());
+
+	for (const auto mol_count : network->molCounts) { mol_counts.push_back(mol_count); }
+
+	return mol_counts;
+}
+
+rust::Vec<std::size_t> scaffold_network_edge_begin_indices(const std::shared_ptr<ScaffoldNetworkClass> &network) {
+	rust::Vec<std::size_t> begin_indices;
+	begin_indices.reserve(network->edges.size());
+
+	for (const auto &edge : network->edges) { begin_indices.push_back(edge.beginIdx); }
+
+	return begin_indices;
+}
+
+rust::Vec<std::size_t> scaffold_network_edge_end_indices(const std::shared_ptr<ScaffoldNetworkClass> &network) {
+	rust::Vec<std::size_t> end_indices;
+	end_indices.reserve(network->edges.size());
+
+	for (const auto &edge : network->edges) { end_indices.push_back(edge.endIdx); }
+
+	return end_indices;
+}
+
+rust::Vec<rust::String> scaffold_network_edge_types(const std::shared_ptr<ScaffoldNetworkClass> &network) {
+	rust::Vec<rust::String> types;
+	types.reserve(network->edges.size());
+
+	// RDKit ships an operator<< for EdgeType; going through it keeps these strings in step with
+	// upstream rather than duplicating the enum names on our side.
+	std::ostringstream type_name;
+	for (const auto &edge : network->edges) {
+		type_name.str("");
+		type_name << edge.type;
+		types.push_back(rust::String(type_name.str()));
+	}
+
+	return types;
+}
+
 } // namespace RDKit

src/lib.rs

@@ -16,6 +16,9 @@ pub use mol_standardize::*;
 mod periodic_table;
 pub use periodic_table::*;
 
+mod scaffold_network;
+pub use scaffold_network::*;
+
 mod substruct_match;
 pub use substruct_match::*;