molecule.cpp

ncbi源码

/*
 * ===========================================================================
 * PRODUCTION $Log: molecule.cpp,v $
 * PRODUCTION Revision 1000.2  2004/06/01 18:28:41  gouriano
 * PRODUCTION PRODUCTION: UPGRADED [GCC34_MSVC7] Dev-tree R1.43
 * PRODUCTION
 * ===========================================================================
 */

/*  $Id: molecule.cpp,v 1000.2 2004/06/01 18:28:41 gouriano Exp $
* ===========================================================================
*
*                            PUBLIC DOMAIN NOTICE
*               National Center for Biotechnology Information
*
*  This software/database is a "United States Government Work" under the
*  terms of the United States Copyright Act.  It was written as part of
*  the author's official duties as a United States Government employee and
*  thus cannot be copyrighted.  This software/database is freely available
*  to the public for use. The National Library of Medicine and the U.S.
*  Government have not placed any restriction on its use or reproduction.
*
*  Although all reasonable efforts have been taken to ensure the accuracy
*  and reliability of the software and data, the NLM and the U.S.
*  Government do not and cannot warrant the performance or results that
*  may be obtained by using this software or data. The NLM and the U.S.
*  Government disclaim all warranties, express or implied, including
*  warranties of performance, merchantability or fitness for any particular
*  purpose.
*
*  Please cite the author in any work or product based on this material.
*
* ===========================================================================
*
* Authors:  Paul Thiessen
*
* File Description:
*      Classes to hold molecules
*
* ===========================================================================
*/

#ifdef _MSC_VER
#pragma warning(disable:4018)   // disable signed/unsigned mismatch warning in MSVC
#endif

#include <ncbi_pch.hpp>
#include <corelib/ncbistd.hpp>

#include <objects/mmdb1/Residue_graph.hpp>
#include <objects/mmdb1/Molecule_id.hpp>
#include <objects/seqloc/Seq_id.hpp>
#include <objects/seqloc/PDB_seq_id.hpp>
#include <objects/general/Object_id.hpp>
#include <objects/mmdb1/Residue_id.hpp>
#include <objects/seqloc/PDB_mol_id.hpp>

#include "molecule.hpp"
#include "residue.hpp"
#include "bond.hpp"
#include "style_manager.hpp"
#include "structure_set.hpp"
#include "coord_set.hpp"
#include "atom_set.hpp"
#include "chemical_graph.hpp"
#include "molecule_identifier.hpp"
#include "show_hide_manager.hpp"
#include "cn3d_tools.hpp"
#include "opengl_renderer.hpp"
#include "cn3d_colors.hpp"

USING_NCBI_SCOPE;
USING_SCOPE(objects);


BEGIN_SCOPE(Cn3D)

const int Molecule::NO_DOMAIN_SET = -1;

Molecule::Molecule(ChemicalGraph *parentGraph,
    const CMolecule_graph& graph,
    const ResidueGraphList& standardDictionary,
    const ResidueGraphList& localDictionary) :
    StructureBase(parentGraph), type(eOther), sequence(NULL), nDomains(0)
{
    int gi = MoleculeIdentifier::VALUE_NOT_SET, pdbChain = MoleculeIdentifier::VALUE_NOT_SET;
    string pdbID, accession;

    // get ID, name, and type
    id = graph.GetId().Get();
    CMolecule_graph::TDescr::const_iterator k, ke=graph.GetDescr().end();
    for (k=graph.GetDescr().begin(); k!=ke; ++k) {
        if (k->GetObject().IsMolecule_type()) {
            type = static_cast<eType>(k->GetObject().GetMolecule_type());
        }
        if (k->GetObject().IsName()) {
            name = k->GetObject().GetName();
        }
    }

    // get Seq-id for biopolymer chains (for sequence alignment stuff)
    if (IsProtein() || IsNucleotide()) {
        if (graph.IsSetSeq_id()) {
            if (graph.GetSeq_id().IsGi())
                gi = graph.GetSeq_id().GetGi();
            else if (graph.GetSeq_id().IsPdb()) {
                pdbID = graph.GetSeq_id().GetPdb().GetMol().Get();
                if (graph.GetSeq_id().GetPdb().IsSetChain())
                    pdbChain = graph.GetSeq_id().GetPdb().GetChain();
                else
                    pdbChain = ' ';
            }
            else if (graph.GetSeq_id().IsLocal() && graph.GetSeq_id().GetLocal().IsStr()) {
                accession = graph.GetSeq_id().GetLocal().GetStr();
                // special case where local accession is actually a PDB chain + extra stuff
                if (pdbID.size() == 0 && accession.size() >= 7 &&
                        accession[4] == ' ' && accession[6] == ' ' && isalpha(accession[5])) {
                    pdbID = accession.substr(0, 4);
                    pdbChain = accession[5];
                    accession.erase();
                }
            }
        }
        if (gi == MoleculeIdentifier::VALUE_NOT_SET && pdbID.size() == 0 && accession.size() == 0) {
            ERRORMSG("Molecule::Molecule() - biopolymer molecule, but can't get Seq-id");
            return;
        }
    }

    // if no PDB id assigned to biopolymer, assign default PDB identifier from parent, assuming 'name'
    // is actually the chainID if this is a biopolymer
    const StructureObject *object;
    if (!GetParentOfType(&object)) return;
    if ((IsProtein() || IsNucleotide()) && pdbID.size() == 0) {
        pdbID = object->pdbID;
        if (name.size() == 1) pdbChain = name[0];
    }

    // load residues from SEQUENCE OF Residue, storing virtual bonds along the way
    CMolecule_graph::TInter_residue_bonds::const_iterator j, je, jOrig;
    if (graph.IsSetInter_residue_bonds()) {
        j = graph.GetInter_residue_bonds().begin();
        je = graph.GetInter_residue_bonds().end();
    }
    const Residue *prevResidue = NULL;
    const Bond *prevBond = NULL;
    CMolecule_graph::TResidue_sequence::const_iterator i, ie=graph.GetResidue_sequence().end();
    int nResidues = 0;

    for (i=graph.GetResidue_sequence().begin(); i!=ie; ++i) {

        const Residue *residue = new Residue(this, (*i).GetObject(), id,
            standardDictionary, localDictionary);
        residues[residue->id] = residue;
        ++nResidues;

        // this assumption is frequently made elsewhere, relating seqLocs (numbering from zero)
        // to residue ID's (numbering from one) - so enforce it here
        if (residue->id != nResidues)
            ERRORMSG("Residue ID's must be ordered consecutively starting with one");

        // virtual bonds between successive alphas
        if (prevResidue && prevResidue->alphaID != Residue::NO_ALPHA_ID &&
            residue->alphaID != Residue::NO_ALPHA_ID) {

            bool foundReal = false;

            // see if there's a "real" inter-residue bond between these in the chemical graph
            if (graph.IsSetInter_residue_bonds()) {
                jOrig = j;
                do {
                    if ((j->GetObject().GetAtom_id_1().GetResidue_id().Get() == prevResidue->id &&
                         j->GetObject().GetAtom_id_2().GetResidue_id().Get() == residue->id) ||
                        (j->GetObject().GetAtom_id_2().GetResidue_id().Get() == prevResidue->id &&
                         j->GetObject().GetAtom_id_1().GetResidue_id().Get() == residue->id))
                        foundReal = true;
                    if (++j == je) j = graph.GetInter_residue_bonds().begin();
                } while (!foundReal && j != jOrig);
            }

            // for alpha only models, there are no inter-residue bonds in the
            // chemical graph, so check inter-atomic distances for cases where either one
            // of the two residues is alpha only.
            if (!foundReal && (IsProtein() || IsNucleotide()) &&
                (residue->NAtomsWithAnyCoords() == 1 || prevResidue->NAtomsWithAnyCoords() == 1)) {

                // get atom coordinates
                AtomPntr ap1(id, residue->id, residue->alphaID);
                const AtomCoord* atom1 = object->coordSets.front()->atomSet->GetAtom(ap1, true, true);
                AtomPntr ap2(id, prevResidue->id, prevResidue->alphaID);
                const AtomCoord* atom2 = object->coordSets.front()->atomSet->GetAtom(ap2, true, true);

                // check distance - ok if <= 5.0 Angstroms (or 10.0 for nucleotides)
                if (atom1 && atom2 &&
                        (atom1->site - atom2->site).length() <=
                            (IsProtein() ? 5.0 : 10.0))
                    foundReal = true;
            }

            if (foundReal) {
                const Bond *bond = MakeBond(this,
                    id, prevResidue->id, prevResidue->alphaID,
                    id, residue->id, residue->alphaID,
                    Bond::eVirtual);
                if (bond) {
                    interResidueBonds.push_back(bond);
                    if (prevBond) {
                        (const_cast<Bond *>(prevBond))->nextVirtual = bond;
                        (const_cast<Bond *>(bond))->previousVirtual = prevBond;
                    }
                }
                prevBond = bond;
            } else
                prevBond = NULL;
        } else
            prevBond = NULL;

        prevResidue = residue;
    }

    residueDomains.resize(residues.size(), NO_DOMAIN_SET);

    // keep s.s. maps only for protein chains
    if (IsProtein()) residueSecondaryStructures.resize(residues.size(), eCoil);

    // load inter-residue bonds from SEQUENCE OF Inter-residue-bond OPTIONAL
    if (graph.IsSetInter_residue_bonds()) {
        je = graph.GetInter_residue_bonds().end();
        for (j=graph.GetInter_residue_bonds().begin(); j!=je; ++j) {

            int order = j->GetObject().IsSetBond_order() ?
                j->GetObject().GetBond_order() : Bond::eUnknown;
            const Bond *bond = MakeBond(this,
                j->GetObject().GetAtom_id_1(),
                j->GetObject().GetAtom_id_2(),
                order);
            if (bond) interResidueBonds.push_back(bond);

            // add to disulfide map if virtual disulfide added or this bond is flagged as disulfide
            if (parentGraph->CheckForDisulfide(this,
                    j->GetObject().GetAtom_id_1(), j->GetObject().GetAtom_id_2(),
                    &interResidueBonds, const_cast<Bond*>(bond), this) ||
                (bond && bond->order == Bond::eRealDisulfide)) {
                disulfideMap[j->GetObject().GetAtom_id_1().GetResidue_id().Get()] =
                    j->GetObject().GetAtom_id_2().GetResidue_id().Get();
                disulfideMap[j->GetObject().GetAtom_id_2().GetResidue_id().Get()] =
                    j->GetObject().GetAtom_id_1().GetResidue_id().Get();
            }
        }
    }

    // get identifier
    identifier = MoleculeIdentifier::GetIdentifier(this, pdbID, pdbChain, gi, accession);
}

Vector Molecule::GetResidueColor(int sequenceIndex) const
{
    static const Vector gray(.5,.5,.5);

    // this assumes that the "index" - the position of the residue in the sequence,
    // starting from zero, is always one less than the residueID from the ASN1
    // data, which starts from one.
    ResidueMap::const_iterator r = residues.find(sequenceIndex + 1);

    if (r == residues.end()) return gray;
    const Residue *residue = r->second;

    // if no known alpha atom, just use gray
    if (residue->alphaID == Residue::NO_ALPHA_ID) return gray;
    AtomPntr atom(id, residue->id, residue->alphaID);

    // just use the first AtomSet (e.g. first model) of this object
    const StructureObject *object;
    if (!GetParentOfType(&object)) return gray;
    const AtomCoord *atomCoord = object->coordSets.front()->atomSet->