开通博客赚积分
发布资源赚积分
分类

源码开发语言/平台
当前位置：首页 > 源码/资料 > 数值算法/人工智能 > 生物技术 > 查看源码
chemical_graph.cpp
资源名称：ncbi_cxx_win--Jul_08_2004.zip　[点击查看]
上传用户：yhdzpy8989
上传日期：2007-06-13
资源大小：13604k
文件大小：27k
源码类别：

生物技术

开发平台：

C/C++

chemical_graph.cpp：源码内容
  1. /*
  2.  * ===========================================================================
  3.  * PRODUCTION $Log: chemical_graph.cpp,v $
  4.  * PRODUCTION Revision 1000.2  2004/06/01 18:28:03  gouriano
  5.  * PRODUCTION PRODUCTION: UPGRADED [GCC34_MSVC7] Dev-tree R1.40
  6.  * PRODUCTION
  7.  * ===========================================================================
  8.  */
  10. /*  $Id: chemical_graph.cpp,v 1000.2 2004/06/01 18:28:03 gouriano Exp $
  11. * ===========================================================================
  12. *
  13. *                            PUBLIC DOMAIN NOTICE
  14. *               National Center for Biotechnology Information
  15. *
  16. *  This software/database is a "United States Government Work" under the
  17. *  terms of the United States Copyright Act.  It was written as part of
  18. *  the author's official duties as a United States Government employee and
  19. *  thus cannot be copyrighted.  This software/database is freely available
  20. *  to the public for use. The National Library of Medicine and the U.S.
  21. *  Government have not placed any restriction on its use or reproduction.
  22. *
  23. *  Although all reasonable efforts have been taken to ensure the accuracy
  24. *  and reliability of the software and data, the NLM and the U.S.
  25. *  Government do not and cannot warrant the performance or results that
  26. *  may be obtained by using this software or data. The NLM and the U.S.
  27. *  Government disclaim all warranties, express or implied, including
  28. *  warranties of performance, merchantability or fitness for any particular
  29. *  purpose.
  30. *
  31. *  Please cite the author in any work or product based on this material.
  32. *
  33. * ===========================================================================
  34. *
  35. * Authors:  Paul Thiessen
  36. *
  37. * File Description:
  38. *      Classes to hold the graph of chemical bonds
  39. *
  40. * ===========================================================================
  41. */
  43. #ifdef _MSC_VER
  44. #pragma warning(disable:4018)   // disable signed/unsigned mismatch warning in MSVC
  45. #endif
  47. #include <ncbi_pch.hpp>
  48. #include <corelib/ncbistd.hpp>
  50. #include <objects/mmdb1/Biostruc_residue_graph_set.hpp>
  51. #include <objects/mmdb1/Biostruc_id.hpp>
  52. #include <objects/general/Dbtag.hpp>
  53. #include <objects/mmdb3/Biostruc_feature_set_descr.hpp>
  54. #include <objects/mmdb3/Biostruc_feature.hpp>
  55. #include <objects/mmdb3/Chem_graph_pntrs.hpp>
  56. #include <objects/mmdb3/Residue_pntrs.hpp>
  57. #include <objects/mmdb3/Residue_interval_pntr.hpp>
  58. #include <objects/mmdb3/Biostruc_feature_id.hpp>
  59. #include <objects/mmdb1/Molecule_id.hpp>
  60. #include <objects/mmdb1/Residue_id.hpp>
  61. #include <objects/mmdb1/Atom_id.hpp>
  62. #include <objects/general/Object_id.hpp>
  64. #include "chemical_graph.hpp"
  65. #include "asn_reader.hpp"
  66. #include "molecule.hpp"
  67. #include "bond.hpp"
  68. #include "structure_set.hpp"
  69. #include "opengl_renderer.hpp"
  70. #include "coord_set.hpp"
  71. #include "atom_set.hpp"
  72. #include "object_3d.hpp"
  73. #include "cn3d_tools.hpp"
  74. #include "molecule_identifier.hpp"
  76. USING_NCBI_SCOPE;
  77. USING_SCOPE(objects);
  80. BEGIN_SCOPE(Cn3D)
  82. static const CBiostruc_residue_graph_set* standardDictionary = NULL;
  84. void LoadStandardDictionary(const char *filename)
  85. {
  86.     standardDictionary = new CBiostruc_residue_graph_set;
  88.     // read dictionary file
  89.     string err;
  90.     if (!ReadASNFromFile(filename,
  91.             const_cast<CBiostruc_residue_graph_set*>(standardDictionary), true, &err))
  92.         FATALMSG("Error reading standard residue dictionary: " << err);
  94.     // make sure it's the right thing
  95.     if (!standardDictionary->IsSetId() ||
  96.         !standardDictionary->GetId().front().GetObject().IsOther_database() ||
  97.         standardDictionary->GetId().front().GetObject().GetOther_database().GetDb() != "Standard residue dictionary" ||
  98.         !standardDictionary->GetId().front().GetObject().GetOther_database().GetTag().IsId() ||
  99.         standardDictionary->GetId().front().GetObject().GetOther_database().GetTag().GetId() != 1)
  100.         FATALMSG("file '" << filename << "' does not contain expected dictionary data");
  101. }
  103. void DeleteStandardDictionary(void)
  104. {
  105.     if (standardDictionary) {
  106.         delete standardDictionary;
  107.         standardDictionary = NULL;
  108.     }
  109. }
  112. ChemicalGraph::ChemicalGraph(StructureBase *parent, const CBiostruc_graph& graph,
  113.     const FeatureList& features) :
  114.     StructureBase(parent), displayListOtherStart(OpenGLRenderer::NO_LIST)
  115. {
  116.     if (!standardDictionary) {
  117.         FATALMSG("need to load standard dictionary first");
  118.         return;
  119.     }
  121.     // figure out what models we'll be drawing, based on contents of parent
  122.     // StructureSet and StructureObject
  123.     const StructureObject *object;
  124.     if (!GetParentOfType(&object)) return;
  126.     // if this is the only StructureObject in this StructureSet, and if this
  127.     // StructureObject has only one CoordSet, and if this CoordSet's AtomSet
  128.     // has multiple ensembles, then use multiple altConf ensemble
  129.     int nAlts = 0;
  130.     if (!parentSet->IsMultiStructure() && object->coordSets.size() == 1 &&
  131.         object->coordSets.front()->atomSet->ensembles.size() > 1) {
  133.         AtomSet *atomSet = object->coordSets.front()->atomSet;
  134.         AtomSet::EnsembleList::iterator e, ee=atomSet->ensembles.end();
  135.         for (e=atomSet->ensembles.begin(); e!=ee; ++e) {
  136.             atomSetList.push_back(make_pair(atomSet, *e));
  137.             ++nAlts;
  138.         }
  140.     // otherwise, use all CoordSets using default altConf ensemble for single
  141.     // structure; for multiple structure StructureSet, use only first CoordSet
  142.     } else {
  143.         StructureObject::CoordSetList::const_iterator c, ce=object->coordSets.end();
  144.         for (c=object->coordSets.begin(); c!=ce; ++c) {
  145.             atomSetList.push_back(make_pair((*c)->atomSet,
  146.                 reinterpret_cast<const string *>(NULL)));   // VC++ requires this cast for some reason...
  147.             ++nAlts;
  148.             if (parentSet->IsMultiStructure()) break;
  149.         }
  150.     }
  151.     TRACEMSG("nAlts = " << nAlts);
  152.     if (!nAlts) {
  153.         WARNINGMSG("ChemicalGraph has zero AtomSets!");
  154.         return;
  155.     }
  157.     unsigned int firstNewFrame = parentSet->frameMap.size();
  158.     parentSet->frameMap.resize(firstNewFrame + nAlts);
  160.     // load molecules from SEQUENCE OF Molecule-graph
  161.     CBiostruc_graph::TMolecule_graphs::const_iterator i, ie=graph.GetMolecule_graphs().end();
  162.     for (i=graph.GetMolecule_graphs().begin(); i!=ie; ++i) {
  163.         Molecule *molecule = new Molecule(this,
  164.             i->GetObject(),
  165.             standardDictionary->GetResidue_graphs(),
  166.             graph.GetResidue_graphs());
  168.         if (molecules.find(molecule->id) != molecules.end())
  169.             ERRORMSG("confused by repeated Molecule-graph ID's");
  170.         molecules[molecule->id] = molecule;
  172.         // set molecules' display list; each protein or nucleotide molecule
  173.         // gets its own display list (one display list for each molecule for
  174.         // each set of coordinates), while everything else - hets, solvents,
  175.         // inter-molecule bonds - goes in a single list.
  176.         for (int n=0; n<nAlts; ++n) {
  178.             if (molecule->IsProtein() || molecule->IsNucleotide()) {
  180.                 molecule->displayLists.push_back(++(parentSet->lastDisplayList));
  181.                 parentSet->frameMap[firstNewFrame + n].push_back(parentSet->lastDisplayList);
  182.                 parentSet->transformMap[parentSet->lastDisplayList] = &(object->transformToMaster);
  184.             } else { // het/solvent
  185.                 if (displayListOtherStart == OpenGLRenderer::NO_LIST) {
  186.                     displayListOtherStart = parentSet->lastDisplayList + 1;
  187.                     parentSet->lastDisplayList += nAlts;
  188.                 }
  189.                 molecule->displayLists.push_back(displayListOtherStart + n);
  190.             }
  191.         }
  192.     }
  194.     // load connections from SEQUENCE OF Inter-residue-bond OPTIONAL
  195.     if (graph.IsSetInter_molecule_bonds()) {
  196.         CBiostruc_graph::TInter_molecule_bonds::const_iterator j, je=graph.GetInter_molecule_bonds().end();
  197.         for (j=graph.GetInter_molecule_bonds().begin(); j!=je; ++j) {
  199.             int order = j->GetObject().IsSetBond_order() ?
  200.                 j->GetObject().GetBond_order() : Bond::eUnknown;
  201.             const Bond *bond = MakeBond(this,
  202.                 j->GetObject().GetAtom_id_1(),
  203.                 j->GetObject().GetAtom_id_2(),
  204.                 order);
  205.             if (bond) interMoleculeBonds.push_back(bond);
  207.             if (CheckForDisulfide(NULL,
  208.                     j->GetObject().GetAtom_id_1(), j->GetObject().GetAtom_id_2(),
  209.                     &interMoleculeBonds, const_cast<Bond*>(bond), this) ||
  210.                 (bond && bond->order == Bond::eRealDisulfide)) {
  211.                 // set inter-molecule bonds' display list(s) if real bond or virtual disulfide created
  212.                 if (displayListOtherStart == OpenGLRenderer::NO_LIST) {
  213.                     displayListOtherStart = parentSet->lastDisplayList + 1;
  214.                     parentSet->lastDisplayList += nAlts;
  215.                 }
  216.             }
  217.         }
  218.     }
  220.     // if hets/solvent/i-m bonds present, add display lists to frames
  221.     if (displayListOtherStart != OpenGLRenderer::NO_LIST) {
  222.         for (int n=0; n<nAlts; ++n) {
  223.             parentSet->frameMap[firstNewFrame + n].push_back(displayListOtherStart + n);
  224.             parentSet->transformMap[displayListOtherStart + n] = &(object->transformToMaster);
  225.         }
  226.     }
  228.     // fill out secondary structure and domain maps from NCBI assigments (in feature block)
  229.     FeatureList::const_iterator l, le = features.end();
  230.     for (l=features.begin(); l!=le; ++l) {
  231.         if (l->GetObject().IsSetDescr()) {
  233.             // find and unpack NCBI sec. struc. or domain features
  234.             CBiostruc_feature_set::TDescr::const_iterator d, de = l->GetObject().GetDescr().end();
  235.             for (d=l->GetObject().GetDescr().begin(); d!=de; ++d) {
  236.                 if (d->GetObject().IsName()) {
  237.                     if (d->GetObject().GetName() == "NCBI assigned secondary structure")
  238.                         UnpackSecondaryStructureFeatures(l->GetObject());
  239.                     else if (d->GetObject().GetName() == "NCBI Domains")
  240.                         UnpackDomainFeatures(l->GetObject());
  241.                     break;
  242.                 }
  243.             }
  244.         }
  245.     }
  247.     // Assign a (new) domain ID to all residues in a chain, if there are not any domains
  248.     // already set; in other words, assume a chain with no given domain features is composed
  249.     // of a single domain.
  250.     list < int > moleculeIDs;
  251.     MoleculeMap::iterator m, me = molecules.end();
  252.     for (m=molecules.begin(); m!=me; ++m) {
  253.         if (m->second->IsProtein() || m->second->IsNucleotide())
  254.             moleculeIDs.push_back(m->first);
  255.     }
  256.     moleculeIDs.sort();    // assign in order of ID
  257.     list < int >::const_iterator id, ide = moleculeIDs.end();
  258.     for (id=moleculeIDs.begin(); id!=ide; ++id) {
  259.         int r;
  260.         Molecule *molecule = const_cast<Molecule*>(molecules[*id]);
  261.         for (r=0; r<molecule->residues.size(); ++r)
  262.             if (molecule->residueDomains[r] != Molecule::NO_DOMAIN_SET) break;
  263.         if (r == molecule->residues.size()) {
  264.             int domainID = ++((const_cast<StructureSet*>(parentSet))->nDomains);
  265.             ++(molecule->nDomains);
  266.             (const_cast<StructureObject*>(object))->domainMap[domainID] = molecule;
  267.             (const_cast<StructureObject*>(object))->domainID2MMDB[domainID] = -1;
  268.             for (r=0; r<molecule->residues.size(); ++r)
  269.                 molecule->residueDomains[r] = domainID;
  270.         }
  271.     }
  272. }
  274. void ChemicalGraph::UnpackDomainFeatures(const CBiostruc_feature_set& featureSet)
  275. {
  276.     TRACEMSG("unpacking NCBI domain features");
  277.     CBiostruc_feature_set::TFeatures::const_iterator f, fe = featureSet.GetFeatures().end();
  278.     for (f=featureSet.GetFeatures().begin(); f!=fe; ++f) {
  279.         if (f->GetObject().GetType() == CBiostruc_feature::eType_subgraph &&
  280.             f->GetObject().IsSetLocation() &&
  281.             f->GetObject().GetLocation().IsSubgraph() &&
  282.             f->GetObject().GetLocation().GetSubgraph().IsResidues() &&
  283.             f->GetObject().GetLocation().GetSubgraph().GetResidues().IsInterval()) {
  285.             // assign a "domain ID" successively (from 1) over the whole StructureSet
  286.             int domainID = ++((const_cast<StructureSet*>(parentSet))->nDomains);
  288.             // find molecule and set regions, warning about overlaps
  289.             Molecule *molecule = NULL;
  290.             CResidue_pntrs::TInterval::const_iterator i,
  291.                 ie = f->GetObject().GetLocation().GetSubgraph().GetResidues().GetInterval().end();
  292.             for (i=f->GetObject().GetLocation().GetSubgraph().GetResidues().GetInterval().begin(); i!=ie; ++i) {
  293.                 MoleculeMap::const_iterator m = molecules.find(i->GetObject().GetMolecule_id().Get());
  294.                 if (m == molecules.end() ||
  295.                     m->second->id !=    // check to make sure all intervals are on same molecule
  296.                         f->GetObject().GetLocation().GetSubgraph().GetResidues().GetInterval().
  297.                             front().GetObject().GetMolecule_id().Get()) {
  298.                     WARNINGMSG("Bad moleculeID in domain interval");
  299.                     continue;
  300.                 }
  301.                 molecule = const_cast<Molecule*>(m->second);
  302.                 for (int r=i->GetObject().GetFrom().Get()-1; r<=i->GetObject().GetTo().Get()-1; ++r) {
  303.                     if (r < 0 || r >= molecule->residues.size()) {
  304.                         ERRORMSG("Bad residue range in domain feature for moleculeID "
  305.                             << molecule->id << " residueID " << r+1);
  306.                         break;
  307.                     } else if (molecule->residueDomains[r] != Molecule::NO_DOMAIN_SET) {
  308.                         WARNINGMSG("Overlapping domain feature for moleculeID "
  309.                             << molecule->id << " residueID " << r+1);
  310.                         break;
  311.                     } else {
  312.                         molecule->residueDomains[r] = domainID;
  313.                     }
  314.                 }
  315.             }
  317.             if (molecule) {
  318.                 const StructureObject *object;
  319.                 if (!GetParentOfType(&object)) return;
  320.                 (const_cast<StructureObject*>(object))->domainMap[domainID] = molecule;
  321.                 (const_cast<StructureObject*>(object))->domainID2MMDB[domainID] =
  322.                     f->GetObject().GetId().Get();
  323.                 ++(molecule->nDomains);
  324.             }
  325.         }
  326.     }
  327. }
  329. void ChemicalGraph::UnpackSecondaryStructureFeatures(const CBiostruc_feature_set& featureSet)
  330. {
  331.     TRACEMSG("unpacking NCBI sec. struc. features");
  332.     CBiostruc_feature_set::TFeatures::const_iterator f, fe = featureSet.GetFeatures().end();
  333.     for (f=featureSet.GetFeatures().begin(); f!=fe; ++f) {
  334.         if ((f->GetObject().GetType() == CBiostruc_feature::eType_helix ||
  335.                 f->GetObject().GetType() == CBiostruc_feature::eType_strand) &&
  336.             f->GetObject().IsSetLocation() &&
  337.             f->GetObject().GetLocation().IsSubgraph() &&
  338.             f->GetObject().GetLocation().GetSubgraph().IsResidues() &&
  339.             f->GetObject().GetLocation().GetSubgraph().GetResidues().IsInterval()) {
  341.             // find molecule and set region, warning about overlaps
  342.             if (f->GetObject().GetLocation().GetSubgraph().GetResidues().GetInterval().size() > 1) {
  343.                 WARNINGMSG("Can't deal with multi-interval sec. struc. regions");
  344.                 continue;
  345.             }
  346.             const CResidue_interval_pntr& interval =
  347.                 f->GetObject().GetLocation().GetSubgraph().GetResidues().GetInterval().front().GetObject();
  348.             MoleculeMap::const_iterator m = molecules.find(interval.GetMolecule_id().Get());
  349.             if (m == molecules.end()) {
  350.                 WARNINGMSG("Bad moleculeID in sec. struc. interval");
  351.                 continue;
  352.             }
  353.             Molecule *molecule = const_cast<Molecule*>(m->second);
  354.             for (int r=interval.GetFrom().Get()-1; r<=interval.GetTo().Get()-1; ++r) {
  355.                 if (r < 0 || r >= molecule->residues.size()) {
  356.                     ERRORMSG("Bad residue range in sec. struc. feature for moleculeID "
  357.                         << molecule->id << " residueID " << r+1);
  358.                     break;
  359.                 } if (molecule->residueSecondaryStructures[r] != Molecule::eCoil) {
  360.                     WARNINGMSG("Overlapping sec. struc. feature at moleculeID "
  361.                         << molecule->id << " residueID " << r+1);
  362.                 } else {
  363.                     molecule->residueSecondaryStructures[r] =
  364.                         (f->GetObject().GetType() == CBiostruc_feature::eType_helix) ?
  365.                             Molecule::eHelix : Molecule::eStrand;
  366.                 }
  367.             }
  368.         }
  369.     }
  370. }
  372. static int moleculeToRedraw = -1;
  374. void ChemicalGraph::RedrawMolecule(int moleculeID) const
  375. {
  376.     moleculeToRedraw = moleculeID;
  377.     DrawAll(NULL);
  378.     moleculeToRedraw = -1;
  379. }
  381. // This is where the work of breaking objects up into display lists gets done.
  382. bool ChemicalGraph::DrawAll(const AtomSet *ignored) const
  383. {
  384.     const StructureObject *object;
  385.     if (!GetParentOfType(&object)) return false;
  387.     if (moleculeToRedraw != -1)
  388.         TRACEMSG("drawing molecule " << moleculeToRedraw
  389.             << " of ChemicalGraph of object " << object->pdbID);
  390.     else
  391.         TRACEMSG("drawing ChemicalGraph of object " << object->pdbID);
  393.     // put each protein (with its 3d-objects) or nucleotide chain in its own display list
  394.     bool continueDraw;
  395.     AtomSetList::const_iterator a, ae=atomSetList.end();
  396.     MoleculeMap::const_iterator m, me=molecules.end();
  397.     for (m=molecules.begin(); m!=me; ++m) {
  398.         if (!m->second->IsProtein() && !m->second->IsNucleotide()) continue;
  399.         if (moleculeToRedraw >= 0 && m->second->id != moleculeToRedraw) continue;
  401.         Molecule::DisplayListList::const_iterator md=m->second->displayLists.begin();
  402.         for (a=atomSetList.begin(); a!=ae; ++a, ++md) {
  404.             // start new display list
  405.             //TESTMSG("drawing molecule #" << m->second->id << " + its objects in display list " << *md
  406.             //        << " of " << atomSetList.size());
  407.             parentSet->renderer->StartDisplayList(*md);
  409.             // draw this molecule with all alternative AtomSets (e.g., NMR's or altConfs)
  410.             a->first->SetActiveEnsemble(a->second);
  411.             continueDraw = m->second->DrawAllWithTerminiLabels(a->first);
  413.             if (continueDraw) {
  414.                 // find 3D objects for this molecule/CoordSet
  415.                 const CoordSet *coordSet;
  416.                 if (a->first->GetParentOfType(&coordSet)) {
  417.                     CoordSet::Object3DMap::const_iterator objList = coordSet->objectMap.find(m->second->id);
  418.                     if (objList != coordSet->objectMap.end()) {
  419.                         CoordSet::Object3DList::const_iterator o, oe=objList->second.end();
  420.                         for (o=objList->second.begin(); o!=oe; ++o) {
  421.                             if (!(continueDraw = (*o)->Draw(a->first))) break;
  422.                         }
  423.                     }
  424.                 }
  425.             }
  427.             // end display list
  428.             parentSet->renderer->EndDisplayList();
  430.             if (!continueDraw) return false;
  431.         }
  433.         // we're done if this was the single molecule meant to be redrawn
  434.         if (moleculeToRedraw >= 0) break;
  435.     }
  437.     // then put everything else (solvents, hets, and intermolecule bonds) in a single display list
  438.     if (displayListOtherStart == OpenGLRenderer::NO_LIST) return true;
  439.     //TESTMSG("drawing hets/solvents/i-m bonds");
  441.     // always redraw all these even if only a single molecule is to be redrawn -
  442.     // that way connections can show/hide in cases where a particular residue
  443.     // changes its display
  444.     int n = 0;
  445.     for (a=atomSetList.begin(); a!=ae; ++a, ++n) {
  447.         a->first->SetActiveEnsemble(a->second);
  448.         parentSet->renderer->StartDisplayList(displayListOtherStart + n);
  450.         for (m=molecules.begin(); m!=me; ++m) {
  451.             if (m->second->IsProtein() || m->second->IsNucleotide()) continue;
  452.             if (!(continueDraw = m->second->DrawAll(a->first))) break;
  453.         }
  455.         if (continueDraw) {
  456.             BondList::const_iterator b, be=interMoleculeBonds.end();
  457.             for (b=interMoleculeBonds.begin(); b!=be; ++b) {
  458.                 if (!(continueDraw = (*b)->Draw(a->first))) break;
  459.             }
  460.         }
  462.         parentSet->renderer->EndDisplayList();
  464.         if (!continueDraw) return false;
  465.     }
  467.     return true;
  468. }
  470. bool ChemicalGraph::CheckForDisulfide(const Molecule *molecule,
  471.     const CAtom_pntr& atomPtr1, const CAtom_pntr& atomPtr2,
  472.     list < const Bond * > *bondList, Bond *bond, StructureBase *parent)
  473. {
  474.     if (atomPtr1.GetMolecule_id().Get() == atomPtr2.GetMolecule_id().Get() &&
  475.         atomPtr1.GetResidue_id().Get() == atomPtr2.GetResidue_id().Get()) return false;
  477.     const Molecule *mol1, *mol2;
  478.     if (molecule) { // when called from Molecule::Molecule() on interresidue (intramolecular) bond
  479.         mol1 = mol2 = molecule;
  480.     } else {        // when called from ChemicalGraph::ChemicalGraph() on intermolecule bond
  481.         MoleculeMap::const_iterator
  482.             m1 = molecules.find(atomPtr1.GetMolecule_id().Get()),
  483.             m2 = molecules.find(atomPtr2.GetMolecule_id().Get());
  484.         if (m1 == molecules.end() || m2 == molecules.end()) {
  485.             ERRORMSG("ChemicalGraph::CheckForDisulfide() - bad molecule ID");
  486.             return false;
  487.         }
  488.         mol1 = m1->second;
  489.         mol2 = m2->second;
  490.     }
  492.     Molecule::ResidueMap::const_iterator
  493.         res1 = mol1->residues.find(atomPtr1.GetResidue_id().Get()),
  494.         res2 = mol2->residues.find(atomPtr2.GetResidue_id().Get());
  495.     if (res1 == mol1->residues.end() || res2 == mol2->residues.end()) {
  496.         ERRORMSG("ChemicalGraph::CheckForDisulfide() - bad residue ID");
  497.         return false;
  498.     }
  500.     // check to make sure both residues are cysteine
  501.     if (res1->second->type != Residue::eAminoAcid || res1->second->code != 'C' ||
  502.         res2->second->type != Residue::eAminoAcid || res2->second->code != 'C') return false;
  504.     const Residue::AtomInfo
  505.         *atom1 = res1->second->GetAtomInfo(atomPtr1.GetAtom_id().Get()),
  506.         *atom2 = res2->second->GetAtomInfo(atomPtr2.GetAtom_id().Get());
  507.     if (!atom1 || !atom2) {
  508.         ERRORMSG("ChemicalGraph::CheckForDisulfide() - bad atom ID");
  509.         return false;
  510.     }
  512.     // check to make sure both atoms are sulfur, and that residue has an alpha
  513.     if (atom1->atomicNumber != 16 || res1->second->alphaID == Residue::NO_ALPHA_ID ||
  514.         atom2->atomicNumber != 16 || res2->second->alphaID == Residue::NO_ALPHA_ID) return false;
  516.     TRACEMSG("found disulfide between molecule " << atomPtr1.GetMolecule_id().Get()
  517.         << " residue " << atomPtr1.GetResidue_id().Get()
  518.         << " and molecule " << atomPtr2.GetMolecule_id().Get()
  519.         << " residue " << atomPtr2.GetResidue_id().Get());
  521.     // first flag this bond as "real disulfide", so it's not drawn with connection style
  522.     if (bond) bond->order = Bond::eRealDisulfide;
  524.     // then, make a new virtual disulfide bond between the alphas of these residues
  525.     const Bond *virtualDisulfide = MakeBond(parent,
  526.         atomPtr1.GetMolecule_id().Get(), atomPtr1.GetResidue_id().Get(), res1->second->alphaID,
  527.         atomPtr2.GetMolecule_id().Get(), atomPtr2.GetResidue_id().Get(), res2->second->alphaID,
  528.         Bond::eVirtualDisulfide);
  529.     if (!virtualDisulfide) return false;
  531.     bondList->push_back(virtualDisulfide);
  532.     return true;
  533. }
  535. END_SCOPE(Cn3D)
  538. /*
  539. * ---------------------------------------------------------------------------
  540. * $Log: chemical_graph.cpp,v $
  541. * Revision 1000.2  2004/06/01 18:28:03  gouriano
  542. * PRODUCTION: UPGRADED [GCC34_MSVC7] Dev-tree R1.40
  543. *
  544. * Revision 1.40  2004/05/21 21:41:39  gorelenk
  545. * Added PCH ncbi_pch.hpp
  546. *
  547. * Revision 1.39  2004/03/15 18:19:23  thiessen
  548. * prefer prefix vs. postfix ++/-- operators
  549. *
  550. * Revision 1.38  2004/02/19 17:04:46  thiessen
  551. * remove cn3d/ from include paths; add pragma to disable annoying msvc warning
  552. *
  553. * Revision 1.37  2003/08/21 17:56:29  thiessen
  554. * change header order for Mac compilation
  555. *
  556. * Revision 1.36  2003/06/21 21:04:41  thiessen
  557. * draw per-model 3D objects
  558. *
  559. * Revision 1.35  2003/03/06 19:23:18  thiessen
  560. * minor tweaks
  561. *
  562. * Revision 1.34  2003/02/03 19:20:02  thiessen
  563. * format changes: move CVS Log to bottom of file, remove std:: from .cpp files, and use new diagnostic macros
  564. *
  565. * Revision 1.33  2002/05/31 14:51:28  thiessen
  566. * small tweaks
  567. *
  568. * Revision 1.32  2002/03/18 15:17:29  thiessen
  569. * fix minor omission
  570. *
  571. * Revision 1.31  2001/12/12 14:04:13  thiessen
  572. * add missing object headers after object loader change
  573. *
  574. * Revision 1.30  2001/11/27 16:26:07  thiessen
  575. * major update to data management system
  576. *
  577. * Revision 1.29  2001/08/21 01:10:45  thiessen
  578. * add labeling
  579. *
  580. * Revision 1.28  2001/07/27 13:52:47  thiessen
  581. * make sure domains are assigned in order of molecule id; tweak pattern dialog
  582. *
  583. * Revision 1.27  2001/07/12 17:35:15  thiessen
  584. * change domain mapping ; add preliminary cdd annotation GUI
  585. *
  586. * Revision 1.26  2001/06/21 02:02:33  thiessen
  587. * major update to molecule identification and highlighting ; add toggle highlight (via alt)
  588. *
  589. * Revision 1.25  2001/05/31 18:47:07  thiessen
  590. * add preliminary style dialog; remove LIST_TYPE; add thread single and delete all; misc tweaks
  591. *
  592. * Revision 1.24  2001/05/18 22:58:52  thiessen
  593. * fix display list bug with disulfides
  594. *
  595. * Revision 1.23  2001/03/23 23:31:56  thiessen
  596. * keep atom info around even if coords not all present; mainly for disulfide parsing in virtual models
  597. *
  598. * Revision 1.22  2001/03/23 04:18:52  thiessen
  599. * parse and display disulfides
  600. *
  601. * Revision 1.21  2001/02/13 01:03:56  thiessen
  602. * backward-compatible domain ID's in output; add ability to delete rows
  603. *
  604. * Revision 1.20  2001/02/08 23:01:48  thiessen
  605. * hook up C-toolkit stuff for threading; working PSSM calculation
  606. *
  607. * Revision 1.19  2000/12/20 23:47:47  thiessen
  608. * load CDD's
  609. *
  610. * Revision 1.18  2000/12/19 16:39:08  thiessen
  611. * tweaks to show/hide
  612. *
  613. * Revision 1.17  2000/12/01 19:35:56  thiessen
  614. * better domain assignment; basic show/hide mechanism
  615. *
  616. * Revision 1.16  2000/11/30 15:49:36  thiessen
  617. * add show/hide rows; unpack sec. struc. and domain features
  618. *
  619. * Revision 1.15  2000/11/02 16:56:01  thiessen
  620. * working editor undo; dynamic slave transforms
  621. *
  622. * Revision 1.14  2000/09/14 14:55:34  thiessen
  623. * add row reordering; misc fixes
  624. *
  625. * Revision 1.13  2000/09/11 01:46:14  thiessen
  626. * working messenger for sequence<->structure window communication
  627. *
  628. * Revision 1.12  2000/08/27 18:52:20  thiessen
  629. * extract sequence information
  630. *
  631. * Revision 1.11  2000/08/21 17:22:37  thiessen
  632. * add primitive highlighting for testing
  633. *
  634. * Revision 1.10  2000/08/17 14:24:05  thiessen
  635. * added working StyleManager
  636. *
  637. * Revision 1.9  2000/08/16 14:18:44  thiessen
  638. * map 3-d objects to molecules
  639. *
  640. * Revision 1.8  2000/08/13 02:43:00  thiessen
  641. * added helix and strand objects
  642. *
  643. * Revision 1.7  2000/08/07 14:13:15  thiessen
  644. * added animation frames
  645. *
  646. * Revision 1.6  2000/08/07 00:21:17  thiessen
  647. * add display list mechanism
  648. *
  649. * Revision 1.5  2000/08/03 15:12:23  thiessen
  650. * add skeleton of style and show/hide managers
  651. *
  652. * Revision 1.4  2000/07/27 13:30:51  thiessen
  653. * remove 'using namespace ...' from all headers
  654. *
  655. * Revision 1.3  2000/07/16 23:19:10  thiessen
  656. * redo of drawing system
  657. *
  658. * Revision 1.2  2000/07/12 23:27:49  thiessen
  659. * now draws basic CPK model
  660. *
  661. * Revision 1.1  2000/07/11 13:45:29  thiessen
  662. * add modules to parse chemical graph; many improvements
  663. *
  664. */