update_viewer.cpp

ncbi源码

            if ((*d)->IsName())
                name = (*d)->GetName()[0];
            else if ((*d)->IsMolecule_type() &&
                     (*d)->GetMolecule_type() == CBiomol_descr::eMolecule_type_protein)
                isProtein = true;
            if (isProtein && name) break;
        }

        // get gi
        if ((*m)->IsSetSeq_id())
            sid = &((*m)->GetSeq_id());

        // add protein to list
        if (isProtein && name && sid != NULL) {
            moleculeIDs[sid] = (*m)->GetId().Get();
            chains.push_back(make_pair(sid, name));
        }
    }
    if (chains.size() == 0) {
        ERRORMSG("No protein chains found in this structure!");
        return;
    }

    // get protein name (PDB identifier)
    string pdbID;
    CBiostruc::TDescr::const_iterator n, ne = biostruc->GetDescr().end();
    for (n=biostruc->GetDescr().begin(); n!=ne; ++n) {
        if ((*n)->IsName()) {
            pdbID = (*n)->GetName();    // assume first 'name' is pdb id
            break;
        }
    }

    // which chains to align?
    vector < const CSeq_id * > sids;
    if (chains.size() == 1) {
        sids.push_back(chains[0].first);
    } else {
        wxString *choices = new wxString[chains.size()];
        int choice;
        for (choice=0; choice<chains.size(); ++choice)
            choices[choice].Printf("%s_%c %s",
                pdbID.c_str(), chains[choice].second, chains[choice].first->GetSeqIdString().c_str());
        wxArrayInt selections;
        int nsel = wxGetMultipleChoices(selections, "Which chain do you want to align?",
            "Select Chain", chains.size(), choices, *viewerWindow);
        if (nsel == 0) return;
        for (choice=0; choice<nsel; ++choice)
            sids.push_back(chains[selections[choice]].first);
    }

    SequenceList newSequences;
    SequenceSet::SequenceList::const_iterator s, se = master->parentSet->sequenceSet->sequences.end();
    map < const Sequence * , const CSeq_id * > seq2id;
    for (int j=0; j<sids.size(); ++j) {

        // first check to see if this sequence is already present
        for (s=master->parentSet->sequenceSet->sequences.begin(); s!=se; ++s) {
            if ((*s)->identifier->MatchesSeqId(*(sids[j]))) {
                TRACEMSG("using existing sequence for " << sids[j]->GetSeqIdString());
                newSequences.push_back(*s);
                seq2id[*s] = sids[j];
                break;
            }
        }

        if (s == se) {
            // if not, find the sequence in the list from the structure file
            BioseqRefList::iterator b, be = bioseqs.end();
            for (b=bioseqs.begin(); b!=be; ++b) {
                CBioseq::TId::const_iterator i, ie = (*b)->GetId().end();
                for (i=(*b)->GetId().begin(); i!=ie; ++i) {
                    if ((*i)->Match(*(sids[j])))
                        break;
                }
                if (i != ie) {
                    const Sequence *sequence = master->parentSet->CreateNewSequence(**b);
                    if (sequence) {
                        TRACEMSG("found Bioseq for " << sids[j]->GetSeqIdString());
                        newSequences.push_back(sequence);
                        seq2id[sequence] = sids[j];
                        break;
                    }
                }
            }
            if (b == be) {
                ERRORMSG("ImportStructure() - can't find " << sids[j]->GetSeqIdString() << " in bioseq list!");
//                return;
            }
        }
    }

    SequenceList::const_iterator w, we = newSequences.end();
    for (w=newSequences.begin(); w!=we; ++w) {

        // add MMDB id tag to Bioseq if not present already
        (*w)->AddMMDBAnnotTag(mmdbID);

        // add MMDB and molecule id to identifier if not already set
        if ((*w)->identifier->mmdbID == MoleculeIdentifier::VALUE_NOT_SET) {
            (const_cast<MoleculeIdentifier*>((*w)->identifier))->mmdbID = mmdbID;
        } else {
            if ((*w)->identifier->mmdbID != mmdbID)
                ERRORMSG("MMDB ID mismatch in sequence " << (*w)->identifier->ToString()
                    << "; " << (*w)->identifier->mmdbID << " vs " << mmdbID);
        }
        if (moleculeIDs.find(seq2id[*w]) != moleculeIDs.end()) {
            if ((*w)->identifier->moleculeID == MoleculeIdentifier::VALUE_NOT_SET) {
                (const_cast<MoleculeIdentifier*>((*w)->identifier))->moleculeID =
                    moleculeIDs[seq2id[*w]];
            } else {
                if ((*w)->identifier->moleculeID != moleculeIDs[seq2id[*w]])
                    ERRORMSG("Molecule ID mismatch in sequence " << (*w)->identifier->ToString());
            }
        } else
            ERRORMSG("No matching id for MMDB sequence " << (*w)->identifier->ToString());
    }

    // create null-alignment
    AlignmentList newAlignments;
    int masterFrom = -1, masterTo = -1;
    const BlockMultipleAlignment *multiple = alignmentManager->GetCurrentMultipleAlignment();
    if (multiple) {
        BlockMultipleAlignment::UngappedAlignedBlockList aBlocks;
        multiple->GetUngappedAlignedBlocks(&aBlocks);
        if (aBlocks.size() > 0) {
            masterFrom = aBlocks.front()->GetRangeOfRow(0)->from;
            masterTo = aBlocks.back()->GetRangeOfRow(0)->to;
        }
    }
    GetVASTAlignments(newSequences, master, &newAlignments,
        &pendingStructureAlignments, masterFrom, masterTo);

    // add new alignment to update list
    if (newAlignments.size() == newSequences.size())
        AddAlignments(newAlignments);
    else {
        ERRORMSG("UpdateViewer::ImportStructure() - no new alignments were created");
        DELETE_ALL_AND_CLEAR(newAlignments, AlignmentList);
        return;
    }

    // will add Biostruc and structure alignments to ASN data later on, after merge
    pendingStructures.push_back(biostruc);

    // inform user of success
    wxMessageBox("The structure has been successfully imported! However, it will not appear until you\n"
        "save this data to a file and then re-load it in a new session. And depending on the type\n"
        "of data, it still might not appear unless the corresponding new pairwise alignment has\n"
        "been merged into the multiple alignment.",
        "Structure Added", wxOK | wxICON_INFORMATION, *viewerWindow);
}

void UpdateViewer::SavePendingStructures(void)
{
    TRACEMSG("saving pending imported structures and structure alignments");
    StructureSet *sSet =
        (alignmentManager->GetCurrentMultipleAlignment() &&
         alignmentManager->GetCurrentMultipleAlignment()->GetMaster()) ?
         alignmentManager->GetCurrentMultipleAlignment()->GetMaster()->parentSet : NULL;
    if (!sSet) return;
    while (pendingStructures.size() > 0) {
        if (!sSet->AddBiostrucToASN(pendingStructures.front().GetPointer())) {
            ERRORMSG("UpdateViewer::SavePendingStructures() - error saving Biostruc");
            return;
        }
        pendingStructures.pop_front();
    }
    while (pendingStructureAlignments.size() > 0) {
        sSet->AddStructureAlignment(pendingStructureAlignments.front().structureAlignment.GetPointer(),
            pendingStructureAlignments.front().masterDomainID,
            pendingStructureAlignments.front().slaveDomainID);
        pendingStructureAlignments.pop_front();
    }
}

void UpdateViewer::BlastUpdate(BlockMultipleAlignment *alignment, bool usePSSMFromMultiple)
{
    const BlockMultipleAlignment *multipleForPSSM = alignmentManager->GetCurrentMultipleAlignment();
    if (usePSSMFromMultiple && !multipleForPSSM) {
        ERRORMSG("Can't do BLAST/PSSM when no multiple alignment is present");
        return;
    }

    // find alignment, and replace it with BLAST result
    AlignmentList::iterator a, ae = GetCurrentAlignments().end();
    for (a=GetCurrentAlignments().begin(); a!=ae; ++a) {
        if (*a != alignment) continue;

        // run BLAST between master and first slave (should be only one slave...)
        BLASTer::AlignmentList toRealign;
        toRealign.push_back(alignment);
        BLASTer::AlignmentList newAlignments;
        alignmentManager->blaster->CreateNewPairwiseAlignmentsByBlast(
            multipleForPSSM, toRealign, &newAlignments, usePSSMFromMultiple);
        if (newAlignments.size() != 1) {
            ERRORMSG("UpdateViewer::BlastUpdate() - CreateNewPairwiseAlignmentsByBlast() failed");
            return;
        }
        if (newAlignments.front()->NAlignedBlocks() == 0) {
            WARNINGMSG("alignment unchanged");
            delete newAlignments.front();
            return;
        }

        // replace alignment with BLAST result
        TRACEMSG("BLAST succeeded - replacing alignment");
        delete alignment;
        *a = newAlignments.front();
        break;
    }

    // recreate alignment display with new alignment
    AlignmentList copy = GetCurrentAlignments();
    GetCurrentAlignments().clear();
    GetCurrentDisplay()->Empty();
    AddAlignments(copy);
//    (*viewerWindow)->ScrollToColumn(GetCurrentDisplay()->GetStartingColumn());
}

static void MapSlaveToMaster(const BlockMultipleAlignment *alignment,
    int slaveRow, vector < int > *slave2master)
{
    slave2master->clear();
    slave2master->resize(alignment->GetSequenceOfRow(slaveRow)->Length(), -1);
    BlockMultipleAlignment::UngappedAlignedBlockList uaBlocks;
    alignment->GetUngappedAlignedBlocks(&uaBlocks);
    BlockMultipleAlignment::UngappedAlignedBlockList::const_iterator b, be = uaBlocks.end();
    for (b=uaBlocks.begin(); b!=be; ++b) {
        const Block::Range
            *masterRange = (*b)->GetRangeOfRow(0),
            *slaveRange = (*b)->GetRangeOfRow(slaveRow);
        for (int i=0; i<(*b)->width; ++i)
            (*slave2master)[slaveRange->from + i] = masterRange->from + i;
    }
}

static BlockMultipleAlignment * GetAlignmentByBestNeighbor(
    const BlockMultipleAlignment *multiple,
    const BLASTer::AlignmentList rowAlignments)
{
    if (rowAlignments.size() != multiple->NRows()) {
        ERRORMSG("GetAlignmentByBestNeighbor: wrong # alignments");
        return NULL;
    }

    // find best-scoring aligment above some threshold
    const BlockMultipleAlignment *bestMatchFromMultiple = NULL;
    int b, bestRow = -1;
    BLASTer::AlignmentList::const_iterator p, pe = rowAlignments.end();
    for (b=0, p=rowAlignments.begin(); p!=pe; ++b, ++p) {
        if (!bestMatchFromMultiple || (*p)->GetRowDouble(0) < bestMatchFromMultiple->GetRowDouble(0)) {
            bestMatchFromMultiple = *p;
            bestRow = b;
        }
    }
    if (!bestMatchFromMultiple || bestMatchFromMultiple->GetRowDouble(0) > 0.000001) {
        WARNINGMSG("GetAlignmentByBestNeighbor: no significant hit found");
        return NULL;
    }
    INFOMSG("Closest neighbor from multiple: sequence "
        << bestMatchFromMultiple->GetMaster()->identifier->ToString()
        << ", E-value: " << bestMatchFromMultiple->GetRowDouble(0));
    GlobalMessenger()->RemoveAllHighlights(true);
    GlobalMessenger()->AddHighlights(
        bestMatchFromMultiple->GetMaster(), 0, bestMatchFromMultiple->GetMaster()->Length()-1);

    // if the best match is the multiple's master, then just use that alignment
    if (bestRow == 0) return bestMatchFromMultiple->Clone();

    // otherwise, use best match as a guide alignment to align the slave with the multiple's master
    vector < int > import2slave, slave2master;
    MapSlaveToMaster(bestMatchFromMultiple, 1, &import2slave);
    MapSlaveToMaster(multiple, bestRow, &slave2master);

    const Sequence *importSeq = bestMatchFromMultiple->GetSequenceOfRow(1);
    BlockMultipleAlignment::SequenceList *seqs = new BlockMultipleAlignment::SequenceList(2);
    (*seqs)[0] = multiple->GetSequenceOfRow(0);
    (*seqs)[1] = importSeq;
    BlockMultipleAlignment *newAlignment =
        new BlockMultipleAlignment(seqs, importSeq->parentSet->alignmentManager);

    // create maximally sized blocks
    int masterStart, importStart, importLoc, slaveLoc, masterLoc, len;
    for (importStart=-1, importLoc=0; importLoc<=importSeq->Length(); ++importLoc) {

        // map import -> slave -> master
        slaveLoc = (importLoc<importSeq->Length()) ? import2slave[importLoc] : -1;
        masterLoc = (slaveLoc >= 0) ? slave2master[slaveLoc] : -1;

        // if we're currently inside a block..
        if (importStart >= 0) {

            // add another residue to a continuously aligned block
            if (masterLoc >= 0 && masterLoc-masterStart == importLoc-importStart) {
                ++len;
            }

            // terminate block
            else {
                UngappedAlignedBlock *newBlock = new UngappedAlignedBlock(newAlignment);
                newBlock->SetRangeOfRow(0, masterStart, masterStart + len - 1);
                newBlock->SetRangeOfRow(1, importStart, importStart + len - 1);
                newBlock->width = len;
                newAlignment->AddAlignedBlockAtEnd(newBlock);
                importStart = -1;
            }
        }

        // search for start of block
        if (importStart < 0) {
            if (masterLoc >= 0) {
                masterStart = masterLoc;
                importStart = importLoc;
                len = 1;
            }
        }
    }

    // finalize and and add new alignment to list
    if (!newAlignment->AddUnalignedBlocks() || !newAlignment->UpdateBlockMapAndColors(false)) {
        ERRORMSG("error finalizing alignment");
        delete newAlignment;
        return NULL;
    }

    return newAlignment;
}

void UpdateViewer::BlastNeighbor(BlockMultipleAlignment *update)
{
    const BlockMultipleAlignment *multiple = alignmentManager->GetCurrentMultipleAlignment();
    if (!multiple) {
        ERRORMSG("Can't do BLAST Neighbor when no multiple alignment is present");
        return;
    }
    BlockMultipleAlignment::UngappedAlignedBlockList uaBlocks;
    multiple->GetUngappedAlignedBlocks(&uaBlocks);
    if (uaBlocks.size() == 0) {
        ERRORMSG("Can't do BLAST Neighbor with null multiple alignment");
        return;
    }
    const Sequence *updateSeq = update->GetSequenceOfRow(1);