diff --git a/Code/Geometry/Transform3D.h b/Code/Geometry/Transform3D.h index 355ad41dc..4cccab6f8 100644 --- a/Code/Geometry/Transform3D.h +++ b/Code/Geometry/Transform3D.h @@ -85,7 +85,7 @@ class RDKIT_RDGEOMETRYLIB_EXPORT Transform3D * * The order is important here, on two transforms t1 and t2 * t3 = t1*t2 - * The resulting transform t3 has the folliwng effect + * The resulting transform t3 has the following effect * t3(point) = t1(t2(point)) */ RDKIT_RDGEOMETRYLIB_EXPORT RDGeom::Transform3D operator*( diff --git a/External/pubchem_shape/CMakeLists.txt b/External/pubchem_shape/CMakeLists.txt index 4de07a088..6b88bfd19 100644 --- a/External/pubchem_shape/CMakeLists.txt +++ b/External/pubchem_shape/CMakeLists.txt @@ -37,6 +37,7 @@ if(needDownload) endif() # simple patch for a typo in the pubchem_align library +# This should no longer be necessary, as it has been fixed upstream. Maybe take it out? file(READ ${PUBCHEMSHAPE_DIR}/shape_neighbor.cpp FILE_CONTENTS) string(REPLACE "memcpy( qom, old_quattrans, 7 * sizeof( float ) );" "memcpy( qom, old_quattrans, 7 * sizeof( double ) );" FILE_CONTENTS "${FILE_CONTENTS}") file(WRITE ${PUBCHEMSHAPE_DIR}/shape_neighbor.cpp "${FILE_CONTENTS}") @@ -48,7 +49,7 @@ if((MSVC AND RDK_INSTALL_DLLS_MSVC) OR ((NOT MSVC) AND WIN32)) endif() rdkit_library(PubChemShape PubChemShape.cpp SHARED - LINK_LIBRARIES pubchem_align3d SmilesParse SubstructMatch) + LINK_LIBRARIES pubchem_align3d SmilesParse SubstructMatch MolTransforms) target_compile_definitions(PubChemShape PRIVATE RDKIT_PUBCHEMSHAPE_BUILD) rdkit_headers(PubChemShape.hpp DEST GraphMol) diff --git a/External/pubchem_shape/PubChemShape.cpp b/External/pubchem_shape/PubChemShape.cpp index e5a970403..f36fa8773 100644 --- a/External/pubchem_shape/PubChemShape.cpp +++ b/External/pubchem_shape/PubChemShape.cpp @@ -3,7 +3,10 @@ #include #include +#include +#include #include +#include #include #include @@ -371,14 +374,17 @@ void extractFeatureCoords( } } -void extractCustomFeatureCoords(const unsigned int nSelectedAtoms, - const ShapeInputOptions &shapeOpts, - const RDGeom::Point3D &ave, - unsigned int &numFeatures, ShapeInput &res, - std::vector &rad_vector) { +void extractCustomFeatureCoords( + const unsigned int nSelectedAtoms, const ShapeInputOptions &shapeOpts, + const RDGeom::Point3D &ave, unsigned int &numFeatures, ShapeInput &res, + std::vector &rad_vector, + const std::unique_ptr &trans) { for (const auto &feature : shapeOpts.customFeatures) { unsigned int feature_type = std::get<0>(feature); RDGeom::Point3D floc = std::get<1>(feature); + if (trans) { + trans->TransformPoint(floc); + } double radius = std::get<2>(feature); floc -= ave; DEBUG_MSG("custom feature type " << feature_type << " (" << floc << ")"); @@ -394,7 +400,8 @@ void extractCustomFeatureCoords(const unsigned int nSelectedAtoms, } } // namespace -// The conformer is left where it is, the shape is translated to the origin. +// The conformer is left where it is, the shape is translated to the origin +// and aligned to the principal axes. ShapeInput PrepareConformer(const ROMol &mol, int confId, const ShapeInputOptions &shapeOpts) { Align3D::setUseCutOff(true); @@ -420,24 +427,61 @@ ShapeInput PrepareConformer(const ROMol &mol, int confId, std::vector rad_vector(nAlignmentAtoms); res.atom_type_vector.resize(nAlignmentAtoms, 0); + Conformer confCp(conformer); + std::unique_ptr trans; + RDGeom::Point3D initialCentroid; + if (shapeOpts.normalize) { + // This ignores hydrogens by default, which is what is needed here. + trans.reset(MolTransforms::computeCanonicalTransform(confCp)); + MolTransforms::transformConformer(confCp, *trans); + + // Save the rotation part of the transformation. The coords will be centred + // below and left that way so we don't need to keep the translation. + for (unsigned int i = 0, k = 0; i < 3; ++i) { + for (unsigned int j = 0; j < 3; ++j, ++k) { + res.inertialRot[k] = trans->getValUnchecked(i, j); + } + } + + // Work back from the canonical transformation to get the initial + // centroid by the inverse rotation of the current translation + // but not taking its negative so it can be used directly as the shift. + initialCentroid = RDGeom::Point3D{trans->getValUnchecked(0, 3), + trans->getValUnchecked(1, 3), + trans->getValUnchecked(2, 3)}; + RDGeom::Transform3D invTrans; + invTrans.setValUnchecked(0, 0, trans->getValUnchecked(0, 0)); + invTrans.setValUnchecked(0, 1, trans->getValUnchecked(1, 0)); + invTrans.setValUnchecked(0, 2, trans->getValUnchecked(2, 0)); + invTrans.setValUnchecked(1, 0, trans->getValUnchecked(0, 1)); + invTrans.setValUnchecked(1, 1, trans->getValUnchecked(1, 1)); + invTrans.setValUnchecked(1, 2, trans->getValUnchecked(2, 1)); + invTrans.setValUnchecked(2, 0, trans->getValUnchecked(0, 2)); + invTrans.setValUnchecked(2, 1, trans->getValUnchecked(1, 2)); + invTrans.setValUnchecked(2, 2, trans->getValUnchecked(2, 2)); + invTrans.TransformPoint(initialCentroid); + } + RWMol molCp(mol); + molCp.removeConformer(0); + molCp.addConformer(new Conformer(confCp), true); + RDGeom::Point3D ave; unsigned int nSelectedAtoms = 0; - extractAtomRadii(conformer, nAtoms, shapeOpts, ave, nSelectedAtoms, - rad_vector); - + // Calculates the average coord as by-product + extractAtomRadii(confCp, nAtoms, shapeOpts, ave, nSelectedAtoms, rad_vector); // translate steric center to origin DEBUG_MSG("steric center: (" << ave << ")"); res.shift = {-ave.x, -ave.y, -ave.z}; res.coord.resize(nAlignmentAtoms * 3); - extractAtomCoords(conformer, nAtoms, shapeOpts, ave, res.coord); + extractAtomCoords(confCp, nAtoms, shapeOpts, ave, res.coord); unsigned int numFeatures = 0; if (shapeOpts.customFeatures.empty()) { - extractFeatureCoords(conformer, nAtoms, nSelectedAtoms, feature_idx_type, + extractFeatureCoords(confCp, nAtoms, nSelectedAtoms, feature_idx_type, shapeOpts, ave, numFeatures, res, rad_vector); } else if (shapeOpts.useColors) { extractCustomFeatureCoords(nSelectedAtoms, shapeOpts, ave, numFeatures, res, - rad_vector); + rad_vector, trans); } // Now cut the final vectors down to the actual number of atoms and @@ -467,6 +511,9 @@ ShapeInput PrepareConformer(const ROMol &mol, int confId, res.coord.data(), res.alpha_vector, res.colorAtomType2IndexVectorMap); DEBUG_MSG("sof: " << res.sof); } + + // Finally set the shift to the negative of the original average coords. + res.shift = {initialCentroid.x, initialCentroid.y, initialCentroid.z}; return res; } @@ -499,7 +546,6 @@ std::pair AlignShape(const ShapeInput &refShape, fitShape.volumeAtomIndexVector, fitMapCp, fitShape.sov, fitShape.sof, !jointColorAtomTypeSet.empty(), max_preiters, max_postiters, opt_param, matrix.data(), nbr_st, nbr_ct); - DEBUG_MSG("Done!"); DEBUG_MSG("nbr_st: " << nbr_st); DEBUG_MSG("nbr_ct: " << nbr_ct); @@ -512,36 +558,73 @@ std::pair AlignShape(const ShapeInput &refShape, } void TransformConformer(const std::vector &finalTrans, - const std::vector &matrix, ShapeInput &fitShape, - Conformer &fitConf) { - // we reuse/modify the coord member of fitShape in order to avoid memory - // allocations - const unsigned int nAtoms = fitConf.getOwningMol().getNumAtoms(); - if (nAtoms > fitShape.volumeAtomIndexVector.size()) { - // Hs are missing... make sure we have space for them - fitShape.coord.resize(3 * nAtoms); - } - // initialize the fitShape coords with the starting atomic positions from - // the conformer shifted to the center of "mass" coordinates. - for (unsigned int i = 0; i < nAtoms; ++i) { - const auto &pos = fitConf.getAtomPos(i); - fitShape.coord[i * 3] = pos.x + fitShape.shift[0]; - fitShape.coord[i * 3 + 1] = pos.y + fitShape.shift[1]; - fitShape.coord[i * 3 + 2] = pos.z + fitShape.shift[2]; - } + const std::vector &finalRot, + const std::vector &matrix, + const ShapeInput &fitShape, Conformer &fitConf) { + // Move fitConf to fitShape's initial centroid and principal axes + RDGeom::Transform3D transform0; + transform0.SetTranslation( + RDGeom::Point3D{fitShape.shift[0], fitShape.shift[1], fitShape.shift[2]}); - std::vector transformed(nAtoms * 3); - Align3D::VApplyRotTransMatrix(transformed.data(), fitShape.coord.data(), - nAtoms, matrix.data()); + RDGeom::Transform3D transform1; + transform1.setValUnchecked(0, 0, fitShape.inertialRot[0]); + transform1.setValUnchecked(0, 1, fitShape.inertialRot[1]); + transform1.setValUnchecked(0, 2, fitShape.inertialRot[2]); + transform1.setValUnchecked(1, 0, fitShape.inertialRot[3]); + transform1.setValUnchecked(1, 1, fitShape.inertialRot[4]); + transform1.setValUnchecked(1, 2, fitShape.inertialRot[5]); + transform1.setValUnchecked(2, 0, fitShape.inertialRot[6]); + transform1.setValUnchecked(2, 1, fitShape.inertialRot[7]); + transform1.setValUnchecked(2, 2, fitShape.inertialRot[8]); + transform1.setValUnchecked(2, 3, 0.0); + transform1.setValUnchecked(3, 0, 0.0); + transform1.setValUnchecked(3, 1, 0.0); + transform1.setValUnchecked(3, 2, 0.0); + transform1.setValUnchecked(3, 3, 1.0); - // now set the coordinates in the conformer; undo the shift of the reference - // shape to center of "mass" coordinates - for (unsigned i = 0; i < nAtoms; ++i) { - RDGeom::Point3D &pos = fitConf.getAtomPos(i); - pos.x = transformed[i * 3] - finalTrans[0]; - pos.y = transformed[(i * 3) + 1] - finalTrans[1]; - pos.z = transformed[(i * 3) + 2] - finalTrans[2]; - } + // Apply the overlay matrix. + // The matrix comes in as a linear form of a 3x3 rotation matrix and a + // 3x1 translation. Convert that into a RDGeom::Trans3D + RDGeom::Transform3D transform2; + transform2.setValUnchecked(0, 0, matrix[0]); + transform2.setValUnchecked(0, 1, matrix[1]); + transform2.setValUnchecked(0, 2, matrix[2]); + transform2.setValUnchecked(0, 3, matrix[9]); + transform2.setValUnchecked(1, 0, matrix[3]); + transform2.setValUnchecked(1, 1, matrix[4]); + transform2.setValUnchecked(1, 2, matrix[5]); + transform2.setValUnchecked(1, 3, matrix[10]); + transform2.setValUnchecked(2, 0, matrix[6]); + transform2.setValUnchecked(2, 1, matrix[7]); + transform2.setValUnchecked(2, 2, matrix[8]); + transform2.setValUnchecked(2, 3, matrix[11]); + transform2.setValUnchecked(3, 0, 0.0); + transform2.setValUnchecked(3, 1, 0.0); + transform2.setValUnchecked(3, 2, 0.0); + transform2.setValUnchecked(3, 3, 1.0); + + // Now apply the final rotation and final translation, which + // is to put it into the reference shape's frame of reference. + RDGeom::Transform3D transform3; + transform3.setValUnchecked(0, 0, finalRot[0]); + transform3.setValUnchecked(0, 1, finalRot[1]); + transform3.setValUnchecked(0, 2, finalRot[2]); + transform3.setValUnchecked(0, 3, finalTrans[0]); + transform3.setValUnchecked(1, 0, finalRot[3]); + transform3.setValUnchecked(1, 1, finalRot[4]); + transform3.setValUnchecked(1, 2, finalRot[5]); + transform3.setValUnchecked(1, 3, finalTrans[1]); + transform3.setValUnchecked(2, 0, finalRot[6]); + transform3.setValUnchecked(2, 1, finalRot[7]); + transform3.setValUnchecked(2, 2, finalRot[8]); + transform3.setValUnchecked(2, 3, finalTrans[2]); + transform3.setValUnchecked(3, 0, 0.0); + transform3.setValUnchecked(3, 1, 0.0); + transform3.setValUnchecked(3, 2, 0.0); + transform3.setValUnchecked(3, 3, 1.0); + + auto finalTransform = transform3 * transform2 * transform1 * transform0; + MolTransforms::transformConformer(fitConf, finalTransform); } std::pair AlignMolecule( @@ -555,14 +638,23 @@ std::pair AlignMolecule( auto fitShape = PrepareConformer(fit, fitConfId, shapeOpts); auto tanis = AlignShape(refShape, fitShape, matrix, opt_param, max_preiters, max_postiters); - // transform fit coords Conformer &fit_conformer = fit.getConformer(fitConfId); std::vector finalTrans{0.0, 0.0, 0.0}; + std::vector finalRot{1.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 1.0}; if (applyRefShift) { - finalTrans = refShape.shift; + finalTrans = std::vector{-refShape.shift[0], -refShape.shift[1], + -refShape.shift[2]}; + // set the final rotation to the inverse of the original inertialRot + // of the reference. + finalRot = + std::vector{refShape.inertialRot[0], refShape.inertialRot[3], + refShape.inertialRot[6], refShape.inertialRot[1], + refShape.inertialRot[4], refShape.inertialRot[7], + refShape.inertialRot[2], refShape.inertialRot[5], + refShape.inertialRot[8]}; } - TransformConformer(finalTrans, matrix, fitShape, fit_conformer); + TransformConformer(finalTrans, finalRot, matrix, fitShape, fit_conformer); fit.setProp("shape_align_shape_tanimoto", tanis.first); fit.setProp("shape_align_color_tanimoto", tanis.second); @@ -582,11 +674,11 @@ std::pair AlignMolecule( std::pair AlignMolecule( const ROMol &ref, ROMol &fit, std::vector &matrix, const ShapeInputOptions &refShapeOpts, - const ShapeInputOptions &probeShapeOpts, int refConfId, int fitConfId, + const ShapeInputOptions &fitShapeOpts, int refConfId, int fitConfId, double opt_param, unsigned int max_preiters, unsigned int max_postiters) { Align3D::setUseCutOff(true); - if (refShapeOpts.useColors != probeShapeOpts.useColors) { + if (refShapeOpts.useColors != fitShapeOpts.useColors) { BOOST_LOG(rdWarningLog) << "useColor values inconsistent between the reference and fit molecules. Colors will not be used in the alignment." << std::endl; @@ -596,7 +688,7 @@ std::pair AlignMolecule( auto refShape = PrepareConformer(ref, refConfId, refShapeOpts); bool applyRefShift = true; auto scores = - AlignMolecule(refShape, fit, matrix, probeShapeOpts, fitConfId, opt_param, + AlignMolecule(refShape, fit, matrix, fitShapeOpts, fitConfId, opt_param, max_preiters, max_postiters, applyRefShift); return scores; } diff --git a/External/pubchem_shape/PubChemShape.hpp b/External/pubchem_shape/PubChemShape.hpp index 6748e8366..ef41d6f8b 100644 --- a/External/pubchem_shape/PubChemShape.hpp +++ b/External/pubchem_shape/PubChemShape.hpp @@ -1,6 +1,8 @@ #ifndef RDKIT_PUBCHEMSHAPE_GUARD #define RDKIT_PUBCHEMSHAPE_GUARD +#include "Geometry/Transform3D.h" + #include #include #include @@ -54,6 +56,7 @@ struct RDKIT_PUBCHEMSHAPE_EXPORT ShapeInput { ar & shift; ar & sov; ar & sof; + ar & inertialRot; } #endif @@ -64,6 +67,10 @@ struct RDKIT_PUBCHEMSHAPE_EXPORT ShapeInput { std::map> colorAtomType2IndexVectorMap; std::vector shift; + // If the conformer the shape was created from was rotated into the + // inertial reference frame at the start, this is the rotation that + // did that, assuming it was already centred on the origin. + std::vector inertialRot{1.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 1.0}; double sov{0.0}; double sof{0.0}; }; @@ -98,7 +105,8 @@ struct RDKIT_PUBCHEMSHAPE_EXPORT ShapeInputOptions { \param confId (optional) the conformer to use \param shapeOpts (optional) Change the default behaviour. - \return a ShapeInput object, translated to the origin + \return a ShapeInput object, translated to the origin and aligned along its + principal axes. */ RDKIT_PUBCHEMSHAPE_EXPORT ShapeInput PrepareConformer(const RDKit::ROMol &mol, int confId = -1, @@ -123,18 +131,20 @@ RDKIT_PUBCHEMSHAPE_EXPORT std::pair AlignShape( unsigned int max_preiters = 10u, unsigned int max_postiters = 30u); //! Assuming that fitShape has been overlaid onto a reference shape to give -//! the ! transformation matrix, apply the same transformation to the given +//! the transformation matrix, apply the same transformation to the given //! conformer. /*! \param finalTrans the final translation to apply to the fitConf coords. + \param finalRot the final rotation to apply to the fitConf coords. \param matrix the transformation matrix produced from alignment \param fitShape the shape that was aligned. The coord vector of this is modified \param fitConf the conformation to be transformed */ RDKIT_PUBCHEMSHAPE_EXPORT void TransformConformer( - const std::vector &finalTrans, const std::vector &matrix, - ShapeInput &fitShape, RDKit::Conformer &fitConf); + const std::vector &finalTrans, const std::vector &finalRot, + const std::vector &matrix, const ShapeInput &fitShape, + RDKit::Conformer &fitConf); //! Align a molecule to a reference shape /*! @@ -202,7 +212,7 @@ RDKIT_PUBCHEMSHAPE_EXPORT std::pair AlignMolecule( RDKIT_PUBCHEMSHAPE_EXPORT std::pair AlignMolecule( const RDKit::ROMol &ref, RDKit::ROMol &fit, std::vector &matrix, const ShapeInputOptions &refShapeOpts, - const ShapeInputOptions &probeShapeOpts, int refConfId = -1, + const ShapeInputOptions &fitShapeOpts, int refConfId = -1, int fitConfId = -1, double opt_param = 1.0, unsigned int max_preiters = 10u, unsigned int max_postiters = 30u); diff --git a/External/pubchem_shape/Wrap/cshapealign.cpp b/External/pubchem_shape/Wrap/cshapealign.cpp index 6b2b72ff2..dd1a40486 100644 --- a/External/pubchem_shape/Wrap/cshapealign.cpp +++ b/External/pubchem_shape/Wrap/cshapealign.cpp @@ -96,6 +96,7 @@ python::tuple scoreShape(const ShapeInput &shape1, ShapeInput &shape2, return python::make_tuple(nbr_st, nbr_ct); } void transformConformer(const python::list &pyFinalTrans, + const python::list &pyFinalRot, const python::list &pyMatrix, ShapeInput probeShape, RDKit::Conformer &probeConf) { std::vector matrix; @@ -112,7 +113,14 @@ void transformConformer(const python::list &pyFinalTrans, "The final translation vector must have 3 values. It had " + std::to_string(finalTrans.size()) + "."); } - TransformConformer(finalTrans, matrix, probeShape, probeConf); + std::vector finalRot; + pythonObjectToVect(pyFinalRot, finalRot); + if (finalRot.size() != 9) { + throw_value_error("The final rotation vector must have 9 values. It had " + + std::to_string(finalRot.size()) + "."); + } + + TransformConformer(finalTrans, finalRot, matrix, probeShape, probeConf); } ShapeInput *prepConf(const RDKit::ROMol &mol, int confId, const python::object &py_opts) { @@ -175,6 +183,14 @@ python::list get_shapeShift(const ShapeInput &shp) { return py_list; } +python::list get_inertialRot(const ShapeInput &shp) { + python::list py_list; + for (const auto &val : shp.inertialRot) { + py_list.append(val); + } + return py_list; +} + void set_customFeatures(ShapeInputOptions &shp, const python::object &s) { shp.customFeatures.clear(); auto len = python::len(s); @@ -428,6 +444,9 @@ Returns &ShapeInput::volumeAtomIndexVector) .add_property("shift", &helpers::get_shapeShift, &helpers::set_shapeShift, "Translation of centre of shape coordinates to origin.") + .add_property( + "inertialRot", &helpers::get_inertialRot, + "Rotation applied to put the shape into its principal axes frame of reference.") .def_readwrite("sov", &ShapeInput::sov) .def_readwrite("sof", &ShapeInput::sof); diff --git a/External/pubchem_shape/Wrap/test_rdshapealign.py b/External/pubchem_shape/Wrap/test_rdshapealign.py index 530c57989..d79feda8b 100644 --- a/External/pubchem_shape/Wrap/test_rdshapealign.py +++ b/External/pubchem_shape/Wrap/test_rdshapealign.py @@ -18,7 +18,7 @@ class TestCase(unittest.TestCase): tpl = rdShapeAlign.AlignMol(self.ref, self.probe, opt_param=0.5, max_preiters=3, max_postiters=16) self.assertAlmostEqual(tpl[0], 0.773, places=3) - self.assertAlmostEqual(tpl[1], 0.303, places=3) + self.assertAlmostEqual(tpl[1], 0.305, places=3) def test2_NoColor(self): tpl = rdShapeAlign.AlignMol(self.ref, self.probe, useColors=False, opt_param=0.5, @@ -31,7 +31,7 @@ class TestCase(unittest.TestCase): self.assertTrue(type(shp) == rdShapeAlign.ShapeInput) tpl = rdShapeAlign.AlignMol(shp, self.probe, opt_param=0.5, max_preiters=3, max_postiters=16) self.assertAlmostEqual(tpl[0], 0.773, places=3) - self.assertAlmostEqual(tpl[1], 0.303, places=3) + self.assertAlmostEqual(tpl[1], 0.305, places=3) def test4_ShapeInputOptions(self): opts = rdShapeAlign.ShapeInputOptions() @@ -58,13 +58,17 @@ class TestCase(unittest.TestCase): probeShp = rdShapeAlign.PrepareConformer(self.probe, -1) tpl = rdShapeAlign.AlignShapes(refShp, probeShp) probeCp = Chem.Mol(self.probe) - rdShapeAlign.TransformConformer(refShp.shift, tpl[2], probeShp, probeCp.GetConformer(-1)) + rdShapeAlign.TransformConformer(refShp.shift, refShp.inertialRot, tpl[2], probeShp, probeCp.GetConformer(-1)) + self.assertEqual(len(refShp.inertialRot), 9) + self.assertAlmostEqual(refShp.inertialRot[0], -0.926125, places=6) + self.assertAlmostEqual(refShp.inertialRot[8], -0.852890, places=6) + # Just show it did something. The full test is in the C++ layer. self.assertNotEqual(self.probe.GetConformer().GetAtomPosition(0), probeCp.GetConformer().GetAtomPosition(0)) matrix = tpl[2][:10] with self.assertRaises(ValueError): - rdShapeAlign.TransformConformer(refShp.shift, matrix, probeShp, probeCp.GetConformer(-1)) + rdShapeAlign.TransformConformer(refShp.shift, refShp.inertialRot, matrix, probeShp, probeCp.GetConformer(-1)) def test6_notColorAtoms(self): m1 = Chem.MolFromSmiles("Nc1ccccc1 |(0.392086,-2.22477,0.190651;" @@ -91,8 +95,8 @@ class TestCase(unittest.TestCase): 0), 1.0), (1, Point3D(1.7571, -0.120174, 0.1), 1.0)) shp2 = rdShapeAlign.PrepareConformer(m2, -1, opts2) tpl = rdShapeAlign.AlignShapes(shp, shp2, opt_param=0.5) - self.assertAlmostEqual(tpl[0], 0.997, places=3) - self.assertAlmostEqual(tpl[1], 0.978, places=3) + self.assertAlmostEqual(tpl[0], 1.000, places=3) + self.assertAlmostEqual(tpl[1], 0.997, places=3) tf = tpl[2] self.assertGreater(0.0, tf[0]) self.assertLess(0.0, tf[3 * 3]) @@ -114,8 +118,8 @@ class TestCase(unittest.TestCase): opts2.customFeatures = ((2, Point3D(-1.75978, 0.148897, 0), 1.0), (1, Point3D(1.7571, -0.120174, 0.1), 1.0)) tpl = rdShapeAlign.AlignMol(m1, m2, opts, opts2, opt_param=0.5) - self.assertAlmostEqual(tpl[0], 0.997, places=3) - self.assertAlmostEqual(tpl[1], 0.978, places=3) + self.assertAlmostEqual(tpl[0], 1.000, places=3) + self.assertAlmostEqual(tpl[1], 0.997, places=3) def test9_FixedScore(self): # Just to make sure it's there and returns a value. diff --git a/External/pubchem_shape/test.cpp b/External/pubchem_shape/test.cpp index 574ef6e80..08d639496 100644 --- a/External/pubchem_shape/test.cpp +++ b/External/pubchem_shape/test.cpp @@ -1,16 +1,19 @@ #include #include +#include #include #include #include #include +#include #include #include #include #include "PubChemShape.hpp" +#include "GraphMol/SmilesParse/SmilesWrite.h" using namespace RDKit; @@ -83,11 +86,16 @@ TEST_CASE("basic alignment") { test_opt_param, test_max_preiters, test_max_postiters); CHECK_THAT(nbr_st, Catch::Matchers::WithinAbs(nbr_st2, 0.005)); CHECK_THAT(nbr_ct, Catch::Matchers::WithinAbs(nbr_ct2, 0.005)); - for (auto i = 0u; i < probe->getNumAtoms(); ++i) { - CHECK_THAT(cp.getConformer().getAtomPos(i).x, - Catch::Matchers::WithinAbs(cp2.getConformer().getAtomPos(i).x, - 0.005)); + CHECK_THAT( + cp.getConformer().getAtomPos(i).x, + Catch::Matchers::WithinAbs(cp2.getConformer().getAtomPos(i).x, 0.05)); + CHECK_THAT( + cp.getConformer().getAtomPos(i).y, + Catch::Matchers::WithinAbs(cp2.getConformer().getAtomPos(i).y, 0.05)); + CHECK_THAT( + cp.getConformer().getAtomPos(i).z, + Catch::Matchers::WithinAbs(cp2.getConformer().getAtomPos(i).z, 0.05)); } } } @@ -126,14 +134,15 @@ TEST_CASE("handling molecules with Hs") { REQUIRE(ref); auto probe = suppl[1]; REQUIRE(probe); + ShapeInputOptions mol1Opts, mol2Opts; SECTION("basics") { RWMol cp(*probe); std::vector matrix(12, 0.0); auto [nbr_st, nbr_ct] = AlignMolecule(*ref, cp, matrix, -1, -1, test_use_colors, test_opt_param, test_max_preiters, test_max_postiters); - CHECK_THAT(nbr_st, Catch::Matchers::WithinAbs(0.837, 0.005)); - CHECK_THAT(nbr_ct, Catch::Matchers::WithinAbs(0.694, 0.005)); + CHECK_THAT(nbr_st, Catch::Matchers::WithinAbs(0.840, 0.005)); + CHECK_THAT(nbr_ct, Catch::Matchers::WithinAbs(0.753, 0.005)); for (auto i = 0u; i < cp.getNumAtoms(); ++i) { // the failure mode here was that Hs had HUGE coordinates auto pos = cp.getConformer().getAtomPos(i); @@ -253,17 +262,19 @@ $$$$ auto [nbr_st, nbr_ct] = AlignMolecule(*ref, cp, matrix, refConfId, prbConfId, useColors, opt_param, preiters, postiters); - CHECK_THAT(nbr_st, Catch::Matchers::WithinAbs(0.501, 0.005)); - CHECK_THAT(nbr_ct, Catch::Matchers::WithinAbs(0.107, 0.005)); + + CHECK_THAT(nbr_st, Catch::Matchers::WithinAbs(0.483, 0.005)); + CHECK_THAT(nbr_ct, Catch::Matchers::WithinAbs(0.046, 0.005)); RWMol cp2(cp); auto [nbr_st2, nbr_ct2] = AlignMolecule(*ref, cp2, matrix, refConfId, prbConfId, useColors, opt_param, preiters, postiters); - CHECK_THAT(nbr_st2, Catch::Matchers::WithinAbs(nbr_st2, 0.005)); - CHECK_THAT(nbr_ct2, Catch::Matchers::WithinAbs(nbr_ct2, 0.005)); + + CHECK_THAT(nbr_st2, Catch::Matchers::WithinAbs(nbr_st, 0.005)); + CHECK_THAT(nbr_ct2, Catch::Matchers::WithinAbs(nbr_ct, 0.005)); auto rmsd = MolAlign::CalcRMS(cp, cp2); - CHECK_THAT(rmsd, Catch::Matchers::WithinAbs(0.017, 0.005)); + CHECK_THAT(rmsd, Catch::Matchers::WithinAbs(0.007, 0.005)); } } @@ -437,36 +448,14 @@ TEST_CASE("Shape-Shape alignment") { std::vector matrix(12, 0.0); auto [mol_st, mol_ct] = AlignMolecule(*ref, *probe, matrix, -1, -1); - auto [shape_st, shape_ct] = AlignShape(refShape, probeShape, matrix); // Check that the same results are achieved when overlaying the probe - // molecule onto the reference and the probe shape onto the reference shape - + // molecule onto the reference molecule and the probe shape onto the reference + // shape. The shapes are transformed to their centroid and principal axes, so + // the final coords will not match. Just check the tanimotos. CHECK_THAT(shape_st, Catch::Matchers::WithinAbs(mol_st, 0.001)); CHECK_THAT(shape_ct, Catch::Matchers::WithinAbs(mol_ct, 0.001)); - for (unsigned int i = 0; i < probe->getNumAtoms(); i++) { - const auto &pos = probe->getConformer().getAtomPos(i); - CHECK_THAT(pos.x, Catch::Matchers::WithinAbs( - probeShape.coord[3 * i] - refShape.shift[0], 0.001)); - CHECK_THAT(pos.y, - Catch::Matchers::WithinAbs( - probeShape.coord[3 * i + 1] - refShape.shift[1], 0.001)); - CHECK_THAT(pos.z, - Catch::Matchers::WithinAbs( - probeShape.coord[3 * i + 2] - refShape.shift[2], 0.001)); - } - - // Also check the TransformConformer function - TransformConformer(refShape.shift, matrix, probeShape, - probeCP.getConformer(-1)); - for (unsigned int i = 0; i < probe->getNumAtoms(); i++) { - const auto &pos = probe->getConformer().getAtomPos(i); - const auto &poscp = probeCP.getConformer().getAtomPos(i); - CHECK_THAT(pos.x, Catch::Matchers::WithinAbs(poscp.x, 0.001)); - CHECK_THAT(pos.y, Catch::Matchers::WithinAbs(poscp.y, 0.001)); - CHECK_THAT(pos.z, Catch::Matchers::WithinAbs(poscp.z, 0.001)); - } } TEST_CASE("Atoms excluded from Color features") { @@ -559,13 +548,13 @@ TEST_CASE("custom feature points") { double opt_param = 0.5; std::vector matrix(12, 0.0); auto [st, ct] = AlignShape(shape2, shape3, matrix, opt_param); - CHECK_THAT(st, Catch::Matchers::WithinAbs(0.997, 0.001)); - CHECK_THAT(ct, Catch::Matchers::WithinAbs(0.978, 0.001)); + CHECK_THAT(st, Catch::Matchers::WithinAbs(1.000, 0.001)); + CHECK_THAT(ct, Catch::Matchers::WithinAbs(0.997, 0.001)); CHECK(shape3.coord[0] > 0); // x coord of first atom CHECK(shape3.coord[3 * 3] < 0); // x coord of fourth atom auto conf = m2.getConformer(-1); - TransformConformer(shape2.shift, matrix, shape3, conf); + TransformConformer(shape2.shift, shape2.inertialRot, matrix, shape3, conf); CHECK(conf.getAtomPos(0).x > 0); CHECK(conf.getAtomPos(3).x < 0); } @@ -585,8 +574,8 @@ TEST_CASE("custom feature points") { std::vector matrix(12, 0.0); auto [st, ct] = AlignMolecule(*m1, m2, matrix, opts2, opts3, -2, -2, opt_param); - CHECK_THAT(st, Catch::Matchers::WithinAbs(0.997, 0.001)); - CHECK_THAT(ct, Catch::Matchers::WithinAbs(0.978, 0.001)); + CHECK_THAT(st, Catch::Matchers::WithinAbs(1.000, 0.001)); + CHECK_THAT(ct, Catch::Matchers::WithinAbs(0.997, 0.001)); auto conf = m2.getConformer(-1); CHECK(conf.getAtomPos(0).x > 0); CHECK(conf.getAtomPos(3).x < 0); @@ -656,4 +645,21 @@ TEST_CASE("Score No Overlay") { CHECK_THAT(singleShape, Catch::Matchers::WithinAbs(0.3376, 0.0005)); CHECK_THAT(singleColor, Catch::Matchers::WithinAbs(0.2674, 0.0005)); } -} \ No newline at end of file +} + +TEST_CASE("Iressa onto Tagrisso") { + // Conformations from PubChem produced by Omega. Iressa rotated and translated + // by a random amount. PubChem puts them both in their inertial frame which + // makes things too easy. + auto tagrisso = + "C=CC(=O)Nc1cc(Nc2nccc(-c3cn(C)c4ccccc34)n2)c(OC)cc1N(C)CCN(C)C |(-0.9161,3.8415,-2.9811;0.1848,3.1933,-2.588;0.1064,1.7789,-2.12;-0.9619,1.1797,-2.0847;1.3654,1.2872,-1.7553;1.6841,0.0144,-1.273;0.6638,-0.9235,-1.1146;0.9578,-2.1997,-0.6343;-0.0813,-3.1358,-0.4783;-1.4556,-2.9979,-0.1847;-2.1716,-4.1359,-0.1085;-3.4803,-3.9673,0.173;-4.0689,-2.7353,0.3728;-3.2269,-1.647,0.2676;-3.7311,-0.317,0.4568;-5.0275,0.0291,0.153;-5.1887,1.3569,0.4454;-6.4231,2.0889,0.2595;-4.0141,1.8811,0.9361;-3.7121,3.1796,1.3615;-2.4139,3.4249,1.8179;-1.4588,2.4106,1.8467;-1.7752,1.1164,1.4181;-3.0776,0.8453,0.9537;-1.9103,-1.7423,-0.011;2.2723,-2.5382,-0.3127;2.58,-3.7798,0.1575;2.539,-3.9651,1.571;3.2927,-1.6003,-0.4713;2.9986,-0.324,-0.9514;4.0475,0.61,-1.1047;4.7738,0.6956,-2.3546;4.4021,1.497,-0.0162;5.4401,0.8254,0.8736;5.8294,1.7155,1.9601;4.8213,1.7057,3.0218;7.1361,1.3324,2.4981)|"_smiles; + REQUIRE(tagrisso); + auto iressa = + "COc1cc2ncnc(Nc3ccc(F)c(Cl)c3)c2cc1OCCCN1CCOCC1 |(11.4672,-0.467948,5.63989;12.0133,0.532631,6.49693;11.2039,1.5801,6.81985;11.2014,2.71958,6.00975;10.3926,3.81652,6.29699;10.4038,4.90395,5.50623;9.58889,5.91871,5.85946;8.76443,5.96486,6.91838;8.77814,4.86059,7.68868;7.92337,4.86224,8.81914;7.44878,5.8925,9.64622;8.22182,7.03851,9.85619;7.75051,8.06265,10.6777;6.50441,7.94546,11.2936;6.06567,8.93802,12.0809;5.72932,6.80403,11.0875;4.19056,6.65372,11.8447;6.20047,5.78015,10.2656;9.57161,3.74547,7.43436;9.56851,2.60328,8.25407;10.3868,1.52933,7.93769;10.3797,0.419365,8.74203;11.3064,0.402096,9.81907;10.7104,-0.399165,10.9685;9.40938,0.22121,11.4678;9.64205,1.59049,11.9223;8.38006,2.22985,12.3199;8.64991,3.65266,12.8011;9.56883,3.64192,13.8942;10.8103,3.05101,13.5078;10.5931,1.61394,13.0425)|"_smiles; + REQUIRE(iressa); + std::vector matrix(12, 0.0); + auto sims = + AlignMolecule(*tagrisso, *iressa, matrix, -1, -1, true, 0.5, 10, 30); + CHECK_THAT(sims.first, Catch::Matchers::WithinAbs(0.582, 0.005)); + CHECK_THAT(sims.second, Catch::Matchers::WithinAbs(0.092, 0.005)); +}