Re: [OpenBabel-Devel] Getting corresponding bonds and angles from two molecules

[Naruki Yoshikawa]

Hi Dave and everyone,

Thanks to your advice, the code below worked to get corresponding bonds.

  // mol1: ground truth loaded from SDF
  // mol2: predicted molecule structure loaded from SDF
  OBMol mol1, mol2;

  // Get mapping between two molecules
  OBQuery *query = CompileMoleculeQuery(&mol1);
  OBIsomorphismMapper *mapper = OBIsomorphismMapper::GetInstance(query);
  OBIsomorphismMapper::Mappings mappings;
  mapper->MapAll(&mol2, mappings);

  map<int, int> idxMapping; // mol1 index to mol2 index

  for(size_t m = 0; m < mappings.size(); ++m) {
    bool isCorrect = true;
    for(size_t i = 0; i < mappings[m].size(); ++i) {
      //obmol indices are 1-indexed while the mapper is zero indexed
      unsigned int idx1 = mappings[m][i].first + 1;
      unsigned int idx2 = mappings[m][i].second + 1;
      idxMapping[idx1] = idx2;

      if(mol1.GetAtom(idx1)->GetAtomicNum() !=
mol2.GetAtom(idx2)->GetAtomicNum()) {
        idxMapping.clear();
        isCorrect = false;
        break;
      }
    }
    if(isCorrect) break;
  }

  // Display the length of each bond
  FOR_BONDS_OF_MOL(bond1, mol1) {
    unsigned int begin1 = idxMapping[bond1->GetBeginAtomIdx()];
    unsigned int end1 = idxMapping[bond1->GetEndAtomIdx()];

    FOR_BONDS_OF_MOL(bond2, mol2) {
      unsigned int begin2 = bond2->GetBeginAtomIdx();
      unsigned int end2 = bond2->GetEndAtomIdx();
      if((begin1 == begin2 && end1 == end2)) {
        cout << bond1->GetLength() << "," << bond2->GetLength() << endl;
        break;
      }
    }
  }


I'm wondering why the following code using MapFirst to get idxMapping
did not work.
I assume that OBIsomorphism only considers connection, but is this
behavior intended?

    // Get mapping between two molecules
    OBQuery *query = CompileMoleculeQuery(&mol1);
    OBIsomorphismMapper *mapper = OBIsomorphismMapper::GetInstance(query);
    OBIsomorphismMapper::Mapping mapping;
    mapper->MapFirst(&mol2, mapping);

    map<int, int> idxMapping; // mol1 index to mol2 index

    for(size_t i = 0; i < mapping.size(); ++i) {
      //obmol indices are 1-indexed while the mapper is zero indexed
      unsigned int idx1 = mapping[i].first + 1;
      unsigned int idx2 = mapping[i].second + 1;
      idxMapping[idx1] = idx2;
      // this assertion fails in some molecules
      assert(mol1.GetAtom(idx1)->GetAtomicNum() !=
mol2.GetAtom(idx2)->GetAtomicNum());
    }

Thanks,
Naruki
2018年10月14日(日) 21:33 Koes, David <dk...@pitt.edu>:
>
> Try using OBIsomorphismMapper to find tge right mapping of atoms which will 
> implicitly give you the right matching for bonds.  For example:  
> https://github.com/dkoes/openbabel/blob/master/tools/obrms.cpp
>
> -Dave
>
> On Oct 14, 2018, at 4:01 AM, Naruki Yoshikawa <naruki.yoshik...@gmail.com> 
> wrote:
>
> Hi everyone,
>
> I'm evaluating my fragment-based coordinate generation method.
> I want to measure the error of bond lengths and angles, but I have a
> problem in getting corresponding ones from two molecules.
>
> What I want to do is like this:
>
>    // mol1: ground truth loaded from SDF
>    // mol2: predicted molecule structure loaded from SDF
>    OBMol mol1, mol2;
>
>    // Display the length of each bond
>    for (int idx = 0; idx < mol1.NumBonds(); ++idx) {
>      OBBond *bond1 = mol1.GetBond(idx);
>      OBBond *bond2 = mol2.GetBond(idx);
>
>      // these assertion fails
>      assert(bond1->GetBeginAtom()->GetAtomicNum()
>               == bond2->GetBeginAtom()->GetAtomicNum());
>      assert(bond1->GetEndAtom()->GetAtomicNum()
>               == bond2->GetEndAtom()->GetAtomicNum());
>
>      cout << bond1->GetLength() << "," << bond2->GetLength() << endl;
>    }
>
> The problem is that the bonds are not ordered properly.
> When I execute the code above, the assertions fail.
> So n-th bond in mol1 was different from n-th bond in mol2.
> I also have a similar problem for angles.
>
> If I can get something like "canonical order of bonds", I can compare
> the length of each bond correctly.
> However, I don't know how to do it.
>
> Does anyone have a good idea to get corresponding bonds and angles?
>
> Thanks,
> Naruki
>
>
> _______________________________________________
> OpenBabel-Devel mailing list
> OpenBabel-Devel@lists.sourceforge.net
> https://na01.safelinks.protection.outlook.com/?url=https%3A%2F%2Flists.sourceforge.net%2Flists%2Flistinfo%2Fopenbabel-devel&amp;data=02%7C01%7Cdkoes%40pitt.edu%7Ca4d3f190d0c74200797d08d631ab31a3%7C9ef9f489e0a04eeb87cc3a526112fd0d%7C1%7C0%7C636751008682354300&amp;sdata=6%2FgdvGOjbSNR5zoo20k4MPJVuZn3057ooWnoFpHo1rI%3D&amp;reserved=0


_______________________________________________
OpenBabel-Devel mailing list
OpenBabel-Devel@lists.sourceforge.net
https://lists.sourceforge.net/lists/listinfo/openbabel-devel