structure_set.cpp

ncbi源码

{
    // try to determine the master structure
    int masterMMDBID = MoleculeIdentifier::VALUE_NOT_SET;
    // explicit master structure
    if (dataManager->GetMasterStructure() &&
        dataManager->GetMasterStructure()->GetId().front()->IsMmdb_id())
        masterMMDBID = dataManager->GetMasterStructure()->GetId().front()->GetMmdb_id().Get();
    // master of structure alignments
    else if (dataManager->GetStructureAlignments() &&
                dataManager->GetStructureAlignments()->IsSetId() &&
                dataManager->GetStructureAlignments()->GetId().front()->IsMmdb_id())
        masterMMDBID = dataManager->GetStructureAlignments()->GetId().front()->GetMmdb_id().Get();
    // SPECIAL CASE: if there's no explicit master structure, but there is a structure
    // list that contains a single structure, then assume that it is the master structure
    else if (dataManager->GetStructureList() &&
                dataManager->GetStructureList()->size() == 1 &&
                dataManager->GetStructureList()->front()->GetId().front()->IsMmdb_id())
        masterMMDBID = dataManager->GetStructureList()->front()->GetId().front()->GetMmdb_id().Get();

    // assume data manager has already screened the alignment list
    ASNDataManager::SeqAnnotList *alignments = dataManager->GetSequenceAlignments();
    typedef list < CRef < CSeq_align > > SeqAlignList;
    const SeqAlignList& seqAligns = alignments->front()->GetData().GetAlign();

    // we need to determine the identity of the master sequence; most rigorous way is to look
    // for a Seq-id that is present in all pairwise alignments
    const Sequence *seq1 = NULL, *seq2 = NULL, *master = NULL;
    bool seq1PresentInAll = true, seq2PresentInAll = true;

    // first, find sequences for first pairwise alignment
    const CSeq_id& frontSid = seqAligns.front()->GetSegs().IsDendiag() ?
        seqAligns.front()->GetSegs().GetDendiag().front()->GetIds().front().GetObject() :
        seqAligns.front()->GetSegs().GetDenseg().GetIds().front().GetObject();
    const CSeq_id& backSid = seqAligns.front()->GetSegs().IsDendiag() ?
        seqAligns.front()->GetSegs().GetDendiag().front()->GetIds().back().GetObject() :
        seqAligns.front()->GetSegs().GetDenseg().GetIds().back().GetObject();
    SequenceSet::SequenceList::const_iterator s, se = sequenceSet->sequences.end();
    for (s=sequenceSet->sequences.begin(); s!=se; ++s) {
        if ((*s)->identifier->MatchesSeqId(frontSid)) seq1 = *s;
        if ((*s)->identifier->MatchesSeqId(backSid)) seq2 = *s;
        if (seq1 && seq2) break;
    }
    if (!(seq1 && seq2)) {
        ERRORMSG("Can't match first pair of Seq-ids to Sequences");
        return;
    }

    // now, make sure one of these sequences is present in all the other pairwise alignments
    SeqAlignList::const_iterator a = seqAligns.begin(), ae = seqAligns.end();
    for (++a; a!=ae; ++a) {
        const CSeq_id& frontSid2 = (*a)->GetSegs().IsDendiag() ?
            (*a)->GetSegs().GetDendiag().front()->GetIds().front().GetObject() :
            (*a)->GetSegs().GetDenseg().GetIds().front().GetObject();
        const CSeq_id& backSid2 = (*a)->GetSegs().IsDendiag() ?
            (*a)->GetSegs().GetDendiag().front()->GetIds().back().GetObject() :
            (*a)->GetSegs().GetDenseg().GetIds().back().GetObject();
        if (!seq1->identifier->MatchesSeqId(frontSid2) && !seq1->identifier->MatchesSeqId(backSid2))
            seq1PresentInAll = false;
        if (!seq2->identifier->MatchesSeqId(frontSid2) && !seq2->identifier->MatchesSeqId(backSid2))
            seq2PresentInAll = false;
    }
    if (!seq1PresentInAll && !seq2PresentInAll) {
        ERRORMSG("All pairwise sequence alignments must have a common master sequence");
        return;
    } else if (seq1PresentInAll && !seq2PresentInAll)
        master = seq1;
    else if (seq2PresentInAll && !seq1PresentInAll)
        master = seq2;
    else if (seq1PresentInAll && seq2PresentInAll && seq1 == seq2)
        master = seq1;

    // if still ambiguous, see if master3d is set in CDD data
    if (!master && dataManager->GetCDDMaster3d()) {
        if (seq1->identifier->MatchesSeqId(*(dataManager->GetCDDMaster3d())))
            master = seq1;
        else if (seq2->identifier->MatchesSeqId(*(dataManager->GetCDDMaster3d())))
            master = seq2;
        else
            ERRORMSG("Unable to match CDD's master3d with either sequence in first pairwise alignment");
    }

    // if still ambiguous, try to use master structure info to find master sequence
    if (!master && masterMMDBID != MoleculeIdentifier::VALUE_NOT_SET) {
        // load master - has side affect of matching gi's with PDB/molecule ID during graph evaluation
        if (structureLimit > 0)
            LoadMaster(masterMMDBID);

        // see if there's a sequence in the master structure that matches
        if (seq1->identifier->mmdbID != MoleculeIdentifier::VALUE_NOT_SET &&
            seq1->identifier->mmdbID != seq2->identifier->mmdbID) {
            if (masterMMDBID == seq1->identifier->mmdbID)
                master = seq1;
            else if (masterMMDBID == seq2->identifier->mmdbID)
                master = seq2;
            else {
                ERRORMSG("Structure master does not contain either sequence in first pairwise alignment");
                return;
            }
        }
    }

    // if still ambiguous, just use the first one
    if (!master) {
        WARNINGMSG("Ambiguous master; using " << seq1->identifier->ToString());
        master = seq1;
    }

    TRACEMSG("determined that master sequence is " << master->identifier->ToString());

    // load alignments now that we know the identity of the master
    alignmentSet = new AlignmentSet(this, master, *(dataManager->GetSequenceAlignments()));

    // check mmdb id's and load master if not already present (and if master has structure)
    if (masterMMDBID == MoleculeIdentifier::VALUE_NOT_SET) {
        masterMMDBID = master->identifier->mmdbID;
    } else if (master->identifier->mmdbID != MoleculeIdentifier::VALUE_NOT_SET &&
               master->identifier->mmdbID != masterMMDBID) {
        ERRORMSG("master structure (" << masterMMDBID <<
            ") disagrees with master sequence (apparently from " << master->identifier->mmdbID << ')');
        return;
    }
    if (objects.size() == 0 && structureLimit > 0 && masterMMDBID != MoleculeIdentifier::VALUE_NOT_SET &&
            master->identifier->accession != "consensus")   // special case for looking at "raw" CD's
        LoadMaster(masterMMDBID);

    // cross-match master sequence and structure
    if (objects.size() == 1 && !MatchSequenceToMoleculeInObject(master, objects.front())) {
        ERRORMSG("MatchSequenceToMoleculeInObject() - can't find molecule in object "
            << objects.front()->pdbID << " to match master sequence "
            << master->identifier->ToString());
        return;
    }

    // IFF there's a master structure, then also load slave structures and cross-match sequences
    if (objects.size() == 1 && structureLimit > 1) {
        ASNDataManager::BiostrucList::const_iterator b, be;
        if (dataManager->GetStructureList()) be = dataManager->GetStructureList()->end();
        int row;
        vector < bool > loadedStructureForSlaveRow(alignmentSet->alignments.size(), false);

        // first, load each remaining slave structure, and for each one, find the first slave
        // sequence that matches it (and that doesn't already have structure)
        AlignmentSet::AlignmentList::const_iterator l, le = alignmentSet->alignments.end();
        if (dataManager->GetStructureList()) {
            for (b=dataManager->GetStructureList()->begin(); b!=be && objects.size()<structureLimit; ++b) {

                // load structure
                if (usedStructures.find(b->GetPointer()) != usedStructures.end()) continue;
                StructureObject *object = new StructureObject(this, **b, false);
                objects.push_back(object);
                if (dataManager->GetStructureAlignments())
                    object->SetTransformToMaster(
                        *(dataManager->GetStructureAlignments()), master->identifier->mmdbID);
                usedStructures[b->GetPointer()] = true;

                // find matching unstructured slave sequence
                for (l=alignmentSet->alignments.begin(), row=0; l!=le; ++l, ++row) {
                    if (loadedStructureForSlaveRow[row]) continue;
                    if (MatchSequenceToMoleculeInObject((*l)->slave, object,
                            &((const_cast<MasterSlaveAlignment*>(*l))->slave))) {
                        loadedStructureForSlaveRow[row] = true;
                        break;
                    }
                }
                if (l == le)
                    ERRORMSG("Warning: Structure " << object->pdbID
                        << " doesn't have a matching slave sequence in the multiple alignment");
            }
        }

        // now loop through slave rows of the alignment; if the slave
        // sequence has an MMDB ID but no structure yet, then load it.
        if (objects.size() < structureLimit && (dataManager->IsCDD() || dataManager->IsGeneralMime())) {

            // for CDD's, ask user which structures to load if structureLimit is low
            if (dataManager->IsCDD())
                SetStructureRowFlags(alignmentSet, &structureLimit, &loadedStructureForSlaveRow);

            for (l=alignmentSet->alignments.begin(), row=0; l!=le && objects.size()<structureLimit; ++l, ++row) {

                if ((*l)->slave->identifier->mmdbID != MoleculeIdentifier::VALUE_NOT_SET &&
                    !loadedStructureForSlaveRow[row]) {

                    // first check the biostruc list to see if this structure is present already
                    CRef < CBiostruc > biostruc;
                    if (dataManager->GetStructureList()) {
                        for (b=dataManager->GetStructureList()->begin(); b!=be ; ++b) {
                            if ((*b)->GetId().front()->IsMmdb_id() &&
                                (*b)->GetId().front()->GetMmdb_id().Get() == (*l)->slave->identifier->mmdbID) {
                                biostruc = *b;
                                break;
                            }
                        }
                    }

                    // if not in list, load Biostruc via HTTP/cache
                    if (biostruc.Empty()) {
                        wxString id;
                        id.Printf("%i", (*l)->slave->identifier->mmdbID);
                        if (!LoadStructureViaCache(id.c_str(),
                                dataManager->GetBiostrucModelType(), biostruc, NULL)) {
                            ERRORMSG("Failed to load MMDB #" << (*l)->slave->identifier->mmdbID);
                            continue;
                        }
                    }

                    // create StructureObject and cross-match
                    StructureObject *object = new StructureObject(this, *biostruc, false);
                    objects.push_back(object);
                    if (dataManager->GetStructureAlignments())
                        object->SetTransformToMaster(
                            *(dataManager->GetStructureAlignments()), master->identifier->mmdbID);
                    if (!MatchSequenceToMoleculeInObject((*l)->slave, object,
                            &((const_cast<MasterSlaveAlignment*>(*l))->slave)))
                        ERRORMSG("Failed to match any molecule in structure " << object->pdbID
                            << " with sequence " << (*l)->slave->identifier->ToString());
                    loadedStructureForSlaveRow[row] = true;
                }
            }
        }
    }
}

void StructureSet::Load(int structureLimit)
{
    // member data initialization
    lastAtomName = OpenGLRenderer::NO_NAME;
    lastDisplayList = OpenGLRenderer::NO_LIST;
    sequenceSet = NULL;
    alignmentSet = NULL;
    alignmentManager = NULL;
    nDomains = 0;
    isAlphaOnly = false;
    parentSet = this;
    showHideManager = new ShowHideManager();
    styleManager = new StyleManager(this);
    havePrevPickedAtomCoord = false;
    hasUserStyle = false;

    // if this is a single structure, then there should be one sequence per biopolymer
    if (dataManager->IsSingleStructure()) {
        const CBiostruc *masterBiostruc = dataManager->GetMasterStructure();
        if (!masterBiostruc && dataManager->GetStructureList() && dataManager->GetStructureList()->size() == 1)
            masterBiostruc = dataManager->GetStructureList()->front().GetPointer();
        if (masterBiostruc)
            objects.push_back(new StructureObject(this, *masterBiostruc, true));
        if (dataManager->GetSequences())
            LoadSequencesForSingleStructure();
    }

    // multiple structure: should have exactly one sequence per structure (plus unstructured sequences)
    else {
        if (!dataManager->GetSequences() || !dataManager->GetSequenceAlignments()) {
            ERRORMSG("Data interpreted as multiple alignment, "
                "but missing sequences and/or sequence alignments");
            return;
        }
        sequenceSet = new SequenceSet(this, *(dataManager->GetSequences()));
        LoadAlignmentsAndStructures(structureLimit);
    }

    // find center of coordinates
    SetCenter();

    // create alignment manager
    if (sequenceSet) {
        if (dataManager->GetUpdates())
            // if updates present, alignment manager will load those into update viewer
            alignmentManager = new AlignmentManager(sequenceSet, alignmentSet, *(dataManager->GetUpdates()));
        else
            alignmentManager = new AlignmentManager(sequenceSet, alignmentSet);
    }

#ifdef _DEBUG
    VerifyFrameMap();
#endif

    // setup show/hide items
    showHideManager->ConstructShowHideArray(this);
    // alignments always start with aligned domains only
    if (alignmentSet)
        showHideManager->ShowAlignedDomains(this);

    // load style dictionary and user annotations
    const CCn3d_style_dictionary *styles = dataManager->GetStyleDictionary();
    if (styles) {
        if (!styleManager->LoadFromASNStyleDictionary(*styles) ||
            !styleManager->CheckGlobalStyleSettings())
            ERRORMSG("Error loading style dictionary");
        dataManager->RemoveStyleDictionary();   // remove now; recreated with current settings upon save
        hasUserStyle = true;
    }

    const CCn3d_user_annotations *annots = dataManager->GetUserAnnotations();
    if (annots) {
        if (!styleManager->LoadFromASNUserAnnotations(*annots) ||
            !renderer->LoadFromASNViewSettings(*annots))
            ERRORMSG("Error loading user annotations or camera settings");
        dataManager->RemoveUserAnnotations();   // remove now; recreated with current settings upon save
    }