开通博客赚积分
发布资源赚积分
分类

源码开发语言/平台
当前位置：首页 > 源码/资料 > 数值算法/人工智能 > 生物技术 > 查看源码
update_viewer.cpp
资源名称：ncbi_cxx_win--Jul_08_2004.zip　[点击查看]
上传用户：yhdzpy8989
上传日期：2007-06-13
资源大小：13604k
文件大小：53k
源码类别：

生物技术

开发平台：

C/C++

update_viewer.cpp：源码内容
  1. /*
  2.  * ===========================================================================
  3.  * PRODUCTION $Log: update_viewer.cpp,v $
  4.  * PRODUCTION Revision 1000.3  2004/06/01 18:29:49  gouriano
  5.  * PRODUCTION PRODUCTION: UPGRADED [GCC34_MSVC7] Dev-tree R1.73
  6.  * PRODUCTION
  7.  * ===========================================================================
  8.  */
  10. /*  $Id: update_viewer.cpp,v 1000.3 2004/06/01 18:29:49 gouriano Exp $
  11. * ===========================================================================
  12. *
  13. *                            PUBLIC DOMAIN NOTICE
  14. *               National Center for Biotechnology Information
  15. *
  16. *  This software/database is a "United States Government Work" under the
  17. *  terms of the United States Copyright Act.  It was written as part of
  18. *  the author's official duties as a United States Government employee and
  19. *  thus cannot be copyrighted.  This software/database is freely available
  20. *  to the public for use. The National Library of Medicine and the U.S.
  21. *  Government have not placed any restriction on its use or reproduction.
  22. *
  23. *  Although all reasonable efforts have been taken to ensure the accuracy
  24. *  and reliability of the software and data, the NLM and the U.S.
  25. *  Government do not and cannot warrant the performance or results that
  26. *  may be obtained by using this software or data. The NLM and the U.S.
  27. *  Government disclaim all warranties, express or implied, including
  28. *  warranties of performance, merchantability or fitness for any particular
  29. *  purpose.
  30. *
  31. *  Please cite the author in any work or product based on this material.
  32. *
  33. * ===========================================================================
  34. *
  35. * Authors:  Paul Thiessen
  36. *
  37. * File Description:
  38. *      implementation of non-GUI part of update viewer
  39. *
  40. * ===========================================================================
  41. */
  43. #ifdef _MSC_VER
  44. #pragma warning(disable:4018)   // disable signed/unsigned mismatch warning in MSVC
  45. #endif
  47. #include <ncbi_pch.hpp>
  48. #include <corelib/ncbistd.hpp>
  49. #include <corelib/ncbistl.hpp>
  51. #include <objects/cdd/Cdd.hpp>
  52. #include <objects/cdd/Update_align.hpp>
  53. #include <objects/cdd/Update_comment.hpp>
  54. #include <objects/seq/Seq_annot.hpp>
  55. #include <objects/seq/Bioseq.hpp>
  56. #include <objects/seqset/Seq_entry.hpp>
  57. #include <objects/seqset/Bioseq_set.hpp>
  58. #include <objects/ncbimime/Ncbi_mime_asn1.hpp>
  59. #include <objects/ncbimime/Biostruc_seq.hpp>
  60. #include <objects/mmdb2/Model_type.hpp>
  61. #include <objects/mmdb2/Biostruc_model.hpp>
  62. #include <objects/mmdb1/Biostruc_graph.hpp>
  63. #include <objects/mmdb1/Biostruc_descr.hpp>
  64. #include <objects/general/Dbtag.hpp>
  65. #include <objects/general/Object_id.hpp>
  66. #include <objects/mmdb1/Biostruc_id.hpp>
  67. #include <objects/mmdb1/Mmdb_id.hpp>
  68. #include <objects/mmdb1/Biostruc_annot_set.hpp>
  69. #include <objects/mmdb3/Biostruc_feature_set.hpp>
  70. #include <objects/mmdb3/Chem_graph_alignment.hpp>
  71. #include <objects/mmdb3/Chem_graph_pntrs.hpp>
  72. #include <objects/mmdb3/Residue_pntrs.hpp>
  73. #include <objects/mmdb3/Residue_interval_pntr.hpp>
  74. #include <objects/mmdb1/Molecule_id.hpp>
  75. #include <objects/mmdb1/Residue_id.hpp>
  76. #include <objects/mmdb3/Biostruc_feature_set_id.hpp>
  77. #include <objects/mmdb3/Biostruc_feature_id.hpp>
  79. #include <memory>
  80. #include <algorithm>
  82. #include "update_viewer.hpp"
  83. #include "asn_reader.hpp"
  84. #include "update_viewer_window.hpp"
  85. #include "messenger.hpp"
  86. #include "sequence_display.hpp"
  87. #include "cn3d_colors.hpp"
  88. #include "alignment_manager.hpp"
  89. #include "cn3d_threader.hpp"
  90. #include "structure_set.hpp"
  91. #include "molecule.hpp"
  92. #include "cn3d_tools.hpp"
  93. #include "asn_converter.hpp"
  94. #include "cn3d_blast.hpp"
  95. #include "molecule_identifier.hpp"
  96. #include "cn3d_cache.hpp"
  97. #include "cn3d_ba_interface.hpp"
  99. #include <wx/tokenzr.h>
  101. // C stuff
  102. #include <stdio.h>
  103. #include <tofasta.h>
  104. #include <objseq.h>
  105. #include <objsset.h>
  107. USING_NCBI_SCOPE;
  108. USING_SCOPE(objects);
  111. BEGIN_SCOPE(Cn3D)
  113. UpdateViewer::UpdateViewer(AlignmentManager *alnMgr) :
  114.     // not sure why this cast is necessary, but MSVC requires it...
  115.     ViewerBase(reinterpret_cast<ViewerWindowBase**>(&updateWindow), alnMgr),
  116.     updateWindow(NULL)
  117. {
  118.     // when first created, start with blank display
  119.     AlignmentList emptyAlignmentList;
  120.     SequenceDisplay *blankDisplay = new SequenceDisplay(true, viewerWindow);
  121.     InitData(&emptyAlignmentList, blankDisplay);
  122.     EnableStacks();
  123. }
  125. UpdateViewer::~UpdateViewer(void)
  126. {
  127. }
  129. void UpdateViewer::SetInitialState(void)
  130. {
  131.     KeepCurrent();
  132.     EnableStacks();
  133. }
  135. void UpdateViewer::SaveDialog(bool prompt)
  136. {
  137.     if (updateWindow) updateWindow->SaveDialog(prompt, false);
  138. }
  140. void UpdateViewer::CreateUpdateWindow(void)
  141. {
  142.     if (updateWindow) {
  143.         updateWindow->Show(true);
  144.         GlobalMessenger()->PostRedrawSequenceViewer(this);
  145.     } else {
  146.         SequenceDisplay *display = GetCurrentDisplay();
  147.         if (display) {
  148.             updateWindow = new UpdateViewerWindow(this);
  149.             updateWindow->NewDisplay(display, false);
  150.             updateWindow->ScrollToColumn(display->GetStartingColumn());
  151.             updateWindow->Show(true);
  152.             // ScrollTo causes immediate redraw, so don't need a second one
  153.             GlobalMessenger()->UnPostRedrawSequenceViewer(this);
  154.         }
  155.     }
  156. }
  158. void UpdateViewer::AddAlignments(const AlignmentList& newAlignments)
  159. {
  160.     AlignmentList& alignments = GetCurrentAlignments();
  161.     SequenceDisplay *display = GetCurrentDisplay();
  163.     // populate successive lines of the display with each alignment, with blank lines inbetween
  164.     AlignmentList::const_iterator a, ae = newAlignments.end();
  165.     int nViolations = 0;
  166.     for (a=newAlignments.begin(); a!=ae; ++a) {
  167.         if ((*a)->NRows() != 2) {
  168.             ERRORMSG("UpdateViewer::AddAlignments() - got alignment with "
  169.                 << (*a)->NRows() << " rows");
  170.             continue;
  171.         }
  173.         // add alignment to stack list
  174.         alignments.push_back(*a);
  176.         // add alignment to the display, including block row since editor is always on
  177.         if (display->NRows() != 0) display->AddRowFromString("");
  178.         display->AddBlockBoundaryRow(*a);
  179.         for (int row=0; row<2; ++row)
  180.             display->AddRowFromAlignment(row, *a);
  181.     }
  183.     if (alignments.size() > 0)
  184.         display->SetStartingColumn(alignments.front()->GetFirstAlignedBlockPosition() - 5);
  186.     Save();    // make this an undoable operation
  187. if (updateWindow) updateWindow->UpdateDisplay(display);
  188. }
  190. void UpdateViewer::ReplaceAlignments(const AlignmentList& alignmentList)
  191. {
  192.     // empty out the current alignment list and display (but not the undo stacks!)
  193.     DELETE_ALL_AND_CLEAR(GetCurrentAlignments(), AlignmentList);
  194.     GetCurrentDisplay()->Empty();
  196.     AddAlignments(alignmentList);
  197. }
  199. void UpdateViewer::DeleteAlignment(BlockMultipleAlignment *toDelete)
  200. {
  201.     AlignmentList keepAlignments;
  202.     AlignmentList::const_iterator a, ae = GetCurrentAlignments().end();
  203.     for (a=GetCurrentAlignments().begin(); a!=ae; ++a)
  204.         if (*a != toDelete)
  205.             keepAlignments.push_back((*a)->Clone());
  207.     ReplaceAlignments(keepAlignments);
  208. }
  210. void UpdateViewer::SaveAlignments(void)
  211. {
  212.     SetInitialState();
  214.     // construct a new list of ASN Update-aligns
  215.     CCdd::TPending updates;
  216.     map < CUpdate_align * , bool > usedUpdateAligns;
  218.     AlignmentList::const_iterator a, ae = GetCurrentAlignments().end();
  219.     for (a=GetCurrentAlignments().begin(); a!=ae; ++a) {
  221.         // create a Seq-align (with Dense-diags) out of this update
  222.         if ((*a)->NRows() != 2) {
  223.             ERRORMSG("CreateSeqAlignFromUpdate() - can only save pairwise updates");
  224.             continue;
  225.         }
  226.         BlockMultipleAlignment::UngappedAlignedBlockList blocks;
  227.         (*a)->GetUngappedAlignedBlocks(&blocks);
  228.         CSeq_align *newSeqAlign = CreatePairwiseSeqAlignFromMultipleRow(*a, blocks, 1);
  229.         if (!newSeqAlign) continue;
  231.         // the Update-align and Seq-annot's list of Seq-aligns where this new Seq-align will go
  232.         CUpdate_align *updateAlign = (*a)->updateOrigin.GetPointer();
  233.         CSeq_annot::C_Data::TAlign *seqAlignList = NULL;
  235.         // if this alignment came from an existing Update-align, then replace just the Seq-align
  236.         // so that the rest of the Update-align's comments/annotations are preserved
  237.         if (updateAlign) {
  238.             if (!updateAlign->IsSetSeqannot() || !updateAlign->GetSeqannot().GetData().IsAlign()) {
  239.                 ERRORMSG("UpdateViewer::SaveAlignments() - confused by Update-align format");
  240.                 continue;
  241.             }
  242.         }
  244.         // if this is a new update, create a new Update-align and tag it
  245.         else {
  246.             updateAlign = new CUpdate_align();
  248.             // set type field depending on whether demoted sequence has structure
  249.             updateAlign->SetType((*a)->GetSequenceOfRow(1)->molecule ?
  250.                 CUpdate_align::eType_demoted_3d : CUpdate_align::eType_demoted);
  252.             // add a text comment
  253.             CUpdate_comment *comment = new CUpdate_comment();
  254.             comment->SetComment("Created by demotion or import in Cn3D 4.0");
  255.             updateAlign->SetDescription().resize(updateAlign->GetDescription().size() + 1);
  256.             updateAlign->SetDescription().back().Reset(comment);
  258.             // create a new Seq-annot
  259.             CRef < CSeq_annot > seqAnnotRef(new CSeq_annot());
  260.             seqAnnotRef->SetData().SetAlign();
  261.             updateAlign->SetSeqannot(*seqAnnotRef);
  262.         }
  264.         // get Seq-align list
  265.         if (!updateAlign || !(seqAlignList = &(updateAlign->SetSeqannot().SetData().SetAlign()))) {
  266.             ERRORMSG("UpdateViewer::SaveAlignments() - can't find Update-align and/or Seq-align list");
  267.             continue;
  268.         }
  270.         // if this is the first re-use of this Update-align, then empty out all existing
  271.         // Seq-aligns and push it onto the new update list
  272.         if (usedUpdateAligns.find(updateAlign) == usedUpdateAligns.end()) {
  273.             seqAlignList->clear();
  274.             updates.resize(updates.size() + 1);
  275.             updates.back().Reset(updateAlign);
  276.             usedUpdateAligns[updateAlign] = true;
  277.         }
  279.         // finally, add the new Seq-align to the list
  280.         seqAlignList->resize(seqAlignList->size() + 1);
  281.         seqAlignList->back().Reset(newSeqAlign);
  282.     }
  284.     // save these changes to the ASN data
  285.     alignmentManager->ReplaceUpdatesInASN(updates);
  286. }
  288. const Sequence * UpdateViewer::GetMasterSequence(void) const
  289. {
  290.     const Sequence *master = NULL;
  292.     // if there's a multiple alignment in the sequence viewer, use same master as that
  293.     if (alignmentManager->GetCurrentMultipleAlignment()) {
  294.         master = alignmentManager->GetCurrentMultipleAlignment()->GetMaster();
  295.     }
  296.     // if there's already an update present, use same master as that
  297.     else if (GetCurrentAlignments().size() > 0) {
  298.         master = GetCurrentAlignments().front()->GetMaster();
  299.     }
  300.     // otherwise, this must be a single structure with no updates, so we need to pick one of its
  301.     // chains as the new master
  302.     else {
  303.         vector < const Sequence * > chains;
  304.         if (alignmentManager->GetStructureProteins(&chains)) {
  305.             if (chains.size() == 1) {
  306.                 master = chains[0];
  307.             } else {
  308.                 wxString *titles = new wxString[chains.size()];
  309.                 int choice;
  310.                 for (choice=0; choice<chains.size(); ++choice)
  311.                     titles[choice] = chains[choice]->identifier->ToString().c_str();
  312.                 choice = wxGetSingleChoiceIndex("Align to which protein chain?",
  313.                     "Select Chain", chains.size(), titles);
  314.                 if (choice >= 0)
  315.                     master = chains[choice];
  316.                 else    // cancelled
  317.                     return NULL;
  318.             }
  319.         }
  320.     }
  321.     if (!master)
  322.         ERRORMSG("UpdateViewer::GetMasterSequence() - no master sequence defined");
  323.     return master;
  324. }
  326. void UpdateViewer::FetchSequencesViaHTTP(SequenceList *newSequences, StructureSet *sSet) const
  327. {
  328.     wxString ids =
  329.         wxGetTextFromUser("Enter a list of protein GIs or Accessions:",
  330.             "Input Identifier", "", *viewerWindow);
  331.     if (ids.size() == 0) return;
  333.     wxStringTokenizer tkz(ids, " ,;trn", wxTOKEN_STRTOK);
  334.     while (tkz.HasMoreTokens()) {
  335.         wxString id = tkz.GetNextToken();
  336.         CRef < CBioseq > bioseq = FetchSequenceViaHTTP(id.c_str());
  337.         const Sequence *sequence = sSet->CreateNewSequence(*bioseq);
  338.         if (sequence) {
  339.             if (sequence->isProtein)
  340.                 newSequences->push_back(sequence);
  341.             else
  342.                 ERRORMSG("The sequence must be a protein");
  343.         }
  344.     }
  345. }
  347. void UpdateViewer::ReadSequencesFromFile(SequenceList *newSequences, StructureSet *sSet) const
  348. {
  349.     newSequences->clear();
  351.     wxString fastaFile = wxFileSelector("Choose a FASTA file from which to import",
  352.         "", "", "", "*.*", wxOPEN | wxFILE_MUST_EXIST, *viewerWindow);
  353.     if (fastaFile.size() > 0) {
  354.         FILE *fp = FileOpen(fastaFile.c_str(), "r");
  355.         if (fp) {
  356.             SeqEntry *sep = NULL;
  357.             string err;
  358.             while ((sep=FastaToSeqEntry(fp, FALSE)) != NULL) {
  360.                 // convert C to C++ SeqEntry
  361.                 CSeq_entry se;
  362.                 if (!ConvertAsnFromCToCPP(sep, (AsnWriteFunc) SeqEntryAsnWrite, &se, &err)) {
  363.                     ERRORMSG("UpdateViewer::ImportSequence() - error converting to C++ object: "
  364.                         << err);
  365.                 } else {
  366.                     // create Sequence - just one from each Seq-entry for now
  367.                     if (se.IsSeq()) {
  368.                         const Sequence *sequence = sSet->CreateNewSequence(se.SetSeq());
  369.                         if (sequence)
  370.                             newSequences->push_back(sequence);
  371.                     } else {
  372.                         ERRORMSG("FastaToSeqEntry() returned Bioseq-set in Seq-entry");
  373.                     }
  374.                 }
  375.                 SeqEntryFree(sep);
  376.             }
  377.             FileClose(fp);
  378.         } else
  379.             ERRORMSG("Unable to open " << fastaFile.c_str());
  380.     }
  381. }
  383. static BlockMultipleAlignment * MakeEmptyAlignment(const Sequence *master, const Sequence *slave)
  384. {
  385.     BlockMultipleAlignment::SequenceList *seqs = new BlockMultipleAlignment::SequenceList(2);
  386.     (*seqs)[0] = master;
  387.     (*seqs)[1] = slave;
  388.     BlockMultipleAlignment *newAlignment =
  389.         new BlockMultipleAlignment(seqs, master->parentSet->alignmentManager);
  390.     if (!newAlignment->AddUnalignedBlocks() || !newAlignment->UpdateBlockMapAndColors(false)) {
  391.         ERRORMSG("MakeEmptyAlignment() - error finalizing alignment");
  392.         delete newAlignment;
  393.         return NULL;
  394.     }
  395.     return newAlignment;
  396. }
  398. void UpdateViewer::MakeEmptyAlignments(const SequenceList& newSequences,
  399.     const Sequence *master, AlignmentList *newAlignments) const
  400. {
  401.     newAlignments->clear();
  402.     SequenceList::const_iterator s, se = newSequences.end();
  403.     for (s=newSequences.begin(); s!=se; ++s) {
  404.         BlockMultipleAlignment *newAlignment = MakeEmptyAlignment(master, *s);
  405.         if (newAlignment) newAlignments->push_back(newAlignment);
  406.     }
  407. }
  409. void UpdateViewer::FetchSequences(StructureSet *sSet, SequenceList *newSequences) const
  410. {
  411.     // choose import type
  412.     static const wxString choiceStrings[] = { "Network via GI/Accession", "From a FASTA File" };
  413.     enum choiceValues { FROM_GI=0, FROM_FASTA, N_CHOICES };
  414.     int importFrom = wxGetSingleChoiceIndex("From what source would you like to import sequences?",
  415.         "Select Import Source", N_CHOICES, choiceStrings, *viewerWindow);
  416.     if (importFrom < 0) return;     // cancelled
  418.     // network import
  419.     if (importFrom == FROM_GI)
  420.         FetchSequencesViaHTTP(newSequences, sSet);
  422.     // FASTA import
  423.     else if (importFrom == FROM_FASTA)
  424.         ReadSequencesFromFile(newSequences, sSet);
  425. }
  427. void UpdateViewer::ImportSequences(void)
  428. {
  429.     // determine the master sequence for new alignments
  430.     const Sequence *master = GetMasterSequence();
  431.     if (!master) return;
  433.     // import the new sequence(s)
  434.     SequenceList newSequences;
  435.     FetchSequences(master->parentSet, &newSequences);
  436.     if (newSequences.size() == 0) {
  437.         WARNINGMSG("UpdateViewer::ImportSequence() - no sequences were imported");
  438.         return;
  439.     }
  441.     // create null-alignments for each sequence
  442.     AlignmentList newAlignments;
  443.     MakeEmptyAlignments(newSequences, master, &newAlignments);
  445.     // add new alignments to update list
  446.     if (newAlignments.size() > 0)
  447.         AddAlignments(newAlignments);
  448.     else
  449.         ERRORMSG("UpdateViewer::ImportSequence() - no new alignments were created");
  450. }
  452. void UpdateViewer::GetVASTAlignments(const SequenceList& newSequences,
  453.     const Sequence *master, AlignmentList *newAlignments,
  454.     PendingStructureAlignments *structureAlignments,
  455.     int masterFrom, int masterTo) const
  456. {
  457.     if (master->identifier->pdbID.size() == 0) {
  458.         WARNINGMSG("UpdateViewer::GetVASTAlignments() - "
  459.             "can't be called with non-MMDB master " << master->identifier->ToString());
  460.         return;
  461.     }
  463.     SequenceList::const_iterator s, se = newSequences.end();
  464.     for (s=newSequences.begin(); s!=se; ++s) {
  465.         if ((*s)->identifier->pdbID.size() == 0) {
  466.             WARNINGMSG("UpdateViewer::GetVASTAlignments() - "
  467.                 "can't be called with non-MMDB slave " << (*s)->identifier->ToString());
  468.             BlockMultipleAlignment *newAlignment = MakeEmptyAlignment(master, *s);
  469.             if (newAlignment) newAlignments->push_back(newAlignment);
  470.             continue;
  471.         }
  473.         // set up URL
  474.         string
  475.             host = "www.ncbi.nlm.nih.gov",
  476.             path = "/Structure/VA/vastalign.cgi", err;
  477.         CNcbiOstrstream argstr;
  478.         argstr << "master=" << master->identifier->ToString()
  479.             << "&slave=" << (*s)->identifier->ToString();
  480.         if (masterFrom >= 0 && masterTo >= 0 && masterFrom <= masterTo &&
  481.             masterFrom < master->Length() && masterTo < master->Length())
  482.             argstr << "&from=" << (masterFrom+1) << "&to=" << (masterTo+1); // URL #'s are 1-based
  483.         argstr << '';
  484.         auto_ptr<char> args(argstr.str());
  486.         // connect to vastalign.cgi
  487.         CBiostruc_annot_set structureAlignment;
  488.         INFOMSG("trying to load VAST alignment data from " << host << path << '?' << args.get());
  489.         if (!GetAsnDataViaHTTP(host, path, args.get(), &structureAlignment, &err)) {
  490.             ERRORMSG("Error calling vastalign.cgi: " << err);
  491.             BlockMultipleAlignment *newAlignment = MakeEmptyAlignment(master, *s);
  492.             if (newAlignment) newAlignments->push_back(newAlignment);
  493.             continue;
  494.         }
  495.         INFOMSG("successfully loaded data from vastalign.cgi");
  496. #ifdef _DEBUG
  497.         WriteASNToFile("vastalign.dat.txt", structureAlignment, false, &err);
  498. #endif
  500.         // create initially empty alignment
  501.         BlockMultipleAlignment::SequenceList *seqs = new BlockMultipleAlignment::SequenceList(2);
  502.         (*seqs)[0] = master;
  503.         (*seqs)[1] = *s;
  504.         BlockMultipleAlignment *newAlignment =
  505.             new BlockMultipleAlignment(seqs, master->parentSet->alignmentManager);
  507.         // skip if no VAST alignment found
  508.         if (structureAlignment.IsSetId() && structureAlignment.GetId().front()->IsMmdb_id() &&
  509.             structureAlignment.GetId().front()->GetMmdb_id().Get() == 0)
  510.         {
  511.             WARNINGMSG("VAST found no alignment for these chains");
  512.         }
  514.         else try {
  516.             if (!structureAlignment.IsSetId() || !structureAlignment.GetId().front()->IsMmdb_id() ||
  517.                 structureAlignment.GetFeatures().size() != 1 ||
  518.                 structureAlignment.GetFeatures().front()->GetFeatures().size() != 1 ||
  519.                 !structureAlignment.GetFeatures().front()->GetFeatures().front()->IsSetLocation() ||
  520.                 !structureAlignment.GetFeatures().front()->GetFeatures().front()->GetLocation().IsAlignment())
  521.             {
  522.                 throw "VAST data does not contain exactly one alignment of recognized format - "
  523.                     "possibly a problem with vastalign.cgi";
  524.             }
  526.             if (structureAlignment.GetId().front()->GetMmdb_id().Get() != master->identifier->mmdbID)
  527.                 throw "Master structure MMDB ID mismatch - check to see if this structure been updated";
  529.             // load in alignment, check format
  530.             const CChem_graph_alignment& alignment =
  531.                 structureAlignment.GetFeatures().front()->GetFeatures().front()->GetLocation().GetAlignment();
  532.             if (alignment.GetDimension() != 2 || alignment.GetAlignment().size() != 2 ||
  533.                 alignment.GetBiostruc_ids().size() != 2 ||
  534.                 !alignment.GetBiostruc_ids().front()->IsMmdb_id() ||
  535.                 alignment.GetBiostruc_ids().front()->GetMmdb_id().Get() != master->identifier->mmdbID ||
  536.                 !alignment.GetBiostruc_ids().back()->IsMmdb_id() ||
  537.                 alignment.GetBiostruc_ids().back()->GetMmdb_id().Get() != (*s)->identifier->mmdbID ||
  538.                 !alignment.GetAlignment().front()->IsResidues() ||
  539.                 !alignment.GetAlignment().front()->GetResidues().IsInterval() ||
  540.                 !alignment.GetAlignment().back()->IsResidues() ||
  541.                 !alignment.GetAlignment().back()->GetResidues().IsInterval() ||
  542.                 alignment.GetAlignment().front()->GetResidues().GetInterval().size() !=
  543.                     alignment.GetAlignment().back()->GetResidues().GetInterval().size())
  544.             {
  545.                 throw "Unrecognized VAST data format";
  546.             }
  548.             // construct alignment from residue intervals
  549.             CResidue_pntrs::TInterval::const_iterator i, j,
  550.                 ie = alignment.GetAlignment().front()->GetResidues().GetInterval().end();
  551.             for (i=alignment.GetAlignment().front()->GetResidues().GetInterval().begin(),
  552.                  j=alignment.GetAlignment().back()->GetResidues().GetInterval().begin(); i!=ie; ++i, ++j)
  553.             {
  554.                 if ((*i)->GetMolecule_id().Get() != master->identifier->moleculeID ||
  555.                     (*j)->GetMolecule_id().Get() != (*s)->identifier->moleculeID)
  556.                 {
  557.                     throw "Mismatch in molecule ids in alignment interval block";
  558.                 }
  559.                 UngappedAlignedBlock *newBlock = new UngappedAlignedBlock(newAlignment);
  560.                 newBlock->SetRangeOfRow(0, (*i)->GetFrom().Get() - 1, (*i)->GetTo().Get() - 1);
  561.                 newBlock->SetRangeOfRow(1, (*j)->GetFrom().Get() - 1, (*j)->GetTo().Get() - 1);
  562.                 newBlock->width = (*i)->GetTo().Get() - (*i)->GetFrom().Get() + 1;
  563.                 newAlignment->AddAlignedBlockAtEnd(newBlock);
  564.             }
  566.             // add structure alignment to list
  567.             if (alignment.GetTransform().size() == 1)
  568.             {
  569.                 structureAlignments->resize(structureAlignments->size() + 1);
  570.                 structureAlignments->back().structureAlignment =
  571.                     structureAlignment.SetFeatures().front()->SetFeatures().front();
  572.                 structureAlignments->back().masterDomainID =
  573.                     structureAlignment.GetFeatures().front()->GetId().Get();
  574.                 structureAlignments->back().slaveDomainID =
  575.                     structureAlignment.GetFeatures().front()->GetFeatures().front()->GetId().Get();
  576.             } else
  577.                 throw "No structure alignment in VAST data blob";
  579.         } catch (const char *err) {
  580.             ERRORMSG("Failed to import VAST alignment: " << err);
  581.         }
  583.         // finalize alignment
  584.         if (!newAlignment->AddUnalignedBlocks() || !newAlignment->UpdateBlockMapAndColors(false)) {
  585.             ERRORMSG("MakeEmptyAlignment() - error finalizing alignment");
  586.             delete newAlignment;
  587.             continue;
  588.         }
  589.         newAlignments->push_back(newAlignment);
  590.     }
  591. }
  593. void UpdateViewer::ImportStructure(void)
  594. {
  595.     // determine the master sequence for new alignments
  596.     const Sequence *master = GetMasterSequence();
  597.     if (!master) return;
  598.     if (!master->molecule) {
  599.         ERRORMSG("Can't import another structure unless master sequence has structure");
  600.         return;
  601.     }
  602.     if (master->parentSet->objects.size() == 1 && master->parentSet->frameMap.size() != 1) {
  603.         ERRORMSG("Can't import another structure when current structure has multiple models");
  604.         return;
  605.     }
  607.     // choose import type for structure
  608.     static const wxString choiceStrings[] = { "Via Network", "From a File" };
  609.     enum choiceValues { FROM_NETWORK=0, FROM_FILE, N_CHOICES };
  610.     int importFrom = wxGetSingleChoiceIndex(
  611.         "From what source would you like to import the structure?", "Select Import Source",
  612.         N_CHOICES, choiceStrings, *viewerWindow);
  613.     if (importFrom < 0) return;     // cancelled
  615.     CRef < CBiostruc > biostruc;
  616.     BioseqRefList bioseqs;
  618.     if (importFrom == FROM_NETWORK) {
  619.         wxString id = wxGetTextFromUser("Enter a PDB or MMDB ID:", "Input Identifier", "", *viewerWindow);
  620.         biostruc.Reset(new CBiostruc());
  621.         if (!LoadStructureViaCache(id.c_str(),
  622.                 (master->parentSet->isAlphaOnly ? eModel_type_ncbi_backbone : eModel_type_ncbi_all_atom),
  623.                 biostruc, &bioseqs)) {
  624.             ERRORMSG("Failed to load structure " << id.c_str());
  625.             return;
  626.         }
  627.     }
  629.     else if (importFrom == FROM_FILE) {
  630.         string filename = wxFileSelector("Choose a single-structure file:",
  631.             GetUserDir().c_str(), "", "", "*.*", wxOPEN | wxFILE_MUST_EXIST, *viewerWindow).c_str();
  632.         if (filename.size() == 0) return;
  633.         bool readOK = false;
  634.         string err;
  635.         TRACEMSG("trying to read file '" << filename << "' as binary mime");
  636.         CRef < CNcbi_mime_asn1 > mime(new CNcbi_mime_asn1());
  637.         SetDiagPostLevel(eDiag_Fatal); // ignore all but Fatal errors while reading data
  638.         readOK = ReadASNFromFile(filename.c_str(), mime.GetPointer(), true, &err);
  639.         SetDiagPostLevel(eDiag_Info);
  640.         if (!readOK) {
  641.             TRACEMSG("error: " << err);
  642.             TRACEMSG("trying to read file '" << filename << "' as ascii mime");
  643.             mime.Reset(new CNcbi_mime_asn1());
  644.             SetDiagPostLevel(eDiag_Fatal); // ignore all but Fatal errors while reading data
  645.             readOK = ReadASNFromFile(filename.c_str(), mime.GetPointer(), false, &err);
  646.             SetDiagPostLevel(eDiag_Info);
  647.         }
  648.         if (!readOK) {
  649.             TRACEMSG("error: " << err);
  650.             ERRORMSG("Couldn't read structure from " << filename);
  651.             return;
  652.         }
  654.         ExtractBiostrucAndBioseqs(*mime, biostruc, &bioseqs);
  655.     }
  657.     int mmdbID;
  658.     if (biostruc->GetId().size() == 0 || !biostruc->GetId().front()->IsMmdb_id()) {
  659.         ERRORMSG("Can't get MMDB ID from loaded structure");
  660.         return;
  661.     }
  662.     mmdbID = biostruc->GetId().front()->GetMmdb_id().Get();
  664.     // make sure Biostruc is of correct type
  665.     if (biostruc->GetModel().size() != 1 ||
  666.         (master->parentSet->isAlphaOnly &&
  667.             biostruc->GetModel().front()->GetType() != eModel_type_ncbi_backbone) ||
  668.         (!master->parentSet->isAlphaOnly &&
  669.             biostruc->GetModel().front()->GetType() != eModel_type_ncbi_all_atom)) {
  670.         ERRORMSG("Biostruc does not match current data - should be "
  671.             << (master->parentSet->isAlphaOnly ? "alpha-only" : "one-coordinate-per-atom") << " model");
  672.         return;
  673.     }
  675.     // make list of protein chains in this structure
  676.     vector < pair < const CSeq_id * , char > > chains;  // holds Seq-id and chain name
  677.     map < const CSeq_id * , int > moleculeIDs;          // maps Seq-id -> molecule ID within MMDB object
  678.     CBiostruc_graph::TDescr::const_iterator d, de;
  679.     CBiostruc_graph::TMolecule_graphs::const_iterator
  680.         m, me = biostruc->GetChemical_graph().GetMolecule_graphs().end();
  681.     for (m=biostruc->GetChemical_graph().GetMolecule_graphs().begin(); m!=me; ++m) {
  682.         bool isProtein = false;
  683.         const CSeq_id *sid = NULL;
  684.         char name = 0;
  686.         // check descr for chain name/type
  687.         de = (*m)->GetDescr().end();
  688.         for (d=(*m)->GetDescr().begin(); d!=de; ++d) {
  689.             if ((*d)->IsName())
  690.                 name = (*d)->GetName()[0];
  691.             else if ((*d)->IsMolecule_type() &&
  692.                      (*d)->GetMolecule_type() == CBiomol_descr::eMolecule_type_protein)
  693.                 isProtein = true;
  694.             if (isProtein && name) break;
  695.         }
  697.         // get gi
  698.         if ((*m)->IsSetSeq_id())
  699.             sid = &((*m)->GetSeq_id());
  701.         // add protein to list
  702.         if (isProtein && name && sid != NULL) {
  703.             moleculeIDs[sid] = (*m)->GetId().Get();
  704.             chains.push_back(make_pair(sid, name));
  705.         }
  706.     }
  707.     if (chains.size() == 0) {
  708.         ERRORMSG("No protein chains found in this structure!");
  709.         return;
  710.     }
  712.     // get protein name (PDB identifier)
  713.     string pdbID;
  714.     CBiostruc::TDescr::const_iterator n, ne = biostruc->GetDescr().end();
  715.     for (n=biostruc->GetDescr().begin(); n!=ne; ++n) {
  716.         if ((*n)->IsName()) {
  717.             pdbID = (*n)->GetName();    // assume first 'name' is pdb id
  718.             break;
  719.         }
  720.     }
  722.     // which chains to align?
  723.     vector < const CSeq_id * > sids;
  724.     if (chains.size() == 1) {
  725.         sids.push_back(chains[0].first);
  726.     } else {
  727.         wxString *choices = new wxString[chains.size()];
  728.         int choice;
  729.         for (choice=0; choice<chains.size(); ++choice)
  730.             choices[choice].Printf("%s_%c %s",
  731.                 pdbID.c_str(), chains[choice].second, chains[choice].first->GetSeqIdString().c_str());
  732.         wxArrayInt selections;
  733.         int nsel = wxGetMultipleChoices(selections, "Which chain do you want to align?",
  734.             "Select Chain", chains.size(), choices, *viewerWindow);
  735.         if (nsel == 0) return;
  736.         for (choice=0; choice<nsel; ++choice)
  737.             sids.push_back(chains[selections[choice]].first);
  738.     }
  740.     SequenceList newSequences;
  741.     SequenceSet::SequenceList::const_iterator s, se = master->parentSet->sequenceSet->sequences.end();
  742.     map < const Sequence * , const CSeq_id * > seq2id;
  743.     for (int j=0; j<sids.size(); ++j) {
  745.         // first check to see if this sequence is already present
  746.         for (s=master->parentSet->sequenceSet->sequences.begin(); s!=se; ++s) {
  747.             if ((*s)->identifier->MatchesSeqId(*(sids[j]))) {
  748.                 TRACEMSG("using existing sequence for " << sids[j]->GetSeqIdString());
  749.                 newSequences.push_back(*s);
  750.                 seq2id[*s] = sids[j];
  751.                 break;
  752.             }
  753.         }
  755.         if (s == se) {
  756.             // if not, find the sequence in the list from the structure file
  757.             BioseqRefList::iterator b, be = bioseqs.end();
  758.             for (b=bioseqs.begin(); b!=be; ++b) {
  759.                 CBioseq::TId::const_iterator i, ie = (*b)->GetId().end();
  760.                 for (i=(*b)->GetId().begin(); i!=ie; ++i) {
  761.                     if ((*i)->Match(*(sids[j])))
  762.                         break;
  763.                 }
  764.                 if (i != ie) {
  765.                     const Sequence *sequence = master->parentSet->CreateNewSequence(**b);
  766.                     if (sequence) {
  767.                         TRACEMSG("found Bioseq for " << sids[j]->GetSeqIdString());
  768.                         newSequences.push_back(sequence);
  769.                         seq2id[sequence] = sids[j];
  770.                         break;
  771.                     }
  772.                 }
  773.             }
  774.             if (b == be) {
  775.                 ERRORMSG("ImportStructure() - can't find " << sids[j]->GetSeqIdString() << " in bioseq list!");
  776. //                return;
  777.             }
  778.         }
  779.     }
  781.     SequenceList::const_iterator w, we = newSequences.end();
  782.     for (w=newSequences.begin(); w!=we; ++w) {
  784.         // add MMDB id tag to Bioseq if not present already
  785.         (*w)->AddMMDBAnnotTag(mmdbID);
  787.         // add MMDB and molecule id to identifier if not already set
  788.         if ((*w)->identifier->mmdbID == MoleculeIdentifier::VALUE_NOT_SET) {
  789.             (const_cast<MoleculeIdentifier*>((*w)->identifier))->mmdbID = mmdbID;
  790.         } else {
  791.             if ((*w)->identifier->mmdbID != mmdbID)
  792.                 ERRORMSG("MMDB ID mismatch in sequence " << (*w)->identifier->ToString()
  793.                     << "; " << (*w)->identifier->mmdbID << " vs " << mmdbID);
  794.         }
  795.         if (moleculeIDs.find(seq2id[*w]) != moleculeIDs.end()) {
  796.             if ((*w)->identifier->moleculeID == MoleculeIdentifier::VALUE_NOT_SET) {
  797.                 (const_cast<MoleculeIdentifier*>((*w)->identifier))->moleculeID =
  798.                     moleculeIDs[seq2id[*w]];
  799.             } else {
  800.                 if ((*w)->identifier->moleculeID != moleculeIDs[seq2id[*w]])
  801.                     ERRORMSG("Molecule ID mismatch in sequence " << (*w)->identifier->ToString());
  802.             }
  803.         } else
  804.             ERRORMSG("No matching id for MMDB sequence " << (*w)->identifier->ToString());
  805.     }
  807.     // create null-alignment
  808.     AlignmentList newAlignments;
  809.     int masterFrom = -1, masterTo = -1;
  810.     const BlockMultipleAlignment *multiple = alignmentManager->GetCurrentMultipleAlignment();
  811.     if (multiple) {
  812.         BlockMultipleAlignment::UngappedAlignedBlockList aBlocks;
  813.         multiple->GetUngappedAlignedBlocks(&aBlocks);
  814.         if (aBlocks.size() > 0) {
  815.             masterFrom = aBlocks.front()->GetRangeOfRow(0)->from;
  816.             masterTo = aBlocks.back()->GetRangeOfRow(0)->to;
  817.         }
  818.     }
  819.     GetVASTAlignments(newSequences, master, &newAlignments,
  820.         &pendingStructureAlignments, masterFrom, masterTo);
  822.     // add new alignment to update list
  823.     if (newAlignments.size() == newSequences.size())
  824.         AddAlignments(newAlignments);
  825.     else {
  826.         ERRORMSG("UpdateViewer::ImportStructure() - no new alignments were created");
  827.         DELETE_ALL_AND_CLEAR(newAlignments, AlignmentList);
  828.         return;
  829.     }
  831.     // will add Biostruc and structure alignments to ASN data later on, after merge
  832.     pendingStructures.push_back(biostruc);
  834.     // inform user of success
  835.     wxMessageBox("The structure has been successfully imported! However, it will not appear until youn"
  836.         "save this data to a file and then re-load it in a new session. And depending on the typen"
  837.         "of data, it still might not appear unless the corresponding new pairwise alignment hasn"
  838.         "been merged into the multiple alignment.",
  839.         "Structure Added", wxOK | wxICON_INFORMATION, *viewerWindow);
  840. }
  842. void UpdateViewer::SavePendingStructures(void)
  843. {
  844.     TRACEMSG("saving pending imported structures and structure alignments");
  845.     StructureSet *sSet =
  846.         (alignmentManager->GetCurrentMultipleAlignment() &&
  847.          alignmentManager->GetCurrentMultipleAlignment()->GetMaster()) ?
  848.          alignmentManager->GetCurrentMultipleAlignment()->GetMaster()->parentSet : NULL;
  849.     if (!sSet) return;
  850.     while (pendingStructures.size() > 0) {
  851.         if (!sSet->AddBiostrucToASN(pendingStructures.front().GetPointer())) {
  852.             ERRORMSG("UpdateViewer::SavePendingStructures() - error saving Biostruc");
  853.             return;
  854.         }
  855.         pendingStructures.pop_front();
  856.     }
  857.     while (pendingStructureAlignments.size() > 0) {
  858.         sSet->AddStructureAlignment(pendingStructureAlignments.front().structureAlignment.GetPointer(),
  859.             pendingStructureAlignments.front().masterDomainID,
  860.             pendingStructureAlignments.front().slaveDomainID);
  861.         pendingStructureAlignments.pop_front();
  862.     }
  863. }
  865. void UpdateViewer::BlastUpdate(BlockMultipleAlignment *alignment, bool usePSSMFromMultiple)
  866. {
  867.     const BlockMultipleAlignment *multipleForPSSM = alignmentManager->GetCurrentMultipleAlignment();
  868.     if (usePSSMFromMultiple && !multipleForPSSM) {
  869.         ERRORMSG("Can't do BLAST/PSSM when no multiple alignment is present");
  870.         return;
  871.     }
  873.     // find alignment, and replace it with BLAST result
  874.     AlignmentList::iterator a, ae = GetCurrentAlignments().end();
  875.     for (a=GetCurrentAlignments().begin(); a!=ae; ++a) {
  876.         if (*a != alignment) continue;
  878.         // run BLAST between master and first slave (should be only one slave...)
  879.         BLASTer::AlignmentList toRealign;
  880.         toRealign.push_back(alignment);
  881.         BLASTer::AlignmentList newAlignments;
  882.         alignmentManager->blaster->CreateNewPairwiseAlignmentsByBlast(
  883.             multipleForPSSM, toRealign, &newAlignments, usePSSMFromMultiple);
  884.         if (newAlignments.size() != 1) {
  885.             ERRORMSG("UpdateViewer::BlastUpdate() - CreateNewPairwiseAlignmentsByBlast() failed");
  886.             return;
  887.         }
  888.         if (newAlignments.front()->NAlignedBlocks() == 0) {
  889.             WARNINGMSG("alignment unchanged");
  890.             delete newAlignments.front();
  891.             return;
  892.         }
  894.         // replace alignment with BLAST result
  895.         TRACEMSG("BLAST succeeded - replacing alignment");
  896.         delete alignment;
  897.         *a = newAlignments.front();
  898.         break;
  899.     }
  901.     // recreate alignment display with new alignment
  902.     AlignmentList copy = GetCurrentAlignments();
  903.     GetCurrentAlignments().clear();
  904.     GetCurrentDisplay()->Empty();
  905.     AddAlignments(copy);
  906. //    (*viewerWindow)->ScrollToColumn(GetCurrentDisplay()->GetStartingColumn());
  907. }
  909. static void MapSlaveToMaster(const BlockMultipleAlignment *alignment,
  910.     int slaveRow, vector < int > *slave2master)
  911. {
  912.     slave2master->clear();
  913.     slave2master->resize(alignment->GetSequenceOfRow(slaveRow)->Length(), -1);
  914.     BlockMultipleAlignment::UngappedAlignedBlockList uaBlocks;
  915.     alignment->GetUngappedAlignedBlocks(&uaBlocks);
  916.     BlockMultipleAlignment::UngappedAlignedBlockList::const_iterator b, be = uaBlocks.end();
  917.     for (b=uaBlocks.begin(); b!=be; ++b) {
  918.         const Block::Range
  919.             *masterRange = (*b)->GetRangeOfRow(0),
  920.             *slaveRange = (*b)->GetRangeOfRow(slaveRow);
  921.         for (int i=0; i<(*b)->width; ++i)
  922.             (*slave2master)[slaveRange->from + i] = masterRange->from + i;
  923.     }
  924. }
  926. static BlockMultipleAlignment * GetAlignmentByBestNeighbor(
  927.     const BlockMultipleAlignment *multiple,
  928.     const BLASTer::AlignmentList rowAlignments)
  929. {
  930.     if (rowAlignments.size() != multiple->NRows()) {
  931.         ERRORMSG("GetAlignmentByBestNeighbor: wrong # alignments");
  932.         return NULL;
  933.     }
  935.     // find best-scoring aligment above some threshold
  936.     const BlockMultipleAlignment *bestMatchFromMultiple = NULL;
  937.     int b, bestRow = -1;
  938.     BLASTer::AlignmentList::const_iterator p, pe = rowAlignments.end();
  939.     for (b=0, p=rowAlignments.begin(); p!=pe; ++b, ++p) {
  940.         if (!bestMatchFromMultiple || (*p)->GetRowDouble(0) < bestMatchFromMultiple->GetRowDouble(0)) {
  941.             bestMatchFromMultiple = *p;
  942.             bestRow = b;
  943.         }
  944.     }
  945.     if (!bestMatchFromMultiple || bestMatchFromMultiple->GetRowDouble(0) > 0.000001) {
  946.         WARNINGMSG("GetAlignmentByBestNeighbor: no significant hit found");
  947.         return NULL;
  948.     }
  949.     INFOMSG("Closest neighbor from multiple: sequence "
  950.         << bestMatchFromMultiple->GetMaster()->identifier->ToString()
  951.         << ", E-value: " << bestMatchFromMultiple->GetRowDouble(0));
  952.     GlobalMessenger()->RemoveAllHighlights(true);
  953.     GlobalMessenger()->AddHighlights(
  954.         bestMatchFromMultiple->GetMaster(), 0, bestMatchFromMultiple->GetMaster()->Length()-1);
  956.     // if the best match is the multiple's master, then just use that alignment
  957.     if (bestRow == 0) return bestMatchFromMultiple->Clone();
  959.     // otherwise, use best match as a guide alignment to align the slave with the multiple's master
  960.     vector < int > import2slave, slave2master;
  961.     MapSlaveToMaster(bestMatchFromMultiple, 1, &import2slave);
  962.     MapSlaveToMaster(multiple, bestRow, &slave2master);
  964.     const Sequence *importSeq = bestMatchFromMultiple->GetSequenceOfRow(1);
  965.     BlockMultipleAlignment::SequenceList *seqs = new BlockMultipleAlignment::SequenceList(2);
  966.     (*seqs)[0] = multiple->GetSequenceOfRow(0);
  967.     (*seqs)[1] = importSeq;
  968.     BlockMultipleAlignment *newAlignment =
  969.         new BlockMultipleAlignment(seqs, importSeq->parentSet->alignmentManager);
  971.     // create maximally sized blocks
  972.     int masterStart, importStart, importLoc, slaveLoc, masterLoc, len;
  973.     for (importStart=-1, importLoc=0; importLoc<=importSeq->Length(); ++importLoc) {
  975.         // map import -> slave -> master
  976.         slaveLoc = (importLoc<importSeq->Length()) ? import2slave[importLoc] : -1;
  977.         masterLoc = (slaveLoc >= 0) ? slave2master[slaveLoc] : -1;
  979.         // if we're currently inside a block..
  980.         if (importStart >= 0) {
  982.             // add another residue to a continuously aligned block
  983.             if (masterLoc >= 0 && masterLoc-masterStart == importLoc-importStart) {
  984.                 ++len;
  985.             }
  987.             // terminate block
  988.             else {
  989.                 UngappedAlignedBlock *newBlock = new UngappedAlignedBlock(newAlignment);
  990.                 newBlock->SetRangeOfRow(0, masterStart, masterStart + len - 1);
  991.                 newBlock->SetRangeOfRow(1, importStart, importStart + len - 1);
  992.                 newBlock->width = len;
  993.                 newAlignment->AddAlignedBlockAtEnd(newBlock);
  994.                 importStart = -1;
  995.             }
  996.         }
  998.         // search for start of block
  999.         if (importStart < 0) {
  1000.             if (masterLoc >= 0) {
  1001.                 masterStart = masterLoc;
  1002.                 importStart = importLoc;
  1003.                 len = 1;
  1004.             }
  1005.         }
  1006.     }
  1008.     // finalize and and add new alignment to list
  1009.     if (!newAlignment->AddUnalignedBlocks() || !newAlignment->UpdateBlockMapAndColors(false)) {
  1010.         ERRORMSG("error finalizing alignment");
  1011.         delete newAlignment;
  1012.         return NULL;
  1013.     }
  1015.     return newAlignment;
  1016. }
  1018. void UpdateViewer::BlastNeighbor(BlockMultipleAlignment *update)
  1019. {
  1020.     const BlockMultipleAlignment *multiple = alignmentManager->GetCurrentMultipleAlignment();
  1021.     if (!multiple) {
  1022.         ERRORMSG("Can't do BLAST Neighbor when no multiple alignment is present");
  1023.         return;
  1024.     }
  1025.     BlockMultipleAlignment::UngappedAlignedBlockList uaBlocks;
  1026.     multiple->GetUngappedAlignedBlocks(&uaBlocks);
  1027.     if (uaBlocks.size() == 0) {
  1028.         ERRORMSG("Can't do BLAST Neighbor with null multiple alignment");
  1029.         return;
  1030.     }
  1031.     const Sequence *updateSeq = update->GetSequenceOfRow(1);
  1033.     // find alignment, to replace it with BLAST result
  1034.     AlignmentList::iterator a, ae = GetCurrentAlignments().end();
  1035.     for (a=GetCurrentAlignments().begin(); a!=ae; ++a)
  1036.         if (*a == update) break;
  1037.     if (a == GetCurrentAlignments().end()) return;
  1039.     // set up BLAST-2-sequences between update slave and each sequence from the multiple
  1040.     BLASTer::AlignmentList toRealign;
  1041.     for (int row=0; row<multiple->NRows(); ++row) {
  1042.         BlockMultipleAlignment::SequenceList *seqs = new BlockMultipleAlignment::SequenceList(2);
  1043.         (*seqs)[0] = multiple->GetSequenceOfRow(row);
  1044.         (*seqs)[1] = updateSeq;
  1045.         BlockMultipleAlignment *newAlignment =
  1046.             new BlockMultipleAlignment(seqs, updateSeq->parentSet->alignmentManager);
  1047.         if (newAlignment->AddUnalignedBlocks() && newAlignment->UpdateBlockMapAndColors(false))
  1048.         {
  1049.             int excess = 0;
  1050.             if (!RegistryGetInteger(REG_ADVANCED_SECTION, REG_FOOTPRINT_RES, &excess))
  1051.                 WARNINGMSG("Can't get footprint excess residues from registry");
  1052.             newAlignment->alignMasterFrom = uaBlocks.front()->GetRangeOfRow(row)->from - excess;
  1053.             if (newAlignment->alignMasterFrom < 0)
  1054.                 newAlignment->alignMasterFrom = 0;
  1055.             newAlignment->alignMasterTo = uaBlocks.back()->GetRangeOfRow(row)->to + excess;
  1056.             if (newAlignment->alignMasterTo >= (*seqs)[0]->Length())
  1057.                 newAlignment->alignMasterTo = (*seqs)[0]->Length() - 1;
  1058.             newAlignment->alignSlaveFrom = update->alignSlaveFrom;
  1059.             newAlignment->alignSlaveTo = update->alignSlaveTo;
  1060.             toRealign.push_back(newAlignment);
  1061.         } else {
  1062.             ERRORMSG("error finalizing alignment");
  1063.             delete newAlignment;
  1064.         }
  1065.     }
  1067.     // actually do BLAST alignments
  1068.     BLASTer::AlignmentList newAlignments;
  1069.     SetDiagPostLevel(eDiag_Error); // ignore all but Errors while reading data
  1070.     alignmentManager->blaster->
  1071.         CreateNewPairwiseAlignmentsByBlast(NULL, toRealign, &newAlignments, false);
  1072.     SetDiagPostLevel(eDiag_Info);
  1073.     DELETE_ALL_AND_CLEAR(toRealign, BLASTer::AlignmentList);
  1074.     if (newAlignments.size() != multiple->NRows()) {
  1075.         ERRORMSG("UpdateViewer::BlastUpdate() - CreateNewPairwiseAlignmentsByBlast() failed");
  1076.         return;
  1077.     }
  1079.     // replace alignment with result
  1080.     BlockMultipleAlignment *alignmentByNeighbor = GetAlignmentByBestNeighbor(multiple, newAlignments);
  1081.     DELETE_ALL_AND_CLEAR(newAlignments, BLASTer::AlignmentList);
  1082.     if (!alignmentByNeighbor) {
  1083.         WARNINGMSG("alignment unchanged");
  1084.         return;
  1085.     }
  1086.     TRACEMSG("BLAST Neighbor succeeded - replacing alignment");
  1087.     delete update;
  1088.     *a = alignmentByNeighbor;
  1090.     // recreate alignment display with new alignment
  1091.     AlignmentList copy = GetCurrentAlignments();
  1092.     GetCurrentAlignments().clear();
  1093.     GetCurrentDisplay()->Empty();
  1094.     AddAlignments(copy);
  1095. //    (*viewerWindow)->ScrollToColumn(GetCurrentDisplay()->GetStartingColumn());
  1096. }
  1098. // comparison function: if CompareRows(a, b) == true, then row a moves up
  1099. typedef bool (*CompareUpdates)(BlockMultipleAlignment *a, BlockMultipleAlignment *b);
  1101. static bool CompareUpdatesByIdentifier(BlockMultipleAlignment *a, BlockMultipleAlignment *b)
  1102. {
  1103.     return MoleculeIdentifier::CompareIdentifiers(
  1104.         a->GetSequenceOfRow(1)->identifier, // sort by first slave row
  1105.         b->GetSequenceOfRow(1)->identifier);
  1106. }
  1108. static CompareUpdates updateComparisonFunction = NULL;
  1110. void UpdateViewer::SortByIdentifier(void)
  1111. {
  1112.     TRACEMSG("sorting updates by identifier");
  1113.     updateComparisonFunction = CompareUpdatesByIdentifier;
  1114.     SortUpdates();
  1115. }
  1117. void UpdateViewer::SortUpdates(void)
  1118. {
  1119.     if (!updateComparisonFunction) {
  1120.         ERRORMSG("UpdateViewer::SortUpdates() - must first set comparison function");
  1121.         return;
  1122.     }
  1124.     // make vector of alignments
  1125.     AlignmentList& currentAlignments = GetCurrentAlignments();
  1126.     if (currentAlignments.size() < 2) return;
  1127.     vector < BlockMultipleAlignment * > sortedVector(currentAlignments.size());
  1128.     AlignmentList::const_iterator a, ae = currentAlignments.end();
  1129.     int i = 0;
  1130.     for (a=currentAlignments.begin(); a!=ae; ++a) sortedVector[i++] = *a;
  1132.     // sort them
  1133.     stable_sort(sortedVector.begin(), sortedVector.end(), updateComparisonFunction);
  1134.     updateComparisonFunction = NULL;
  1136.     // replace window contents with sorted list
  1137.     currentAlignments.clear();
  1138.     GetCurrentDisplay()->Empty();
  1140.     AlignmentList sortedList;
  1141.     for (i=0; i<sortedVector.size(); ++i) sortedList.push_back(sortedVector[i]);
  1142.     AddAlignments(sortedList);
  1143. }
  1145. END_SCOPE(Cn3D)
  1148. /*
  1149. * ---------------------------------------------------------------------------
  1150. * $Log: update_viewer.cpp,v $
  1151. * Revision 1000.3  2004/06/01 18:29:49  gouriano
  1152. * PRODUCTION: UPGRADED [GCC34_MSVC7] Dev-tree R1.73
  1153. *
  1154. * Revision 1.73  2004/05/21 21:41:40  gorelenk
  1155. * Added PCH ncbi_pch.hpp
  1156. *
  1157. * Revision 1.72  2004/05/21 17:29:51  thiessen
  1158. * allow conversion of mime to cdd data
  1159. *
  1160. * Revision 1.71  2004/03/15 18:38:52  thiessen
  1161. * prefer prefix vs. postfix ++/-- operators
  1162. *
  1163. * Revision 1.70  2004/02/19 17:05:20  thiessen
  1164. * remove cn3d/ from include paths; add pragma to disable annoying msvc warning
  1165. *
  1166. * Revision 1.69  2004/01/17 00:17:32  thiessen
  1167. * add Biostruc and network structure load
  1168. *
  1169. * Revision 1.68  2004/01/05 17:09:16  thiessen
  1170. * abort import and warn if same accession different gi
  1171. *
  1172. * Revision 1.67  2003/11/06 18:52:32  thiessen
  1173. * make geometry violations shown on/off; allow multiple pmid entry in ref dialog
  1174. *
  1175. * Revision 1.66  2003/09/25 14:11:43  thiessen
  1176. * don't assume gi Seq-id's for imported structure's sequences
  1177. *
  1178. * Revision 1.65  2003/08/21 18:27:40  thiessen
  1179. * change header order for Mac compilation
  1180. *
  1181. * Revision 1.64  2003/07/14 18:37:08  thiessen
  1182. * change GetUngappedAlignedBlocks() param types; other syntax changes
  1183. *
  1184. * Revision 1.63  2003/04/04 14:02:22  thiessen
  1185. * more informative error messages on structure import
  1186. *
  1187. * Revision 1.62  2003/04/02 18:03:16  thiessen
  1188. * fix wxString/string confusion
  1189. *
  1190. * Revision 1.61  2003/04/02 17:49:18  thiessen
  1191. * allow pdb id's in structure import dialog
  1192. *
  1193. * Revision 1.60  2003/03/14 19:57:07  thiessen
  1194. * adjust wxStringTokenizer mode
  1195. *
  1196. * Revision 1.59  2003/03/14 19:48:51  thiessen
  1197. * allow multiple gi's in network sequence import dialog
  1198. *
  1199. * Revision 1.58  2003/03/13 18:55:17  thiessen
  1200. * tweak file load error reporting
  1201. *
  1202. * Revision 1.57  2003/02/03 19:20:08  thiessen
  1203. * format changes: move CVS Log to bottom of file, remove std:: from .cpp files, and use new diagnostic macros
  1204. *
  1205. * Revision 1.56  2003/01/31 17:18:58  thiessen
  1206. * many small additions and changes...
  1207. *
  1208. * Revision 1.55  2003/01/23 20:03:05  thiessen
  1209. * add BLAST Neighbor algorithm
  1210. *
  1211. * Revision 1.54  2002/11/19 21:19:44  thiessen
  1212. * more const changes for objects; fix user vs default style bug
  1213. *
  1214. * Revision 1.53  2002/11/06 00:18:10  thiessen
  1215. * fixes for new CRef/const rules in objects
  1216. *
  1217. * Revision 1.52  2002/10/27 22:23:51  thiessen
  1218. * save structure alignments from vastalign.cgi imports
  1219. *
  1220. * Revision 1.51  2002/10/25 19:00:02  thiessen
  1221. * retrieve VAST alignment from vastalign.cgi on structure import
  1222. *
  1223. * Revision 1.50  2002/10/15 22:04:09  thiessen
  1224. * fix geom vltns bug
  1225. *
  1226. * Revision 1.49  2002/10/13 22:58:08  thiessen
  1227. * add redo ability to editor
  1228. *
  1229. * Revision 1.48  2002/09/30 17:13:02  thiessen
  1230. * change structure import to do sequences as well; change cache to hold mimes; change block aligner vocabulary; fix block aligner dialog bugs
  1231. *
  1232. * Revision 1.47  2002/09/26 18:31:24  thiessen
  1233. * allow simultaneous import of multiple chains from single PDB
  1234. *
  1235. * Revision 1.46  2002/09/16 21:24:58  thiessen
  1236. * add block freezing to block aligner
  1237. *
  1238. * Revision 1.45  2002/09/09 13:38:23  thiessen
  1239. * separate save and save-as
  1240. *
  1241. * Revision 1.44  2002/08/15 22:13:18  thiessen
  1242. * update for wx2.3.2+ only; add structure pick dialog; fix MultitextDialog bug
  1243. *
  1244. * Revision 1.43  2002/08/13 20:46:37  thiessen
  1245. * add global block aligner
  1246. *
  1247. * Revision 1.42  2002/07/27 12:29:52  thiessen
  1248. * fix block aligner crash
  1249. *
  1250. * Revision 1.41  2002/07/26 15:28:48  thiessen
  1251. * add Alejandro's block alignment algorithm
  1252. *
  1253. * Revision 1.40  2002/07/26 13:07:01  thiessen
  1254. * fix const object problem
  1255. *
  1256. * Revision 1.39  2002/07/01 23:17:04  thiessen
  1257. * skip warning if master choice canceled
  1258. *
  1259. * Revision 1.38  2002/06/06 01:30:02  thiessen
  1260. * fixes for linux/gcc
  1261. *
  1262. * Revision 1.37  2002/06/05 14:28:41  thiessen
  1263. * reorganize handling of window titles
  1264. *
  1265. * Revision 1.36  2002/06/04 12:48:56  thiessen
  1266. * tweaks for release ; fill out help menu
  1267. *
  1268. * Revision 1.35  2002/05/22 17:17:09  thiessen
  1269. * progress on BLAST interface ; change custom spin ctrl implementation
  1270. *
  1271. * Revision 1.34  2002/05/17 19:10:27  thiessen
  1272. * preliminary range restriction for BLAST/PSSM
  1273. *
  1274. * Revision 1.33  2002/04/26 19:01:00  thiessen
  1275. * fix display delete bug
  1276. *
  1277. * Revision 1.32  2002/03/28 14:06:02  thiessen
  1278. * preliminary BLAST/PSSM ; new CD startup style
  1279. *
  1280. * Revision 1.31  2002/03/07 19:16:04  thiessen
  1281. * don't auto-show sequence windows
  1282. *
  1283. * Revision 1.30  2002/03/04 15:52:15  thiessen
  1284. * hide sequence windows instead of destroying ; add perspective/orthographic projection choice
  1285. *
  1286. * Revision 1.29  2002/02/27 16:29:42  thiessen
  1287. * add model type flag to general mime type
  1288. *
  1289. * Revision 1.28  2002/02/22 14:24:01  thiessen
  1290. * sort sequences in reject dialog ; general identifier comparison
  1291. *
  1292. * Revision 1.27  2002/02/13 14:53:30  thiessen
  1293. * add update sort
  1294. *
  1295. * Revision 1.26  2002/02/12 17:19:22  thiessen
  1296. * first working structure import
  1297. *
  1298. * Revision 1.25  2002/02/01 00:41:21  thiessen
  1299. * tweaks
  1300. *
  1301. * Revision 1.24  2002/01/24 20:07:57  thiessen
  1302. * read multiple FAST sequences
  1303. *
  1304. * Revision 1.23  2002/01/02 02:08:29  thiessen
  1305. * go back to viewer.cgi to test http/301 forwarding
  1306. *
  1307. * Revision 1.22  2001/12/15 03:15:59  thiessen
  1308. * adjustments for slightly changed object loader Set...() API
  1309. *
  1310. * Revision 1.21  2001/12/12 14:58:10  thiessen
  1311. * change URL to viewer.fcgi
  1312. *
  1313. * Revision 1.20  2001/12/06 23:13:47  thiessen
  1314. * finish import/align new sequences into single-structure data; many small tweaks
  1315. *
  1316. * Revision 1.19  2001/11/30 14:02:05  thiessen
  1317. * progress on sequence imports to single structures
  1318. *
  1319. * Revision 1.18  2001/11/27 16:26:10  thiessen
  1320. * major update to data management system
  1321. *
  1322. * Revision 1.17  2001/10/01 16:04:25  thiessen
  1323. * make CDD annotation window non-modal; add SetWindowTitle to viewers
  1324. *
  1325. * Revision 1.16  2001/09/27 15:38:00  thiessen
  1326. * decouple sequence import and BLAST
  1327. *
  1328. * Revision 1.15  2001/09/20 19:31:30  thiessen
  1329. * fixes for SGI and wxWin 2.3.2
  1330. *
  1331. * Revision 1.14  2001/09/19 22:55:39  thiessen
  1332. * add preliminary net import and BLAST
  1333. *
  1334. * Revision 1.13  2001/09/18 03:10:46  thiessen
  1335. * add preliminary sequence import pipeline
  1336. *
  1337. * Revision 1.12  2001/06/02 17:22:46  thiessen
  1338. * fixes for GCC
  1339. *
  1340. * Revision 1.11  2001/05/17 18:34:26  thiessen
  1341. * spelling fixes; change dialogs to inherit from wxDialog
  1342. *
  1343. * Revision 1.10  2001/05/02 13:46:28  thiessen
  1344. * major revision of stuff relating to saving of updates; allow stored null-alignments
  1345. *
  1346. * Revision 1.9  2001/04/20 18:02:41  thiessen
  1347. * don't open update viewer right away
  1348. *
  1349. * Revision 1.8  2001/04/19 12:58:32  thiessen
  1350. * allow merge and delete of individual updates
  1351. *
  1352. * Revision 1.7  2001/04/05 22:55:36  thiessen
  1353. * change bg color handling ; show geometry violations
  1354. *
  1355. * Revision 1.6  2001/04/04 00:27:15  thiessen
  1356. * major update - add merging, threader GUI controls
  1357. *
  1358. * Revision 1.5  2001/03/30 03:07:34  thiessen
  1359. * add threader score calculation & sorting
  1360. *
  1361. * Revision 1.4  2001/03/22 00:33:17  thiessen
  1362. * initial threading working (PSSM only); free color storage in undo stack
  1363. *
  1364. * Revision 1.3  2001/03/17 14:06:49  thiessen
  1365. * more workarounds for namespace/#define conflicts
  1366. *
  1367. * Revision 1.2  2001/03/13 01:25:06  thiessen
  1368. * working undo system for >1 alignment (e.g., update window)
  1369. *
  1370. * Revision 1.1  2001/03/09 15:49:06  thiessen
  1371. * major changes to add initial update viewer
  1372. *
  1373. */