update_viewer.cpp

ncbi源码

}

void UpdateViewer::ReadSequencesFromFile(SequenceList *newSequences, StructureSet *sSet) const
{
    newSequences->clear();

    wxString fastaFile = wxFileSelector("Choose a FASTA file from which to import",
        "", "", "", "*.*", wxOPEN | wxFILE_MUST_EXIST, *viewerWindow);
    if (fastaFile.size() > 0) {
        FILE *fp = FileOpen(fastaFile.c_str(), "r");
        if (fp) {
            SeqEntry *sep = NULL;
            string err;
            while ((sep=FastaToSeqEntry(fp, FALSE)) != NULL) {

                // convert C to C++ SeqEntry
                CSeq_entry se;
                if (!ConvertAsnFromCToCPP(sep, (AsnWriteFunc) SeqEntryAsnWrite, &se, &err)) {
                    ERRORMSG("UpdateViewer::ImportSequence() - error converting to C++ object: "
                        << err);
                } else {
                    // create Sequence - just one from each Seq-entry for now
                    if (se.IsSeq()) {
                        const Sequence *sequence = sSet->CreateNewSequence(se.SetSeq());
                        if (sequence)
                            newSequences->push_back(sequence);
                    } else {
                        ERRORMSG("FastaToSeqEntry() returned Bioseq-set in Seq-entry");
                    }
                }
                SeqEntryFree(sep);
            }
            FileClose(fp);
        } else
            ERRORMSG("Unable to open " << fastaFile.c_str());
    }
}

static BlockMultipleAlignment * MakeEmptyAlignment(const Sequence *master, const Sequence *slave)
{
    BlockMultipleAlignment::SequenceList *seqs = new BlockMultipleAlignment::SequenceList(2);
    (*seqs)[0] = master;
    (*seqs)[1] = slave;
    BlockMultipleAlignment *newAlignment =
        new BlockMultipleAlignment(seqs, master->parentSet->alignmentManager);
    if (!newAlignment->AddUnalignedBlocks() || !newAlignment->UpdateBlockMapAndColors(false)) {
        ERRORMSG("MakeEmptyAlignment() - error finalizing alignment");
        delete newAlignment;
        return NULL;
    }
    return newAlignment;
}

void UpdateViewer::MakeEmptyAlignments(const SequenceList& newSequences,
    const Sequence *master, AlignmentList *newAlignments) const
{
    newAlignments->clear();
    SequenceList::const_iterator s, se = newSequences.end();
    for (s=newSequences.begin(); s!=se; ++s) {
        BlockMultipleAlignment *newAlignment = MakeEmptyAlignment(master, *s);
        if (newAlignment) newAlignments->push_back(newAlignment);
    }
}

void UpdateViewer::FetchSequences(StructureSet *sSet, SequenceList *newSequences) const
{
    // choose import type
    static const wxString choiceStrings[] = { "Network via GI/Accession", "From a FASTA File" };
    enum choiceValues { FROM_GI=0, FROM_FASTA, N_CHOICES };
    int importFrom = wxGetSingleChoiceIndex("From what source would you like to import sequences?",
        "Select Import Source", N_CHOICES, choiceStrings, *viewerWindow);
    if (importFrom < 0) return;     // cancelled

    // network import
    if (importFrom == FROM_GI)
        FetchSequencesViaHTTP(newSequences, sSet);

    // FASTA import
    else if (importFrom == FROM_FASTA)
        ReadSequencesFromFile(newSequences, sSet);
}

void UpdateViewer::ImportSequences(void)
{
    // determine the master sequence for new alignments
    const Sequence *master = GetMasterSequence();
    if (!master) return;

    // import the new sequence(s)
    SequenceList newSequences;
    FetchSequences(master->parentSet, &newSequences);
    if (newSequences.size() == 0) {
        WARNINGMSG("UpdateViewer::ImportSequence() - no sequences were imported");
        return;
    }

    // create null-alignments for each sequence
    AlignmentList newAlignments;
    MakeEmptyAlignments(newSequences, master, &newAlignments);

    // add new alignments to update list
    if (newAlignments.size() > 0)
        AddAlignments(newAlignments);
    else
        ERRORMSG("UpdateViewer::ImportSequence() - no new alignments were created");
}

void UpdateViewer::GetVASTAlignments(const SequenceList& newSequences,
    const Sequence *master, AlignmentList *newAlignments,
    PendingStructureAlignments *structureAlignments,
    int masterFrom, int masterTo) const
{
    if (master->identifier->pdbID.size() == 0) {
        WARNINGMSG("UpdateViewer::GetVASTAlignments() - "
            "can't be called with non-MMDB master " << master->identifier->ToString());
        return;
    }

    SequenceList::const_iterator s, se = newSequences.end();
    for (s=newSequences.begin(); s!=se; ++s) {
        if ((*s)->identifier->pdbID.size() == 0) {
            WARNINGMSG("UpdateViewer::GetVASTAlignments() - "
                "can't be called with non-MMDB slave " << (*s)->identifier->ToString());
            BlockMultipleAlignment *newAlignment = MakeEmptyAlignment(master, *s);
            if (newAlignment) newAlignments->push_back(newAlignment);
            continue;
        }

        // set up URL
        string
            host = "www.ncbi.nlm.nih.gov",
            path = "/Structure/VA/vastalign.cgi", err;
        CNcbiOstrstream argstr;
        argstr << "master=" << master->identifier->ToString()
            << "&slave=" << (*s)->identifier->ToString();
        if (masterFrom >= 0 && masterTo >= 0 && masterFrom <= masterTo &&
            masterFrom < master->Length() && masterTo < master->Length())
            argstr << "&from=" << (masterFrom+1) << "&to=" << (masterTo+1); // URL #'s are 1-based
        argstr << '\0';
        auto_ptr<char> args(argstr.str());

        // connect to vastalign.cgi
        CBiostruc_annot_set structureAlignment;
        INFOMSG("trying to load VAST alignment data from " << host << path << '?' << args.get());
        if (!GetAsnDataViaHTTP(host, path, args.get(), &structureAlignment, &err)) {
            ERRORMSG("Error calling vastalign.cgi: " << err);
            BlockMultipleAlignment *newAlignment = MakeEmptyAlignment(master, *s);
            if (newAlignment) newAlignments->push_back(newAlignment);
            continue;
        }
        INFOMSG("successfully loaded data from vastalign.cgi");
#ifdef _DEBUG
        WriteASNToFile("vastalign.dat.txt", structureAlignment, false, &err);
#endif

        // create initially empty alignment
        BlockMultipleAlignment::SequenceList *seqs = new BlockMultipleAlignment::SequenceList(2);
        (*seqs)[0] = master;
        (*seqs)[1] = *s;
        BlockMultipleAlignment *newAlignment =
            new BlockMultipleAlignment(seqs, master->parentSet->alignmentManager);

        // skip if no VAST alignment found
        if (structureAlignment.IsSetId() && structureAlignment.GetId().front()->IsMmdb_id() &&
            structureAlignment.GetId().front()->GetMmdb_id().Get() == 0)
        {
            WARNINGMSG("VAST found no alignment for these chains");
        }

        else try {

            if (!structureAlignment.IsSetId() || !structureAlignment.GetId().front()->IsMmdb_id() ||
                structureAlignment.GetFeatures().size() != 1 ||
                structureAlignment.GetFeatures().front()->GetFeatures().size() != 1 ||
                !structureAlignment.GetFeatures().front()->GetFeatures().front()->IsSetLocation() ||
                !structureAlignment.GetFeatures().front()->GetFeatures().front()->GetLocation().IsAlignment())
            {
                throw "VAST data does not contain exactly one alignment of recognized format - "
                    "possibly a problem with vastalign.cgi";
            }

            if (structureAlignment.GetId().front()->GetMmdb_id().Get() != master->identifier->mmdbID)
                throw "Master structure MMDB ID mismatch - check to see if this structure been updated";

            // load in alignment, check format
            const CChem_graph_alignment& alignment =
                structureAlignment.GetFeatures().front()->GetFeatures().front()->GetLocation().GetAlignment();
            if (alignment.GetDimension() != 2 || alignment.GetAlignment().size() != 2 ||
                alignment.GetBiostruc_ids().size() != 2 ||
                !alignment.GetBiostruc_ids().front()->IsMmdb_id() ||
                alignment.GetBiostruc_ids().front()->GetMmdb_id().Get() != master->identifier->mmdbID ||
                !alignment.GetBiostruc_ids().back()->IsMmdb_id() ||
                alignment.GetBiostruc_ids().back()->GetMmdb_id().Get() != (*s)->identifier->mmdbID ||
                !alignment.GetAlignment().front()->IsResidues() ||
                !alignment.GetAlignment().front()->GetResidues().IsInterval() ||
                !alignment.GetAlignment().back()->IsResidues() ||
                !alignment.GetAlignment().back()->GetResidues().IsInterval() ||
                alignment.GetAlignment().front()->GetResidues().GetInterval().size() !=
                    alignment.GetAlignment().back()->GetResidues().GetInterval().size())
            {
                throw "Unrecognized VAST data format";
            }

            // construct alignment from residue intervals
            CResidue_pntrs::TInterval::const_iterator i, j,
                ie = alignment.GetAlignment().front()->GetResidues().GetInterval().end();
            for (i=alignment.GetAlignment().front()->GetResidues().GetInterval().begin(),
                 j=alignment.GetAlignment().back()->GetResidues().GetInterval().begin(); i!=ie; ++i, ++j)
            {
                if ((*i)->GetMolecule_id().Get() != master->identifier->moleculeID ||
                    (*j)->GetMolecule_id().Get() != (*s)->identifier->moleculeID)
                {
                    throw "Mismatch in molecule ids in alignment interval block";
                }
                UngappedAlignedBlock *newBlock = new UngappedAlignedBlock(newAlignment);
                newBlock->SetRangeOfRow(0, (*i)->GetFrom().Get() - 1, (*i)->GetTo().Get() - 1);
                newBlock->SetRangeOfRow(1, (*j)->GetFrom().Get() - 1, (*j)->GetTo().Get() - 1);
                newBlock->width = (*i)->GetTo().Get() - (*i)->GetFrom().Get() + 1;
                newAlignment->AddAlignedBlockAtEnd(newBlock);
            }

            // add structure alignment to list
            if (alignment.GetTransform().size() == 1)
            {
                structureAlignments->resize(structureAlignments->size() + 1);
                structureAlignments->back().structureAlignment =
                    structureAlignment.SetFeatures().front()->SetFeatures().front();
                structureAlignments->back().masterDomainID =
                    structureAlignment.GetFeatures().front()->GetId().Get();
                structureAlignments->back().slaveDomainID =
                    structureAlignment.GetFeatures().front()->GetFeatures().front()->GetId().Get();
            } else
                throw "No structure alignment in VAST data blob";

        } catch (const char *err) {
            ERRORMSG("Failed to import VAST alignment: " << err);
        }

        // finalize alignment
        if (!newAlignment->AddUnalignedBlocks() || !newAlignment->UpdateBlockMapAndColors(false)) {
            ERRORMSG("MakeEmptyAlignment() - error finalizing alignment");
            delete newAlignment;
            continue;
        }
        newAlignments->push_back(newAlignment);
    }
}

void UpdateViewer::ImportStructure(void)
{
    // determine the master sequence for new alignments
    const Sequence *master = GetMasterSequence();
    if (!master) return;
    if (!master->molecule) {
        ERRORMSG("Can't import another structure unless master sequence has structure");
        return;
    }
    if (master->parentSet->objects.size() == 1 && master->parentSet->frameMap.size() != 1) {
        ERRORMSG("Can't import another structure when current structure has multiple models");
        return;
    }

    // choose import type for structure
    static const wxString choiceStrings[] = { "Via Network", "From a File" };
    enum choiceValues { FROM_NETWORK=0, FROM_FILE, N_CHOICES };
    int importFrom = wxGetSingleChoiceIndex(
        "From what source would you like to import the structure?", "Select Import Source",
        N_CHOICES, choiceStrings, *viewerWindow);
    if (importFrom < 0) return;     // cancelled

    CRef < CBiostruc > biostruc;
    BioseqRefList bioseqs;

    if (importFrom == FROM_NETWORK) {
        wxString id = wxGetTextFromUser("Enter a PDB or MMDB ID:", "Input Identifier", "", *viewerWindow);
        biostruc.Reset(new CBiostruc());
        if (!LoadStructureViaCache(id.c_str(),
                (master->parentSet->isAlphaOnly ? eModel_type_ncbi_backbone : eModel_type_ncbi_all_atom),
                biostruc, &bioseqs)) {
            ERRORMSG("Failed to load structure " << id.c_str());
            return;
        }
    }

    else if (importFrom == FROM_FILE) {
        string filename = wxFileSelector("Choose a single-structure file:",
            GetUserDir().c_str(), "", "", "*.*", wxOPEN | wxFILE_MUST_EXIST, *viewerWindow).c_str();
        if (filename.size() == 0) return;
        bool readOK = false;
        string err;
        TRACEMSG("trying to read file '" << filename << "' as binary mime");
        CRef < CNcbi_mime_asn1 > mime(new CNcbi_mime_asn1());
        SetDiagPostLevel(eDiag_Fatal); // ignore all but Fatal errors while reading data
        readOK = ReadASNFromFile(filename.c_str(), mime.GetPointer(), true, &err);
        SetDiagPostLevel(eDiag_Info);
        if (!readOK) {
            TRACEMSG("error: " << err);
            TRACEMSG("trying to read file '" << filename << "' as ascii mime");
            mime.Reset(new CNcbi_mime_asn1());
            SetDiagPostLevel(eDiag_Fatal); // ignore all but Fatal errors while reading data
            readOK = ReadASNFromFile(filename.c_str(), mime.GetPointer(), false, &err);
            SetDiagPostLevel(eDiag_Info);
        }
        if (!readOK) {
            TRACEMSG("error: " << err);
            ERRORMSG("Couldn't read structure from " << filename);
            return;
        }

        ExtractBiostrucAndBioseqs(*mime, biostruc, &bioseqs);
    }

    int mmdbID;
    if (biostruc->GetId().size() == 0 || !biostruc->GetId().front()->IsMmdb_id()) {
        ERRORMSG("Can't get MMDB ID from loaded structure");
        return;
    }
    mmdbID = biostruc->GetId().front()->GetMmdb_id().Get();

    // make sure Biostruc is of correct type
    if (biostruc->GetModel().size() != 1 ||
        (master->parentSet->isAlphaOnly &&
            biostruc->GetModel().front()->GetType() != eModel_type_ncbi_backbone) ||
        (!master->parentSet->isAlphaOnly &&
            biostruc->GetModel().front()->GetType() != eModel_type_ncbi_all_atom)) {
        ERRORMSG("Biostruc does not match current data - should be "
            << (master->parentSet->isAlphaOnly ? "alpha-only" : "one-coordinate-per-atom") << " model");
        return;
    }

    // make list of protein chains in this structure
    vector < pair < const CSeq_id * , char > > chains;  // holds Seq-id and chain name
    map < const CSeq_id * , int > moleculeIDs;          // maps Seq-id -> molecule ID within MMDB object
    CBiostruc_graph::TDescr::const_iterator d, de;
    CBiostruc_graph::TMolecule_graphs::const_iterator
        m, me = biostruc->GetChemical_graph().GetMolecule_graphs().end();
    for (m=biostruc->GetChemical_graph().GetMolecule_graphs().begin(); m!=me; ++m) {
        bool isProtein = false;
        const CSeq_id *sid = NULL;
        char name = 0;

        // check descr for chain name/type
        de = (*m)->GetDescr().end();
        for (d=(*m)->GetDescr().begin(); d!=de; ++d) {