Bioisostere Examples

The following table lists the currently available Bioisostere TK examples:

Program

Description

database

Building Brood Database (CHOMP)

query

Creating Brood Query

hitlist

Generating Brood Hits

matches

Generating Brood Matches

fragment_overlay

Overlay between Fragments

replace_fragment

Replacing a Fragment in a Molecule

cpddb

Building and Using an External Compound Database

2molLinker_query

Creating a Brood Query by Linking Two Molecules (Bridging)

Building Brood Database (CHOMP)

The following code example shows how to build a Brood database from a library of molecules.

Listing 1: Building Brood Database (CHOMP)

/*
 (C) 2022 Cadence Design Systems, Inc. (Cadence)
 All rights reserved.
 TERMS FOR USE OF SAMPLE CODE The software below ("Sample Code") is
 provided to current licensees or subscribers of Cadence products or
 SaaS offerings (each a "Customer").
 Customer is hereby permitted to use, copy, and modify the Sample Code,
 subject to these terms. Cadence claims no rights to Customer's
 modifications. Modification of Sample Code is at Customer's sole and
 exclusive risk. Sample Code may require Customer to have a then
 current license or subscription to the applicable Cadence offering.
 THE SAMPLE CODE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
 EXPRESS OR IMPLIED.  OPENEYE DISCLAIMS ALL WARRANTIES, INCLUDING, BUT
 NOT LIMITED TO, WARRANTIES OF MERCHANTABILITY, FITNESS FOR A
 PARTICULAR PURPOSE AND NONINFRINGEMENT. In no event shall Cadence be
 liable for any damages or liability in connection with the Sample Code
 or its use.
 */
#include <cstdio>
#include <string>
#include <openeye.h>

#include <oeplatform.h>
#include <oesystem.h>
#include <oechem.h>
#include <oebioisostere.h>
#include <oemolprop.h>


class ChompOptions : public OESystem::OEOptions
{
public:
  ChompOptions(std::string name = "ChompOptions")
    :OESystem::OEOptions(name)
  {
    OESystem::OEStringParameter pDBName("-out");
    pDBName.SetRequired(true);
    pDBName.SetKeyless(2);
    pDBName.SetVisibility(OESystem::OEParamVisibility::Simple);
    pDBName.SetBrief("Output database folder name");
    m_dbParam = AddParameter(pDBName);
    
    m_FragmentOptions = static_cast<OEBioisostere::OEFragmentOptions*>(AddOption(OEBioisostere::OEFragmentOptions()));
    m_ScreenOptions = static_cast<OEBioisostere::OEDBScreenOptions*>(AddOption(OEBioisostere::OEDBScreenOptions()));
  }
  
  ChompOptions(const ChompOptions&) = default;
  ChompOptions& operator=(const ChompOptions&) = default;
  ~ChompOptions() override =default;
  ChompOptions* CreateCopy() const override { return new ChompOptions(*this);}
  
  std::string GetDBName() const { return m_dbParam->GetStringValue(); }
  OEBioisostere::OEFragmentOptions* GetFragOpts() const { return m_FragmentOptions; }
  OEBioisostere::OEDBScreenOptions* GetScreenOpts() const { return m_ScreenOptions; }

private:
  OESystem::OEParameter* m_dbParam;
  OEBioisostere::OEFragmentOptions* m_FragmentOptions;
  OEBioisostere::OEDBScreenOptions* m_ScreenOptions;
};


int main(int argc, char* argv[])
{
  ChompOptions chompOpts;
  OEChem::OESimpleAppOptions opts(chompOpts, "BroodDataBase", OEChem::OEFileStringType::Mol);
  
  if (OESystem::OEConfigureOpts(opts, argc, argv, false) == OESystem::OEOptsConfigureStatus::Help)
    return 0;
  chompOpts.UpdateValues(opts);
  
  OEChem::oemolistream ifs;
  if (!ifs.open(opts.GetInFile()))
    OESystem::OEThrow.Fatal("Unable to open %s for reading", opts.GetInFile().c_str());
  
  printf("Generating fragments...\n");
  OEBioisostere::OEDBBuilder builder = OEBioisostere::OEDBBuilder( *chompOpts.GetFragOpts());
  OEChem::OEGraphMol mol;
  while(OEChem::OEReadMolecule(ifs, mol))
    builder.Generate(mol);
  
  printf("Building 2D fragments library...\n");
  builder.Filter(*OEBioisostere::OECreateFragFilter());
  builder.Expand(*OEBioisostere::OECreateFlipperOptions());
  builder.Screen(*chompOpts.GetScreenOpts());
  
  printf("Generating fragment conformers...\n");
  OEBioisostere::OEDBWriter writer;
  writer.Init(chompOpts.GetDBName());
  unsigned count=0;
  OESystem::OEIter<OEChem::OEMCMolBase> frag = builder.GetFrags();
  for(; frag; ++frag)
  {
    count++;
    OEBioisostere::OEGenerateConformers(*frag);
    writer.Write(*frag);
  }
  writer.Finish();
  printf("Generated fragments: %d \n",count);
}

Download code

brood_database.cpp

Creating Brood Query

The following code example shows how to create a Brood query from a molecule that could be used for bioisosteric fragment replacements using Brood.

See also

Listing 2: Creating Brood Query

/*
 (C) 2022 Cadence Design Systems, Inc. (Cadence)
 All rights reserved.
 TERMS FOR USE OF SAMPLE CODE The software below ("Sample Code") is
 provided to current licensees or subscribers of Cadence products or
 SaaS offerings (each a "Customer").
 Customer is hereby permitted to use, copy, and modify the Sample Code,
 subject to these terms. Cadence claims no rights to Customer's
 modifications. Modification of Sample Code is at Customer's sole and
 exclusive risk. Sample Code may require Customer to have a then
 current license or subscription to the applicable Cadence offering.
 THE SAMPLE CODE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
 EXPRESS OR IMPLIED.  OPENEYE DISCLAIMS ALL WARRANTIES, INCLUDING, BUT
 NOT LIMITED TO, WARRANTIES OF MERCHANTABILITY, FITNESS FOR A
 PARTICULAR PURPOSE AND NONINFRINGEMENT. In no event shall Cadence be
 liable for any damages or liability in connection with the Sample Code
 or its use.
 */
#include <cstdio>
#include <string>
#include <openeye.h>

#include <oeplatform.h>
#include <oesystem.h>
#include <oechem.h>
#include <oebioisostere.h>


class QueryOptions : public OESystem::OEOptions
{
public:
  QueryOptions(std::string name = "QueryOptions")
  :OESystem::OEOptions(name)
  {
    OESystem::OEUIntParameter idxParam("-atomIndices");
    idxParam.SetIsList(true);
    idxParam.SetRequired(true);
    idxParam.SetBrief("Index of atoms in fragment");
    m_idxParam = AddParameter(idxParam);
    
    OEChem::OEFileStringParameter duParam("-du",OEChem::OEFileStringType::DU);
    duParam.SetBrief("Design unit containg protein target for bump check");
    m_duParam = AddParameter(duParam);
    
    OESystem::OEUIntParameter maskParam("-proteinMask", OEBio::OEDesignUnitComponents::TargetComplexNoSolvent);
    maskParam.SetBrief("Design unit mask to identify protein target");
    m_maskParam = AddParameter(maskParam);
    
    OEChem::OEFileStringParameter selectDUParam("-selectDU",OEChem::OEFileStringType::DU);
    selectDUParam.SetBrief("Design unit containg select protein for bump check");
    m_selectDUParam = AddParameter(selectDUParam);
    
    OESystem::OEUIntParameter maskSelectParam("-proteinSelectMask", OEBio::OEDesignUnitComponents::TargetComplexNoSolvent);
    maskSelectParam.SetBrief("Design unit mask to identify select protein target");
    m_maskSelectParam = AddParameter(maskSelectParam);
  }
  
  QueryOptions(const QueryOptions&) = default;
  QueryOptions& operator=(const QueryOptions&) = default;
  ~QueryOptions() override =default;
  QueryOptions* CreateCopy() const override { return new QueryOptions(*this);}
  
  std::vector<int> GetIndices()
  {
    std::vector<int> indices;
    OESystem::OEIter<const std::string> Value = m_idxParam->GetStringValues();
    for (;Value;++Value)
    {
      indices.push_back(std::stoi(Value));
    }
    return indices;
  }
  
  bool GetDU(OEChem::OESimpleAppOptions opts, OEBio::OEDesignUnit& du)
  {
    OEPlatform::oeifstream ifs;
    if(!m_duParam->GetHasValue())
      return false;
    if(!ifs.open(m_duParam->GetStringValue()))
      OESystem::OEThrow.Fatal("Unable to open %s for reading", opts.GetInFile().c_str());

    if(!OEBio::OEReadDesignUnit(ifs, du))
      OESystem::OEThrow.Fatal("Unable to read design unit");
    return true;
  }
  
  int GetProteinMask()
  {
    std::string mask = m_maskParam->GetStringValue();
    if (mask == "")
      mask = m_maskParam->GetStringDefault();
    return std::stoi(mask);
  }

  bool GetSelectDU(OEChem::OESimpleAppOptions opts, OEBio::OEDesignUnit& selectDU)
  {
    OEPlatform::oeifstream ifs;
    if(!m_selectDUParam->GetHasValue())
      return false;
    if(!ifs.open(m_selectDUParam->GetStringValue()))
      OESystem::OEThrow.Fatal("Unable to open %s for reading", opts.GetInFile().c_str());
    
    if(!OEBio::OEReadDesignUnit(ifs, selectDU))
      OESystem::OEThrow.Fatal("Unable to read design unit");
    return true;
  }
  
  int GetSelectProteinMask()
  {
    std::string mask = m_maskSelectParam->GetStringValue();
    if (mask == "")
      mask = m_maskSelectParam->GetStringDefault();
    return std::stoi(mask);
  }
private:
  OESystem::OEParameter* m_idxParam;
  OESystem::OEParameter* m_duParam;
  OESystem::OEParameter* m_maskParam;
  OESystem::OEParameter* m_selectDUParam;
  OESystem::OEParameter* m_maskSelectParam;
};

int main(int argc, char* argv[])
{
  QueryOptions queryOpts;
  OEChem::OESimpleAppOptions opts(queryOpts, "BroodQuery", OEChem::OEFileStringType::Mol3D, "oeb");
  
  if (OESystem::OEConfigureOpts(opts, argc, argv, false) == OESystem::OEOptsConfigureStatus::Help)
    return 0;
  queryOpts.UpdateValues(opts);
  
  OEChem::oemolistream ifs;
  if (!ifs.open(opts.GetInFile()))
    OESystem::OEThrow.Fatal("Unable to open %s for reading", opts.GetInFile().c_str());
  
  OEChem::oemolostream ofs;
  if (!ofs.open(opts.GetOutFile()))
    OESystem::OEThrow.Fatal("Unable to open %s for writing", opts.GetOutFile().c_str());
  
  OEChem::OEMol queryMol;
  if (!OEChem::OEReadMolecule(ifs,queryMol))
    OESystem::OEThrow.Fatal("Unable to load molecule");
  
  std::vector<int> indices = queryOpts.GetIndices();
  std::vector<OEChem::OEAtomBase*> atoms;
  for (const auto& idx:indices)
  {
    OEChem::OEAtomBase* atom;
    atom = queryMol.GetAtom(OEChem::OEHasAtomIdx(idx));
    if (atom == nullptr)
      OESystem::OEThrow.Fatal("Invalid atom index %d", idx);
    atoms.push_back(std::move(atom));
  }

  OEChem::OEAtomBondSet selection;
  selection.AddAtoms(atoms);
  
  OEBioisostere::OEBroodQuery query;
  OEBio::OEDesignUnit du;
  OEBio::OEDesignUnit selectDU;
  unsigned retCode;
  
  if (!queryOpts.GetDU(opts, du) && !queryOpts.GetSelectDU(opts, selectDU))
    retCode = OEBioisostere::OECreateBroodQuery(query, queryMol, selection);
  else if (!queryOpts.GetDU(opts, du) && queryOpts.GetSelectDU(opts, selectDU))
  {
    const bool passingProteinSelect = true;
    retCode = OEBioisostere::OECreateBroodQuery(query, queryMol, selection, selectDU, queryOpts.GetSelectProteinMask(), passingProteinSelect);
  }
  else if (!queryOpts.GetSelectDU(opts, selectDU))
    retCode = OEBioisostere::OECreateBroodQuery(query, queryMol, selection, du, queryOpts.GetProteinMask());
  else
    retCode = OEBioisostere::OECreateBroodQuery(query, queryMol, selection, du, queryOpts.GetProteinMask(), selectDU, queryOpts.GetSelectProteinMask());
  
  if(retCode != OEBioisostere::OEBroodStatusCode::Success)
    OESystem::OEThrow.Fatal("%s", OEBioisostere::OEGetBroodStatus(retCode).c_str());
  
  OEBioisostere::OEWriteBroodQuery(ofs, query);
}

Download code

brood_query.cpp

Generating Brood Hits

The following code example shows how to perform bioisosteric fragment replacements on a Brood query and generate a hit list.

Listing 3: Generating Brood Hits

/*
 (C) 2022 Cadence Design Systems, Inc. (Cadence)
 All rights reserved.
 TERMS FOR USE OF SAMPLE CODE The software below ("Sample Code") is
 provided to current licensees or subscribers of Cadence products or
 SaaS offerings (each a "Customer").
 Customer is hereby permitted to use, copy, and modify the Sample Code,
 subject to these terms. Cadence claims no rights to Customer's
 modifications. Modification of Sample Code is at Customer's sole and
 exclusive risk. Sample Code may require Customer to have a then
 current license or subscription to the applicable Cadence offering.
 THE SAMPLE CODE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
 EXPRESS OR IMPLIED.  OPENEYE DISCLAIMS ALL WARRANTIES, INCLUDING, BUT
 NOT LIMITED TO, WARRANTIES OF MERCHANTABILITY, FITNESS FOR A
 PARTICULAR PURPOSE AND NONINFRINGEMENT. In no event shall Cadence be
 liable for any damages or liability in connection with the Sample Code
 or its use.
 */
#include <cstdio>
#include <openeye.h>

#include <oeplatform.h>
#include <oesystem.h>
#include <oechem.h>
#include <oebioisostere.h>


class BroodOptions : public OESystem::OEOptions
{
public:
  BroodOptions(std::string name = "BroodOptions")
  :OESystem::OEOptions(name)
  {
    OESystem::OEStringParameter pDBName("-db");
    pDBName.SetRequired(true);
    pDBName.SetVisibility(OESystem::OEParamVisibility::Simple);
    pDBName.SetBrief("Database Folder");
    m_dbParam = AddParameter(pDBName);
    
    m_GeneralOptions = static_cast<OEBioisostere::OEBroodGeneralOptions*>(AddOption(OEBioisostere::OEBroodGeneralOptions()));
    m_ScoreOptions = static_cast<OEBioisostere::OEBroodScoreOptions*>(AddOption(OEBioisostere::OEBroodScoreOptions()));
    m_HitlistOptions = static_cast<OEBioisostere::OEBroodHitlistOptions*>(AddOption(OEBioisostere::OEBroodHitlistOptions()));
  }
  
  BroodOptions(const BroodOptions&) = default;
  BroodOptions& operator=(const BroodOptions&) = default;
  ~BroodOptions() override =default;
  BroodOptions* CreateCopy() const override { return new BroodOptions(*this);}
  
  std::string GetDataBase() const { return m_dbParam->GetStringValue(); }
  OEBioisostere::OEBroodGeneralOptions* GetGenOpts() const { return m_GeneralOptions; }
  OEBioisostere::OEBroodScoreOptions* GetScoreOpts() const { return m_ScoreOptions; }
  OEBioisostere::OEBroodHitlistOptions* GetHitlistOpts() const { return m_HitlistOptions; }
  
private:
  OESystem::OEParameter* m_dbParam;
  OEBioisostere::OEBroodGeneralOptions* m_GeneralOptions;
  OEBioisostere::OEBroodScoreOptions* m_ScoreOptions;
  OEBioisostere::OEBroodHitlistOptions* m_HitlistOptions;
};

int main(int argc, char* argv[])
{
  BroodOptions broodOpts;
  OEChem::OESimpleAppOptions opts(broodOpts, "BroodHitlist", OEChem::OEFileStringType::Mol3D, OEChem::OEFileStringType::Mol3D);
  
  if (OESystem::OEConfigureOpts(opts, argc, argv, false) == OESystem::OEOptsConfigureStatus::Help)
    return 0;
  broodOpts.UpdateValues(opts);
  
  OEChem::oemolistream ifs;
  if (!ifs.open(opts.GetInFile()))
    OESystem::OEThrow.Fatal("Unable to open %s for reading", opts.GetInFile().c_str());
  
  OEChem::oemolostream ofs;
  if (!ofs.open(opts.GetOutFile()))
    OESystem::OEThrow.Fatal("Unable to open %s for writing", opts.GetOutFile().c_str());
  
  OEBioisostere::OEBroodQuery query;
  bool retCode = OEBioisostere::OEReadBroodQuery(ifs,query);
  if (retCode != OEBioisostere::OEBroodStatusCode::Success)
    OESystem::OEThrow.Fatal("Failed: %", OEBioisostere::OEGetBroodStatus(retCode).c_str());
  
  OEBioisostere::OEDBReader reader;
  unsigned retValue = reader.Init(broodOpts.GetDataBase(), query, *broodOpts.GetGenOpts());
  if (retValue != OEBioisostere::OEBroodStatusCode::Success)
    OESystem::OEThrow.Fatal("Unable to load Brood database");
  
  OEBioisostere::OEBroodOverlay overlay = OEBioisostere::OEBroodOverlay(*broodOpts.GetGenOpts(), *broodOpts.GetScoreOpts());
  overlay.SetupRef(query);
  
  OEBioisostere::OEHitlistBuilder hlist = OEBioisostere::OEHitlistBuilder(query, *broodOpts.GetGenOpts(), *broodOpts.GetHitlistOpts());
  
  unsigned packetCount = 0;
  OEBioisostere::OEBroodDBPacket packet;
  while(reader.GetNextPacket(packet))
  {
    packetCount++;
    printf("Processing packet %d with %d fragments \n", packetCount, packet.GetFragCount());
    std::vector<OEBioisostere::OEBroodScore> vecScores = overlay.Overlay(packet);
    hlist.AddScores(vecScores);
  }
  printf("Generating hitlist... \n");
  hlist.Build();
  printf("Total number of fragments overlayed: %d \n", hlist.GetAddCount());
  printf("Number of final hits: %d \n", hlist.GetHitCount());
  
  std::vector<OEBioisostere::OEBroodHit> vecHits = hlist.GetHits();
  unsigned index = 0;
  for(const auto& hit:vecHits)
  {
    OEChem::OEWriteConstMolecule(ofs, hit.GetMol());
    double comboScore = hit.GetComboScore();
    if ((*broodOpts.GetGenOpts()).GetScoreType() == OEBioisostere::OEBroodScoreType::ET)
      comboScore = hit.GetETComboScore();
    printf("Hit: %d %s Combo Score: %2f Belief Score: %.2f Complexity: %2f \n", index+1, hit.GetMol().GetTitle(), comboScore, hit.GetBeliefScore(), hit.GetComplexity());
    index++;
  }
}

Download code

brood_hitlist.cpp

Generating Brood Matches

The following code example shows how to perform bioisosteric fragment replacements on a Brood query and generate all possible matches. This could be a use case when generating all possible design ideas using BROOD and postprocessing them with other tools.

Listing 4: Generating Brood Matches

/*
 (C) 2022 Cadence Design Systems, Inc. (Cadence)
 All rights reserved.
 TERMS FOR USE OF SAMPLE CODE The software below ("Sample Code") is
 provided to current licensees or subscribers of Cadence products or
 SaaS offerings (each a "Customer").
 Customer is hereby permitted to use, copy, and modify the Sample Code,
 subject to these terms. Cadence claims no rights to Customer's
 modifications. Modification of Sample Code is at Customer's sole and
 exclusive risk. Sample Code may require Customer to have a then
 current license or subscription to the applicable Cadence offering.
 THE SAMPLE CODE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
 EXPRESS OR IMPLIED.  OPENEYE DISCLAIMS ALL WARRANTIES, INCLUDING, BUT
 NOT LIMITED TO, WARRANTIES OF MERCHANTABILITY, FITNESS FOR A
 PARTICULAR PURPOSE AND NONINFRINGEMENT. In no event shall Cadence be
 liable for any damages or liability in connection with the Sample Code
 or its use.
 */
#include <cstdio>
#include <string>
#include <openeye.h>

#include <oeplatform.h>
#include <oesystem.h>
#include <oechem.h>
#include <oebioisostere.h>


class BroodOptions : public OESystem::OEOptions
{
public:
  BroodOptions(std::string name = "BroodOptions")
  :OESystem::OEOptions(name)
  {
    OESystem::OEStringParameter pDBName("-db");
    pDBName.SetRequired(true);
    pDBName.SetVisibility(OESystem::OEParamVisibility::Simple);
    pDBName.SetBrief("Database Folder");
    m_dbParam = AddParameter(pDBName);
    
    m_GeneralOptions = static_cast<OEBioisostere::OEBroodGeneralOptions*>(AddOption(OEBioisostere::OEBroodGeneralOptions()));
    m_ScoreOptions = static_cast<OEBioisostere::OEBroodScoreOptions*>(AddOption(OEBioisostere::OEBroodScoreOptions()));
    m_BuildOptions = static_cast<OEBioisostere::OEBroodBuildOptions*>(AddOption(OEBioisostere::OEBroodBuildOptions()));
  }
  
  BroodOptions(const BroodOptions&) = default;
  BroodOptions& operator=(const BroodOptions&) = default;
  ~BroodOptions() override =default;
  BroodOptions* CreateCopy() const override { return new BroodOptions(*this);}
  
  std::string GetDataBase() const { return m_dbParam->GetStringValue(); }
  OEBioisostere::OEBroodGeneralOptions* GetGenOpts() const { return m_GeneralOptions; }
  OEBioisostere::OEBroodScoreOptions* GetScoreOpts() const { return m_ScoreOptions; }
  OEBioisostere::OEBroodBuildOptions* GetBuildOpts() const { return m_BuildOptions; }

private:
  OESystem::OEParameter* m_dbParam;
  OEBioisostere::OEBroodGeneralOptions* m_GeneralOptions;
  OEBioisostere::OEBroodScoreOptions* m_ScoreOptions;
  OEBioisostere::OEBroodBuildOptions* m_BuildOptions;
};


int main(int argc, char* argv[])
{
  BroodOptions broodOpts;
  OEChem::OESimpleAppOptions opts(broodOpts, "BroodMatching", OEChem::OEFileStringType::Mol3D, OEChem::OEFileStringType::Mol3D);
  
  if (OESystem::OEConfigureOpts(opts, argc, argv, false) == OESystem::OEOptsConfigureStatus::Help)
    return 0;
  broodOpts.UpdateValues(opts);
  
  OEChem::oemolistream ifs;
  if (!ifs.open(opts.GetInFile()))
    OESystem::OEThrow.Fatal("Unable to open %s for reading", opts.GetInFile().c_str());
  
  OEChem::oemolostream ofs;
  if (!ofs.open(opts.GetOutFile()))
    OESystem::OEThrow.Fatal("Unable to open %s for writing", opts.GetOutFile().c_str());
  
  OEBioisostere::OEBroodQuery query;
  bool retCode = OEBioisostere::OEReadBroodQuery(ifs,query);
  if (retCode != OEBioisostere::OEBroodStatusCode::Success)
    OESystem::OEThrow.Fatal("Failed: %s", OEBioisostere::OEGetBroodStatus(retCode).c_str());
  
  OEBioisostere::OEDBReader reader;
  unsigned retValue = reader.Init(broodOpts.GetDataBase(), query, *broodOpts.GetGenOpts());
  if (retValue != OEBioisostere::OEBroodStatusCode::Success)
    OESystem::OEThrow.Fatal("Unable to load Brood database");
  
  OEBioisostere::OEBroodOverlay overlay = OEBioisostere::OEBroodOverlay(*broodOpts.GetGenOpts(), *broodOpts.GetScoreOpts());
  overlay.SetupRef(query);
  
  OEBioisostere::OEBroodMolBuilder builder = OEBioisostere::OEBroodMolBuilder(query, *broodOpts.GetGenOpts(), *broodOpts.GetBuildOpts());
  
  unsigned packetCount = 0;
  unsigned totalCount = 0;
  unsigned successCount = 0;
  OEBioisostere::OEBroodDBPacket packet;
  while(reader.GetNextPacket(packet))
  {
    packetCount++;
    printf("Processing packet %d with %d fragments \n", packetCount, packet.GetFragCount());
    std::vector<OEBioisostere::OEBroodScore> vecScores = overlay.Overlay(packet);
    for (const auto& score:vecScores)
    {
      totalCount++;
      if (score.GetStatus() == OEBioisostere::OEBroodStatusCode::Success)
      {
        OEBioisostere::OEBroodHit hit;
        if(builder.Build(score,hit)==OEBioisostere::OEBroodStatusCode::Success)
        {
          OEChem::OEWriteConstMolecule(ofs, hit.GetMol());
          successCount++;
        }
      }
    }
  }
  ofs.close();
  printf("Total number of fragments overlayed: %d \n",totalCount);
  printf("Number of successful matches: %d \n",successCount);
}

Download code

brood_matching.cpp

Overlay between Fragments

The following code example shows how to overlay a fragment against a query fragment. This could be a use case when working with synthons and trying to find similar synthons based on 3D similarity.

Listing 5: Overlay between Fragments

/*
 (C) 2022 Cadence Design Systems, Inc. (Cadence)
 All rights reserved.
 TERMS FOR USE OF SAMPLE CODE The software below ("Sample Code") is
 provided to current licensees or subscribers of Cadence products or
 SaaS offerings (each a "Customer").
 Customer is hereby permitted to use, copy, and modify the Sample Code,
 subject to these terms. Cadence claims no rights to Customer's
 modifications. Modification of Sample Code is at Customer's sole and
 exclusive risk. Sample Code may require Customer to have a then
 current license or subscription to the applicable Cadence offering.
 THE SAMPLE CODE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
 EXPRESS OR IMPLIED.  OPENEYE DISCLAIMS ALL WARRANTIES, INCLUDING, BUT
 NOT LIMITED TO, WARRANTIES OF MERCHANTABILITY, FITNESS FOR A
 PARTICULAR PURPOSE AND NONINFRINGEMENT. In no event shall Cadence be
 liable for any damages or liability in connection with the Sample Code
 or its use.
 */
#include <cstdio>
#include <string>
#include <openeye.h>

#include <oeplatform.h>
#include <oesystem.h>
#include <oechem.h>
#include <oebioisostere.h>


int main(int argc, char* argv[])
{
  OEBioisostere::OEBroodGeneralOptions genOpts;
  OEChem::OERefInputAppOptions opts(genOpts, "FragmentOverlay", OEChem::OEFileStringType::Mol3D, OEChem::OEFileStringType::Mol3D, OEChem::OEFileStringType::Mol3D, "-queryFrag");

  if (OESystem::OEConfigureOpts(opts, argc, argv, false) == OESystem::OEOptsConfigureStatus::Help)
    return 0;
  genOpts.UpdateValues(opts);

  OEChem::oemolistream ifs;
  if (!ifs.open(opts.GetInFile()))
    OESystem::OEThrow.Fatal("Unable to open %s for reading", opts.GetInFile().c_str());

  OEChem::oemolistream rfs;
  if (!rfs.open(opts.GetRefFile()))
    OESystem::OEThrow.Fatal("Unable to open %s for reading", opts.GetInFile().c_str());

  OEChem::oemolostream ofs;
  if (!ofs.open(opts.GetOutFile()))
    OESystem::OEThrow.Fatal("Unable to open %s for writing", opts.GetOutFile().c_str());

  OEChem::OEMol query;
  if (!OEChem::OEReadMolecule(rfs, query))
    OESystem::OEThrow.Fatal("Failed to read Query fragment");

  OEBioisostere::OEBroodFragPrep prep;
  OEBioisostere::OEFragOverlay overlay = OEBioisostere::OEFragOverlay(genOpts,OEBioisostere::OEBroodScoreOptions());
  overlay.SetupRef(query);
  
  OEChem::OEMol mol;
  while (OEChem::OEReadMolecule(ifs, mol))
  {
    printf("Overlaying %s \n",mol.GetTitle());
    prep.Prep(mol);
    OEBioisostere::OEBroodScore score;
    unsigned retCode = overlay.Overlay(mol, score);
    if (retCode == OEBioisostere::OEBroodStatusCode::Success)
    {
      OEChem::OEGraphMol outmol;
      score.GetFragMol(outmol);
      overlay.Transform(outmol);
      double comboScore = score.GetComboScore();
      if(genOpts.GetScoreType() == OEBioisostere::OEBroodScoreType::ET)
        comboScore = score.GetETComboScore();
      printf("Fragment: %s Combo Score: %.5f \n",mol.GetTitle(),comboScore);
      OEChem::OEWriteConstMolecule(ofs,outmol);
    }
    else
    {
      OESystem::OEThrow.Warning("%s: %s",mol.GetTitle(), OEBioisostere::OEGetBroodStatus(retCode).c_str());
    }
  }
  return 0;
}

Download code

fragment_overlay.cpp

Replacing a Fragment in a Molecule

The following code example shows how to replace a fragment in a molecule defined in the form of a BROOD query.

Listing 6: Replacing a Fragment in a Molecule

/*
 (C) 2022 Cadence Design Systems, Inc. (Cadence)
 All rights reserved.
 TERMS FOR USE OF SAMPLE CODE The software below ("Sample Code") is
 provided to current licensees or subscribers of Cadence products or
 SaaS offerings (each a "Customer").
 Customer is hereby permitted to use, copy, and modify the Sample Code,
 subject to these terms. Cadence claims no rights to Customer's
 modifications. Modification of Sample Code is at Customer's sole and
 exclusive risk. Sample Code may require Customer to have a then
 current license or subscription to the applicable Cadence offering.
 THE SAMPLE CODE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
 EXPRESS OR IMPLIED.  OPENEYE DISCLAIMS ALL WARRANTIES, INCLUDING, BUT
 NOT LIMITED TO, WARRANTIES OF MERCHANTABILITY, FITNESS FOR A
 PARTICULAR PURPOSE AND NONINFRINGEMENT. In no event shall Cadence be
 liable for any damages or liability in connection with the Sample Code
 or its use.
 */
#include <cstdio>
#include <string>
#include <openeye.h>

#include <oeplatform.h>
#include <oesystem.h>
#include <oechem.h>
#include <oebioisostere.h>

class BroodOptions : public OESystem::OEOptions
{
public:
  BroodOptions(std::string name = "BroodOptions")
  :OESystem::OEOptions(name)
  {
    m_GeneralOptions = static_cast<OEBioisostere::OEBroodGeneralOptions*>(AddOption(OEBioisostere::OEBroodGeneralOptions()));
    m_ScoreOptions = static_cast<OEBioisostere::OEBroodScoreOptions*>(AddOption(OEBioisostere::OEBroodScoreOptions()));
    m_BuildOptions = static_cast<OEBioisostere::OEBroodBuildOptions*>(AddOption(OEBioisostere::OEBroodBuildOptions()));
  }
  
  BroodOptions(const BroodOptions&) = default;
  BroodOptions& operator=(const BroodOptions&) = default;
  ~BroodOptions() override =default;
  BroodOptions* CreateCopy() const override { return new BroodOptions(*this);}
  OEBioisostere::OEBroodGeneralOptions* GetGenOpts() const { return m_GeneralOptions; }
  OEBioisostere::OEBroodScoreOptions* GetScoreOpts() const { return m_ScoreOptions; }
  OEBioisostere::OEBroodBuildOptions* GetBuildOpts() const { return m_BuildOptions; }
  
private:
  OEBioisostere::OEBroodGeneralOptions* m_GeneralOptions;
  OEBioisostere::OEBroodScoreOptions* m_ScoreOptions;
  OEBioisostere::OEBroodBuildOptions* m_BuildOptions;
};


int main(int argc, char* argv[])
{
  BroodOptions broodOpts;
  OEChem::OERefInputAppOptions opts(broodOpts, "ReplaceFragment", OEChem::OEFileStringType::Mol3D,OEChem::OEFileStringType::Mol3D, OEChem::OEFileStringType::Mol3D, "-query");
  
  if (OESystem::OEConfigureOpts(opts, argc, argv, false) == OESystem::OEOptsConfigureStatus::Help)
    return 0;
  broodOpts.UpdateValues(opts);
  
  OEChem::oemolistream ifs;
  if (!ifs.open(opts.GetInFile()))
    OESystem::OEThrow.Fatal("Unable to open %s for reading", opts.GetInFile().c_str());
  
  OEChem::oemolistream rfs;
  if (!rfs.open(opts.GetRefFile()))
    OESystem::OEThrow.Fatal("Unable to open %s for reading", opts.GetRefFile().c_str());
  
  OEChem::oemolostream ofs;
  if (!ofs.open(opts.GetOutFile()))
    OESystem::OEThrow.Fatal("Unable to open %s for writing", opts.GetOutFile().c_str());
  
  OEBioisostere::OEBroodQuery query;
  bool retCode = OEBioisostere::OEReadBroodQuery(rfs,query);
  if (retCode != OEBioisostere::OEBroodStatusCode::Success)
    OESystem::OEThrow.Fatal("Failed: %s", OEBioisostere::OEGetBroodStatus(retCode).c_str());
  
  OEBioisostere::OEBroodFragPrep prep;
  OEBioisostere::OEBroodOverlay overlay = OEBioisostere::OEBroodOverlay(*broodOpts.GetGenOpts(), *broodOpts.GetScoreOpts());
  overlay.SetupRef(query);
  
  OEBioisostere::OEBroodMolBuilder builder = OEBioisostere::OEBroodMolBuilder(query, *broodOpts.GetGenOpts(), *broodOpts.GetBuildOpts());
  
  OEChem::OEMol mol;
  while (OEChem::OEReadMolecule(ifs, mol))
  {
    printf("Replacement fragment %s \n",mol.GetTitle());
    prep.Prep(mol);
    OEBioisostere::OEBroodScore score;
    unsigned retCode = overlay.Overlay(mol, score);
    if (retCode == OEBioisostere::OEBroodStatusCode::Success)
    {
      double comboScore = score.GetComboScore();
      if(broodOpts.GetGenOpts()->GetScoreType() == OEBioisostere::OEBroodScoreType::ET)
        comboScore = score.GetETComboScore();
      printf("Fragment: %s Combo Score: %.5f \n",mol.GetTitle(),comboScore);
      OEBioisostere::OEBroodHit hit;
      if(builder.Build(score, hit) == OEBioisostere::OEBroodStatusCode::Success)
        OEChem::OEWriteConstMolecule(ofs, hit.GetMol());
    }
    else
    {
      OESystem::OEThrow.Warning("%s: %s",mol.GetTitle(), OEBioisostere::OEGetBroodStatus(retCode).c_str());
    }
  }
  
  return 0;
}

Download code

replace_fragment.cpp

Building and Using an External Compound Database

The following code example shows how to build and use an external (in-house) compound database to find similar 2D compounds to a generated hit list.

Listing 7: Building and using an external Compound Database

/*
 (C) 2022 Cadence Design Systems, Inc. (Cadence)
 All rights reserved.
 TERMS FOR USE OF SAMPLE CODE The software below ("Sample Code") is
 provided to current licensees or subscribers of Cadence products or
 SaaS offerings (each a "Customer").
 Customer is hereby permitted to use, copy, and modify the Sample Code,
 subject to these terms. Cadence claims no rights to Customer's
 modifications. Modification of Sample Code is at Customer's sole and
 exclusive risk. Sample Code may require Customer to have a then
 current license or subscription to the applicable Cadence offering.
 THE SAMPLE CODE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
 EXPRESS OR IMPLIED.  OPENEYE DISCLAIMS ALL WARRANTIES, INCLUDING, BUT
 NOT LIMITED TO, WARRANTIES OF MERCHANTABILITY, FITNESS FOR A
 PARTICULAR PURPOSE AND NONINFRINGEMENT. In no event shall Cadence be
 liable for any damages or liability in connection with the Sample Code
 or its use.
 */
#include <cstdio>
#include <openeye.h>

#include <oeplatform.h>
#include <oesystem.h>
#include <oechem.h>
#include <oebioisostere.h>


class BroodOptions : public OESystem::OEOptions
{
public:
  BroodOptions(std::string name = "BroodOptions")
  :OESystem::OEOptions(name)
  {
    OESystem::OEStringParameter pDBName("-db");
    pDBName.SetRequired(true);
    pDBName.SetVisibility(OESystem::OEParamVisibility::Simple);
    pDBName.SetBrief("Brood Database Directory");
    m_dbParam = AddParameter(pDBName);
    
    OESystem::OEStringParameter pCPDDBName("-cpddb");
    pCPDDBName.SetRequired(true);
    pCPDDBName.SetVisibility(OESystem::OEParamVisibility::Simple);
    pCPDDBName.SetBrief("CPDDatabase File (database of known compounds to identify available compounds similar to hits)");
    m_cpddbParam = AddParameter(pCPDDBName);
    
    m_GeneralOptions = static_cast<OEBioisostere::OEBroodGeneralOptions*>(AddOption(OEBioisostere::OEBroodGeneralOptions()));
    m_ScoreOptions = static_cast<OEBioisostere::OEBroodScoreOptions*>(AddOption(OEBioisostere::OEBroodScoreOptions()));
    m_HitlistOptions = static_cast<OEBioisostere::OEBroodHitlistOptions*>(AddOption(OEBioisostere::OEBroodHitlistOptions()));
  }
  
  BroodOptions(const BroodOptions&) = default;
  BroodOptions& operator=(const BroodOptions&) = default;
  ~BroodOptions() override =default;
  BroodOptions* CreateCopy() const override { return new BroodOptions(*this);}
  
  std::string GetDataBase() const { return m_dbParam->GetStringValue(); }
  std::string GetCPDDataBase() const { return m_cpddbParam->GetStringValue(); }
  OEBioisostere::OEBroodGeneralOptions* GetGenOpts() const { return m_GeneralOptions; }
  OEBioisostere::OEBroodScoreOptions* GetScoreOpts() const { return m_ScoreOptions; }
  OEBioisostere::OEBroodHitlistOptions* GetHitlistOpts() const { return m_HitlistOptions; }
  
private:
  OESystem::OEParameter* m_dbParam;
  OESystem::OEParameter* m_cpddbParam;
  OEBioisostere::OEBroodGeneralOptions* m_GeneralOptions;
  OEBioisostere::OEBroodScoreOptions* m_ScoreOptions;
  OEBioisostere::OEBroodHitlistOptions* m_HitlistOptions;
};

int main(int argc, char* argv[])
{
  BroodOptions broodOpts;
  OEChem::OESimpleAppOptions opts(broodOpts, "BroodCPDDB", OEChem::OEFileStringType::Mol3D, "oeb");
  
  if (OESystem::OEConfigureOpts(opts, argc, argv, false) == OESystem::OEOptsConfigureStatus::Help)
    return 0;
  broodOpts.UpdateValues(opts);
  
  OEChem::oemolistream ifs;
  if (!ifs.open(opts.GetInFile()))
    OESystem::OEThrow.Fatal("Unable to open %s for reading", opts.GetInFile().c_str());
  
  OEBioisostere::OEBroodQuery query;
  bool retCode = OEBioisostere::OEReadBroodQuery(ifs,query);
  if (retCode != OEBioisostere::OEBroodStatusCode::Success)
    OESystem::OEThrow.Fatal("Failed: %", OEBioisostere::OEGetBroodStatus(retCode).c_str());
  
  OEBioisostere::OEDBReader reader;
  unsigned retValue = reader.Init(broodOpts.GetDataBase(), query, *broodOpts.GetGenOpts());
  if (retValue != OEBioisostere::OEBroodStatusCode::Success)
    OESystem::OEThrow.Fatal("Unable to load Brood database");
  
  OEChem::oemolistream ifsCPDDB;
  if (!ifsCPDDB.open(broodOpts.GetCPDDataBase()))
    OESystem::OEThrow.Fatal("Unable to load Brood CPDDatabase file");
  
  OEBioisostere::OECPDDatabase cpddb;
  retValue = cpddb.Prep(ifsCPDDB);
  OESystem::OEThrow.Info("Prepration Mode: %s", OEBioisostere::OEGetBroodStatus(retValue).c_str());
  
  int numOfCPDDBMolecule = cpddb.GetNumMolecules();
  OESystem::OEThrow.Info("Number of molecules in CPDDatabse: %d", numOfCPDDBMolecule);
  
  if(retValue == OEBioisostere::OEBroodStatusCode::NewFPGenerated){
    OEChem::oemolostream ofsCPDDB;
    if (!ofsCPDDB.open(opts.GetOutFile()))
      OESystem::OEThrow.Fatal("Unable to open %s for writing the newly generated finegrprint" , opts.GetOutFile().c_str());
    
    if(cpddb.Write(ifsCPDDB, ofsCPDDB))
      OESystem::OEThrow.Info("Fingerprints are written in the %s. You can use this file in the future with -cpddb for improved efficiency over current useage.",ofsCPDDB.GetFileName().c_str());
    else
      OESystem::OEThrow.Warning("Fingerprint mismatch in writing updated -cpddb file.",ofsCPDDB.GetFileName().c_str());
  }
  
  OEBioisostere::OEBroodOverlay overlay = OEBioisostere::OEBroodOverlay(*broodOpts.GetGenOpts(), *broodOpts.GetScoreOpts());
  overlay.SetupRef(query);
  
  OEBioisostere::OEHitlistBuilder hlist = OEBioisostere::OEHitlistBuilder(query, *broodOpts.GetGenOpts(), *broodOpts.GetHitlistOpts());
  
  OEBioisostere::OEBroodDBPacket packet;
  while(reader.GetNextPacket(packet))
  {
    std::vector<OEBioisostere::OEBroodScore> vecScores = overlay.Overlay(packet);
    hlist.AddScores(vecScores);
  }
  printf("Generating hitlist... \n");
  hlist.Build();
  printf("Number of final hits: %d \n \n", hlist.GetHitCount());
  
  std::vector<OEBioisostere::OEBroodHit> vecHits = hlist.GetHits();
  
  for(const auto& hit:vecHits)
  {
    if (cpddb.IsPrepared()){
      std::vector<std::string> vecCpddbValues;
      cpddb.GetSimilarMolecules(hit.GetMol(), vecCpddbValues);
      if (!vecCpddbValues[0].empty())
        printf("Hit Mol (SMILES): %s \nAnalog Mol (SMILES): %s \nAnalog Mol Label: %s \n \n", OEChem::OEMolToSmiles(hit.GetMol()).c_str(),vecCpddbValues[0].c_str(),vecCpddbValues[1].c_str());
    }
  }
}

Download code

brood_cpddb.cpp

Creating Brood Query by Linking Two Molecules (Bridging)

The following code example shows how to create a Brood query from two molecules to find a suitable linker between the two using bioisosteric fragment replacements using Brood.

See also

Listing 8: Creating Brood Query by Linking Two Molecules (Bridging)

/*
 (C) 2022 Cadence Design Systems, Inc. (Cadence)
 All rights reserved.
 TERMS FOR USE OF SAMPLE CODE The software below ("Sample Code") is
 provided to current licensees or subscribers of Cadence products or
 SaaS offerings (each a "Customer").
 Customer is hereby permitted to use, copy, and modify the Sample Code,
 subject to these terms. Cadence claims no rights to Customer's
 modifications. Modification of Sample Code is at Customer's sole and
 exclusive risk. Sample Code may require Customer to have a then
 current license or subscription to the applicable Cadence offering.
 THE SAMPLE CODE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
 EXPRESS OR IMPLIED.  OPENEYE DISCLAIMS ALL WARRANTIES, INCLUDING, BUT
 NOT LIMITED TO, WARRANTIES OF MERCHANTABILITY, FITNESS FOR A
 PARTICULAR PURPOSE AND NONINFRINGEMENT. In no event shall Cadence be
 liable for any damages or liability in connection with the Sample Code
 or its use.
 */
#include <cstdio>
#include <string>
#include <openeye.h>

#include <oeplatform.h>
#include <oesystem.h>
#include <oechem.h>
#include <oebioisostere.h>


class QueryOptions : public OESystem::OEOptions
{
public:
  QueryOptions(std::string name = "2MolLinkerQueryOptions")
  :OESystem::OEOptions(name)
  {
    OESystem::OEUIntParameter idxNo1Param("-atomIndicesNo1");
    idxNo1Param.SetIsList(true);
    idxNo1Param.SetRequired(true);
    idxNo1Param.SetBrief("Index of atoms in fragment No.1");
    m_idxNo1Param = AddParameter(idxNo1Param);
    
    OESystem::OEUIntParameter idxNo2Param("-atomIndicesNo2");
    idxNo2Param.SetIsList(true);
    idxNo2Param.SetRequired(true);
    idxNo2Param.SetBrief("Index of atoms in fragment No.2");
    m_idxNo2Param = AddParameter(idxNo2Param);
    
    OEChem::OEFileStringParameter duParam("-du",OEChem::OEFileStringType::DU);
    duParam.SetBrief("Design unit containg protein target for bump check");
    m_duParam = AddParameter(duParam);
    
    OESystem::OEUIntParameter maskParam("-proteinMask", OEBio::OEDesignUnitComponents::TargetComplexNoSolvent);
    maskParam.SetBrief("Design unit mask to identify protein target");
    m_maskParam = AddParameter(maskParam);
    
    OEChem::OEFileStringParameter selectDUParam("-selectDU",OEChem::OEFileStringType::DU);
    selectDUParam.SetBrief("Design unit containg select protein for bump check");
    m_selectDUParam = AddParameter(selectDUParam);
    
    OESystem::OEUIntParameter maskSelectParam("-proteinSelectMask", OEBio::OEDesignUnitComponents::TargetComplexNoSolvent);
    maskSelectParam.SetBrief("Design unit mask to identify select protein target");
    m_maskSelectParam = AddParameter(maskSelectParam);
  }
  
  QueryOptions(const QueryOptions&) = default;
  QueryOptions& operator=(const QueryOptions&) = default;
  ~QueryOptions() override =default;
  QueryOptions* CreateCopy() const override { return new QueryOptions(*this);}
  
  std::vector<int> GetFirstMolIndices()
  {
    std::vector<int> indices;
    OESystem::OEIter<const std::string> Value = m_idxNo1Param->GetStringValues();
    for (;Value;++Value)
    {
      indices.push_back(std::stoi(Value));
    }
    return indices;
  }
  
  std::vector<int> GetSecondMolIndices()
  {
    std::vector<int> indices;
    OESystem::OEIter<const std::string> Value = m_idxNo2Param->GetStringValues();
    for (;Value;++Value)
    {
      indices.push_back(std::stoi(Value));
    }
    return indices;
  }
  
  bool GetDU(OEChem::OERefInputAppOptions opts, OEBio::OEDesignUnit& du)
  {
    OEPlatform::oeifstream ifs;
    if(!m_duParam->GetHasValue())
      return false;
    if(!ifs.open(m_duParam->GetStringValue()))
      OESystem::OEThrow.Fatal("Unable to open %s for reading", opts.GetInFile().c_str());

    if(!OEBio::OEReadDesignUnit(ifs, du))
      OESystem::OEThrow.Fatal("Unable to read design unit");
    return true;
  }
  
  int GetProteinMask()
  {
    std::string mask = m_maskParam->GetStringValue();
    if (mask == "")
      mask = m_maskParam->GetStringDefault();
    return std::stoi(mask);
  }

  bool GetSelectDU(OEChem::OERefInputAppOptions opts, OEBio::OEDesignUnit& selectDU)
  {
    OEPlatform::oeifstream ifs;
    if(!m_selectDUParam->GetHasValue())
      return false;
    if(!ifs.open(m_selectDUParam->GetStringValue()))
      OESystem::OEThrow.Fatal("Unable to open %s for reading", opts.GetInFile().c_str());
    
    if(!OEBio::OEReadDesignUnit(ifs, selectDU))
      OESystem::OEThrow.Fatal("Unable to read design unit");
    return true;
  }
  
  int GetSelectProteinMask()
  {
    std::string mask = m_maskSelectParam->GetStringValue();
    if (mask == "")
      mask = m_maskSelectParam->GetStringDefault();
    return std::stoi(mask);
  }
private:
  OESystem::OEParameter* m_idxNo1Param;
  OESystem::OEParameter* m_idxNo2Param;
  OESystem::OEParameter* m_duParam;
  OESystem::OEParameter* m_maskParam;
  OESystem::OEParameter* m_selectDUParam;
  OESystem::OEParameter* m_maskSelectParam;
};

int main(int argc, char* argv[])
{
  QueryOptions queryOpts;
  OEChem::OERefInputAppOptions opts(queryOpts, "2MolLinkerBroodQuery", OEChem::OEFileStringType::Mol3D, OEChem::OEFileStringType::Mol3D,
                                    OEChem::OEFileStringType::Mol3D, "-in2");
  
  if (OESystem::OEConfigureOpts(opts, argc, argv, false) == OESystem::OEOptsConfigureStatus::Help)
    return 0;
  queryOpts.UpdateValues(opts);
  
  OEChem::oemolistream ifs;
  if (!ifs.open(opts.GetInFile()))
    OESystem::OEThrow.Fatal("Unable to open %s for reading", opts.GetInFile().c_str());
  
  OEChem::oemolistream rfs;
  if (!rfs.open(opts.GetRefFile()))
    OESystem::OEThrow.Fatal("Unable to open %s for reading", opts.GetRefFile().c_str());
  
  OEChem::oemolostream ofs;
  if (!ofs.open(opts.GetOutFile()))
    OESystem::OEThrow.Fatal("Unable to open %s for writing", opts.GetOutFile().c_str());
  
  OEChem::OEMol firstMol;
  if (!OEChem::OEReadMolecule(ifs,firstMol))
    OESystem::OEThrow.Fatal("Unable to load molecule 1");
  
  OEChem::OEMol secondMol;
  if (!OEChem::OEReadMolecule(rfs, secondMol))
    OESystem::OEThrow.Fatal("Unable to load molecule 2");
  
  std::vector<int> firstMolIndices = queryOpts.GetFirstMolIndices();
  std::vector<OEChem::OEAtomBase*> firstMolAtoms;
  for (const auto& idx:firstMolIndices)
  {
    OEChem::OEAtomBase* firstMolAtom;
    firstMolAtom = firstMol.GetAtom(OEChem::OEHasAtomIdx(idx));
    if (firstMolAtom == nullptr)
      OESystem::OEThrow.Fatal("Invalid atom index %d", idx);
    firstMolAtoms.push_back(std::move(firstMolAtom));
  }

  std::vector<int> secondMolIndices = queryOpts.GetSecondMolIndices();
  std::vector<OEChem::OEAtomBase*> secondMolAtoms;
  for (const auto& idx:secondMolIndices)
  {
    OEChem::OEAtomBase* secondMolAtom;
    secondMolAtom = secondMol.GetAtom(OEChem::OEHasAtomIdx(idx));
    if (secondMolAtom == nullptr)
      OESystem::OEThrow.Fatal("Invalid atom index %d", idx);
    secondMolAtoms.push_back(std::move(secondMolAtom));
  }
  
  OEChem::OEAtomBondSet firstMolSelection;
  firstMolSelection.AddAtoms(firstMolAtoms);
  
  OEChem::OEAtomBondSet secondMolSelection;
  secondMolSelection.AddAtoms(secondMolAtoms);
  
  
  OEBioisostere::OEBroodQuery query;
  OEBio::OEDesignUnit du;
  OEBio::OEDesignUnit selectDU;
  
  unsigned retCode;
  
  if (!queryOpts.GetDU(opts, du) && !queryOpts.GetSelectDU(opts, selectDU))
    retCode = OEBioisostere::OECreateBroodQuery(query, firstMol, secondMol, firstMolSelection, secondMolSelection);
  else if (!queryOpts.GetDU(opts, du) && queryOpts.GetSelectDU(opts, selectDU))
  {
    const bool passingProteinSelect = true;
    retCode = OEBioisostere::OECreateBroodQuery(query, firstMol, secondMol, firstMolSelection, secondMolSelection, selectDU, queryOpts.GetSelectProteinMask(), passingProteinSelect);
  }
  else if (!queryOpts.GetSelectDU(opts, selectDU))
    retCode = OEBioisostere::OECreateBroodQuery(query, firstMol, secondMol, firstMolSelection, secondMolSelection, du, queryOpts.GetProteinMask());
  else
    retCode = OEBioisostere::OECreateBroodQuery(query, firstMol, secondMol, firstMolSelection, secondMolSelection, du, queryOpts.GetProteinMask(), selectDU, queryOpts.GetSelectProteinMask());

  if(retCode != OEBioisostere::OEBroodStatusCode::Success)
    OESystem::OEThrow.Fatal("%s", OEBioisostere::OEGetBroodStatus(retCode).c_str());
  
  OEBioisostere::OEWriteBroodQuery(ofs, query);
}