structure_set.cpp

ncbi源码

                        parentSet->usedFeatures.end())
                continue;

            // look for alignment feature
            if (f2->GetObject().IsSetType() &&
				f2->GetObject().GetType() == CBiostruc_feature::eType_alignment &&
                f2->GetObject().IsSetLocation() &&
                f2->GetObject().GetLocation().IsAlignment()) {
                const CChem_graph_alignment& graphAlign =
					f2->GetObject().GetLocation().GetAlignment();

                // find transform alignment of this object with master
                if (graphAlign.GetDimension() == 2 &&
                    graphAlign.GetBiostruc_ids().size() == 2 &&
                    graphAlign.IsSetTransform() &&
                    graphAlign.GetTransform().size() == 1 &&
                    graphAlign.GetBiostruc_ids().front().GetObject().IsMmdb_id() &&
                    graphAlign.GetBiostruc_ids().front().GetObject().GetMmdb_id().Get() == masterMMDBID &&
                    graphAlign.GetBiostruc_ids().back().GetObject().IsMmdb_id() &&
                    graphAlign.GetBiostruc_ids().back().GetObject().GetMmdb_id().Get() == mmdbID) {

                    // mark feature as used
                    if (f2->GetObject().IsSetId())
                        parentSet->usedFeatures[f2->GetObject().GetId().Get()] = true;
                    TRACEMSG("got transform for " << pdbID << "->master");

                    // unpack transform into matrix, moves in reverse order;
                    Matrix xform;
                    transformToMaster = new Matrix();
                    CTransform::TMoves::const_iterator
                        m, me=graphAlign.GetTransform().front().GetObject().GetMoves().end();
                    for (m=graphAlign.GetTransform().front().GetObject().GetMoves().begin(); m!=me; ++m) {
                        Matrix xmat;
                        double scale;
                        if (m->GetObject().IsTranslate()) {
                            const CTrans_matrix& trans = m->GetObject().GetTranslate();
                            scale = 1.0 / trans.GetScale_factor();
                            SetTranslationMatrix(&xmat,
                                Vector(scale * trans.GetTran_1(),
                                       scale * trans.GetTran_2(),
                                       scale * trans.GetTran_3()));
                        } else { // rotate
                            const CRot_matrix& rot = m->GetObject().GetRotate();
                            scale = 1.0 / rot.GetScale_factor();
                            xmat.m[0]=scale*rot.GetRot_11(); xmat.m[1]=scale*rot.GetRot_21(); xmat.m[2]= scale*rot.GetRot_31(); xmat.m[3]=0;
                            xmat.m[4]=scale*rot.GetRot_12(); xmat.m[5]=scale*rot.GetRot_22(); xmat.m[6]= scale*rot.GetRot_32(); xmat.m[7]=0;
                            xmat.m[8]=scale*rot.GetRot_13(); xmat.m[9]=scale*rot.GetRot_23(); xmat.m[10]=scale*rot.GetRot_33(); xmat.m[11]=0;
                            xmat.m[12]=0;                    xmat.m[13]=0;                    xmat.m[14]=0;                     xmat.m[15]=1;
                        }
                        ComposeInto(transformToMaster, xmat, xform);
                        xform = *transformToMaster;
                    }
                    return true;
                }
            }
        }
    }

    WARNINGMSG("Can't get structure alignment for slave " << pdbID
        << " with master " << masterMMDBID << ";\nwill likely require manual realignment");
    return false;
}

static void AddDomain(int *domain, const Molecule *molecule, const Block::Range *range)
{
    const StructureObject *object;
    if (!molecule->GetParentOfType(&object)) return;
    for (int l=range->from; l<=range->to; ++l) {
        if (molecule->residueDomains[l] != Molecule::NO_DOMAIN_SET) {
            if (*domain == NO_DOMAIN)
                *domain = object->domainID2MMDB.find(molecule->residueDomains[l])->second;
            else if (*domain != MULTI_DOMAIN &&
                     *domain != object->domainID2MMDB.find(molecule->residueDomains[l])->second)
                *domain = MULTI_DOMAIN;
        }
    }
}

void StructureObject::RealignStructure(int nCoords,
    const Vector * const *masterCoords, const Vector * const *slaveCoords,
    const double *weights, int slaveRow)
{
    Vector masterCOM, slaveCOM; // centers of mass for master, slave
    Matrix slaveRotation;       // rotation to align slave with master
    if (!transformToMaster) transformToMaster = new Matrix();

    // do the fit
    RigidBodyFit(nCoords, masterCoords, slaveCoords, weights, masterCOM, slaveCOM, slaveRotation);

    // apply the resulting transform elements from the fit to this object's transform Matrix
    Matrix single, combined;
    SetTranslationMatrix(&single, -slaveCOM);
    ComposeInto(transformToMaster, slaveRotation, single);
    combined = *transformToMaster;
    SetTranslationMatrix(&single, masterCOM);
    ComposeInto(transformToMaster, single, combined);

    // print out RMSD
    Vector x;
    double rmsd = 0.0, d;
    int n = 0;
    for (int c=0; c<nCoords; ++c) {
        if (!slaveCoords[c] || !masterCoords[c]) continue;
        x = *(slaveCoords[c]);
        ApplyTransformation(&x, *transformToMaster);
        d = (*(masterCoords[c]) - x).length();
        rmsd += d * d;
        ++n;
    }
    rmsd = sqrt(rmsd / n);
    INFOMSG("RMSD of aligned alpha coordinates between master structure and " << pdbID << ": "
        << setprecision(3) << rmsd << setprecision(6) << " A");

    // create a new Biostruc-feature that contains this alignment
    CBiostruc_feature *feature = new CBiostruc_feature();
    feature->SetType(CBiostruc_feature::eType_alignment);
    CRef<CBiostruc_feature::C_Location> location(new CBiostruc_feature::C_Location());
    feature->SetLocation(*location);
    CChem_graph_alignment *graphAlignment = new CChem_graph_alignment();
    location.GetObject().SetAlignment(*graphAlignment);

    // fill out the Chem-graph-alignment
    graphAlignment->SetDimension(2);
    CMmdb_id
        *masterMID = new CMmdb_id(parentSet->objects.front()->mmdbID),
        *slaveMID = new CMmdb_id(mmdbID);
    CBiostruc_id
        *masterBID = new CBiostruc_id(),
        *slaveBID = new CBiostruc_id();
    masterBID->SetMmdb_id(*masterMID);
    slaveBID->SetMmdb_id(*slaveMID);
    graphAlignment->SetBiostruc_ids().resize(2);
    graphAlignment->SetBiostruc_ids().front().Reset(masterBID);
    graphAlignment->SetBiostruc_ids().back().Reset(slaveBID);
    graphAlignment->SetAlignment().resize(2);

    // fill out sequence alignment intervals, tracking domains in alignment
    const BlockMultipleAlignment *multiple = parentSet->alignmentManager->GetCurrentMultipleAlignment();
    int masterDomain = NO_DOMAIN, slaveDomain = NO_DOMAIN;
    const Molecule
        *masterMolecule = multiple->GetSequenceOfRow(0)->molecule,
        *slaveMolecule = multiple->GetSequenceOfRow(slaveRow)->molecule;
    BlockMultipleAlignment::UngappedAlignedBlockList blocks;
    multiple->GetUngappedAlignedBlocks(&blocks);
    if (blocks.size() > 0) {
        CChem_graph_pntrs
            *masterCGPs = new CChem_graph_pntrs(),
            *slaveCGPs = new CChem_graph_pntrs();
        graphAlignment->SetAlignment().front().Reset(masterCGPs);
        graphAlignment->SetAlignment().back().Reset(slaveCGPs);
        CResidue_pntrs
            *masterRPs = new CResidue_pntrs(),
            *slaveRPs = new CResidue_pntrs();
        masterCGPs->SetResidues(*masterRPs);
        slaveCGPs->SetResidues(*slaveRPs);

        masterRPs->SetInterval().resize(blocks.size());
        slaveRPs->SetInterval().resize(blocks.size());
        BlockMultipleAlignment::UngappedAlignedBlockList::const_iterator b, be = blocks.end();
        CResidue_pntrs::TInterval::iterator
            mi = masterRPs->SetInterval().begin(),
            si = slaveRPs->SetInterval().begin();
        for (b=blocks.begin(); b!=be; ++b, ++mi, ++si) {
            CResidue_interval_pntr
                *masterRIP = new CResidue_interval_pntr(),
                *slaveRIP = new CResidue_interval_pntr();
            mi->Reset(masterRIP);
            si->Reset(slaveRIP);

            masterRIP->SetMolecule_id().Set(masterMolecule->id);
            slaveRIP->SetMolecule_id().Set(slaveMolecule->id);

            const Block::Range *range = (*b)->GetRangeOfRow(0);
            masterRIP->SetFrom().Set(range->from + 1); // +1 to convert seqLoc to residueID
            masterRIP->SetTo().Set(range->to + 1);
            AddDomain(&masterDomain, masterMolecule, range);

            range = (*b)->GetRangeOfRow(slaveRow);
            slaveRIP->SetFrom().Set(range->from + 1);
            slaveRIP->SetTo().Set(range->to + 1);
            AddDomain(&slaveDomain, slaveMolecule, range);
        }
    }

    // fill out structure alignment transform
    CTransform *xform = new CTransform();
    graphAlignment->SetTransform().resize(1);
    graphAlignment->SetTransform().front().Reset(xform);
    xform->SetId(1);
    xform->SetMoves().resize(3);
    CTransform::TMoves::iterator m = xform->SetMoves().begin();
    for (int i=0; i<3; ++i, ++m) {
        CMove *move = new CMove();
        m->Reset(move);
        static const int scaleFactor = 100000;
        if (i == 0) {   // translate slave so its COM is at origin
            CTrans_matrix *trans = new CTrans_matrix();
            move->SetTranslate(*trans);
            trans->SetScale_factor(scaleFactor);
            trans->SetTran_1((int)(-(slaveCOM.x * scaleFactor)));
            trans->SetTran_2((int)(-(slaveCOM.y * scaleFactor)));
            trans->SetTran_3((int)(-(slaveCOM.z * scaleFactor)));
        } else if (i == 1) {
            CRot_matrix *rot = new CRot_matrix();
            move->SetRotate(*rot);
            rot->SetScale_factor(scaleFactor);
            rot->SetRot_11((int)(slaveRotation[0] * scaleFactor));
            rot->SetRot_12((int)(slaveRotation[4] * scaleFactor));
            rot->SetRot_13((int)(slaveRotation[8] * scaleFactor));
            rot->SetRot_21((int)(slaveRotation[1] * scaleFactor));
            rot->SetRot_22((int)(slaveRotation[5] * scaleFactor));
            rot->SetRot_23((int)(slaveRotation[9] * scaleFactor));
            rot->SetRot_31((int)(slaveRotation[2] * scaleFactor));
            rot->SetRot_32((int)(slaveRotation[6] * scaleFactor));
            rot->SetRot_33((int)(slaveRotation[10] * scaleFactor));
        } else if (i == 2) {    // translate slave so its COM is at COM of master
            CTrans_matrix *trans = new CTrans_matrix();
            move->SetTranslate(*trans);
            trans->SetScale_factor(scaleFactor);
            trans->SetTran_1((int)(masterCOM.x * scaleFactor));
            trans->SetTran_2((int)(masterCOM.y * scaleFactor));
            trans->SetTran_3((int)(masterCOM.z * scaleFactor));
        }
    }

    // store the new alignment in the Biostruc-annot-set,
    // setting the feature id depending on the aligned domain(s)
    if (masterDomain == NO_DOMAIN) masterDomain = 0;    // can happen if single domain chain
    if (slaveDomain == NO_DOMAIN) slaveDomain = 0;
    const StructureObject *masterObject;
    if (!masterMolecule->GetParentOfType(&masterObject)) return;
    int
        masterDomainID = masterObject->mmdbID*10000 + masterMolecule->id*100 + masterDomain,
        slaveDomainID = mmdbID*100000 + slaveMolecule->id*1000 + slaveDomain*10 + 1;
    parentSet->AddStructureAlignment(feature, masterDomainID, slaveDomainID);

    // for backward-compatibility with Cn3D 3.5, need name to encode chain/domain
    CNcbiOstrstream oss;
    oss << masterMolecule->identifier->pdbID << ((char) masterMolecule->identifier->pdbChain) << masterDomain << ' '
        << slaveMolecule->identifier->pdbID << ((char) slaveMolecule->identifier->pdbChain) << slaveDomain << ' '
        << "Structure alignment of slave " << multiple->GetSequenceOfRow(slaveRow)->identifier->ToString()
        << " with master " << multiple->GetSequenceOfRow(0)->identifier->ToString()
        << ", as computed by Cn3D" << '\0';
    feature->SetName(string(oss.str()));
    delete oss.str();
}

void StructureObject::SelectByDistance(double cutoff, bool biopolymersOnly, bool otherMoleculesOnly,
    ResidueMap *selectedResidues) const
{
    // first make a list of coordinates of atoms in selected residues,
    typedef vector < const AtomCoord * > CoordList;
    CoordList highlightedAtoms;
    typedef map < const Molecule * , bool > MoleculeList;
    MoleculeList moleculesWithHighlights;

    ChemicalGraph::MoleculeMap::const_iterator m, me = graph->molecules.end();
    for (m=graph->molecules.begin(); m!=me; ++m) {
        Molecule::ResidueMap::const_iterator r, re = m->second->residues.end();
        for (r=m->second->residues.begin(); r!=re; ++r) {

            if (GlobalMessenger()->IsHighlighted(m->second, r->second->id)) {
                const Residue::AtomInfoMap& atomInfos = r->second->GetAtomInfos();
                Residue::AtomInfoMap::const_iterator a, ae = atomInfos.end();
                for (a=atomInfos.begin(); a!=ae; ++a) {
                    const AtomCoord *atomCoord = coordSets.front()->atomSet->
                        GetAtom(AtomPntr(m->second->id, r->second->id, a->first), true, true);
                    if (atomCoord) highlightedAtoms.push_back(atomCoord);
                }
                moleculesWithHighlights[m->second] = true;
            }
        }
    }
    if (highlightedAtoms.size() == 0) return;

    // now make a list of unselected residues to check for proximity
    typedef vector < const Residue * > ResidueList;
    ResidueList unhighlightedResidues;

    for (m=graph->molecules.begin(); m!=me; ++m) {
        Molecule::ResidueMap::const_iterator r, re = m->second->residues.end();

        if (otherMoleculesOnly &&
                moleculesWithHighlights.find(m->second) != moleculesWithHighlights.end())
            continue;

        for (r=m->second->residues.begin(); r!=re; ++r) {

            if (!GlobalMessenger()->IsHighlighted(m->second, r->second->id) &&
                    (!biopolymersOnly || m->second->IsProtein() || m->second->IsNucleotide()))
                unhighlightedResidues.push_back(r->second);
        }
    }
    if (unhighlightedResidues.size() == 0) return;

    // now check all unhighlighted residues, to see if any atoms are within cutoff distance
    // of any highlighted atoms; if so, add to residue selection list
    CoordList::const_iterator h, he = highlightedAtoms.end();
    ResidueLis