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

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

生物技术

开发平台：

C/C++

structure_set.cpp：源码内容
  1. /*
  2.  * ===========================================================================
  3.  * PRODUCTION $Log: structure_set.cpp,v $
  4.  * PRODUCTION Revision 1000.2  2004/06/01 18:29:28  gouriano
  5.  * PRODUCTION PRODUCTION: UPGRADED [GCC34_MSVC7] Dev-tree R1.141
  6.  * PRODUCTION
  7.  * ===========================================================================
  8.  */
  10. /*  $Id: structure_set.cpp,v 1000.2 2004/06/01 18:29:28 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. *      Classes to hold sets of structure data
  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/ncbistre.hpp>
  50. #include <corelib/ncbi_limits.h>
  51. #include <corelib/ncbistl.hpp>
  53. #include <deque>
  55. #include <objects/ncbimime/Biostruc_seq.hpp>
  56. #include <objects/ncbimime/Biostruc_seqs.hpp>
  57. #include <objects/ncbimime/Biostruc_align.hpp>
  58. #include <objects/ncbimime/Entrez_general.hpp>
  59. #include <objects/mmdb1/Biostruc_id.hpp>
  60. #include <objects/mmdb1/Mmdb_id.hpp>
  61. #include <objects/mmdb1/Biostruc_descr.hpp>
  62. #include <objects/mmdb2/Biostruc_model.hpp>
  63. #include <objects/mmdb2/Model_type.hpp>
  64. #include <objects/mmdb3/Biostruc_feature.hpp>
  65. #include <objects/mmdb3/Biostruc_feature_id.hpp>
  66. #include <objects/mmdb3/Chem_graph_alignment.hpp>
  67. #include <objects/mmdb3/Transform.hpp>
  68. #include <objects/mmdb3/Move.hpp>
  69. #include <objects/mmdb3/Trans_matrix.hpp>
  70. #include <objects/mmdb3/Rot_matrix.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/seqset/Bioseq_set.hpp>
  75. #include <objects/cn3d/Cn3d_style_dictionary.hpp>
  76. #include <objects/cn3d/Cn3d_user_annotations.hpp>
  77. #include <objects/seqalign/Dense_diag.hpp>
  78. #include <objects/seqalign/Dense_seg.hpp>
  79. #include <objects/mmdb3/Biostruc_feature_set_id.hpp>
  80. #include <objects/mmdb1/Molecule_id.hpp>
  81. #include <objects/mmdb1/Residue_id.hpp>
  82. #include <objects/cdd/Reject_id.hpp>
  83. #include <objects/cdd/Update_comment.hpp>
  85. #include "structure_set.hpp"
  86. #include "data_manager.hpp"
  87. #include "coord_set.hpp"
  88. #include "chemical_graph.hpp"
  89. #include "atom_set.hpp"
  90. #include "opengl_renderer.hpp"
  91. #include "show_hide_manager.hpp"
  92. #include "style_manager.hpp"
  93. #include "sequence_set.hpp"
  94. #include "alignment_set.hpp"
  95. #include "alignment_manager.hpp"
  96. #include "messenger.hpp"
  97. #include "asn_converter.hpp"
  98. #include "block_multiple_alignment.hpp"
  99. #include "cn3d_tools.hpp"
  100. #include "molecule_identifier.hpp"
  101. #include "cn3d_cache.hpp"
  102. #include "molecule.hpp"
  103. #include "residue.hpp"
  104. #include "show_hide_dialog.hpp"
  106. #include <objseq.h>
  108. USING_NCBI_SCOPE;
  109. USING_SCOPE(objects);
  112. BEGIN_SCOPE(Cn3D)
  114. const unsigned int
  115.     StructureSet::ePSSMData                   = 0x01,
  116.     StructureSet::eRowOrderData               = 0x02,
  117.     StructureSet::eAnyAlignmentData           = 0x04,
  118.     StructureSet::eStructureAlignmentData     = 0x08,
  119.     StructureSet::eSequenceData               = 0x10,
  120.     StructureSet::eUpdateData                 = 0x20,
  121.     StructureSet::eStyleData                  = 0x40,
  122.     StructureSet::eUserAnnotationData         = 0x80,
  123.     StructureSet::eCDDData                    = 0x100,
  124.     StructureSet::eOtherData                  = 0x200;
  126. StructureSet::StructureSet(CNcbi_mime_asn1 *mime, int structureLimit, OpenGLRenderer *r) :
  127.     StructureBase(NULL), renderer(r)
  128. {
  129.     dataManager = new ASNDataManager(mime);
  130.     Load(structureLimit);
  131. }
  133. StructureSet::StructureSet(CCdd *cdd, int structureLimit, OpenGLRenderer *r) :
  134.     StructureBase(NULL), renderer(r)
  135. {
  136.     dataManager = new ASNDataManager(cdd);
  137.     Load(structureLimit);
  138. }
  140. void StructureSet::LoadSequencesForSingleStructure(void)
  141. {
  142.     sequenceSet = new SequenceSet(this, *(dataManager->GetSequences()));
  144.     if (objects.size() > 1)
  145.         ERRORMSG("LoadSequencesForSingleStructure() called, but there is > 1 structure");
  146.     if (objects.size() != 1) return;
  148.     // look for biopolymer molecules
  149.     ChemicalGraph::MoleculeMap::const_iterator m, me = objects.front()->graph->molecules.end();
  150.     SequenceSet::SequenceList::const_iterator s, se = sequenceSet->sequences.end();
  151.     for (m=objects.front()->graph->molecules.begin(); m!=me; ++m) {
  152.         if (!m->second->IsProtein() && !m->second->IsNucleotide()) continue;
  154.         // find matching sequence for each biopolymer
  155.         for (s=sequenceSet->sequences.begin(); s!=se; ++s) {
  156.             if ((*s)->molecule != NULL) continue; // skip already-matched sequences
  158.             if (m->second->identifier == (*s)->identifier) {
  160.                 // verify length
  161.                 if (m->second->residues.size() != (*s)->Length()) {
  162.                     ERRORMSG(
  163.                         "LoadSequencesForSingleStructure() - length mismatch between sequence gi "
  164.                         << "and matching molecule " << m->second->identifier->ToString());
  165.                     continue;
  166.                 }
  167.                 TRACEMSG("matched sequence " << " gi " << (*s)->identifier->gi << " with object "
  168.                     << objects.front()->pdbID << " moleculeID " << m->second->id);
  170.                 (const_cast<Molecule*>(m->second))->sequence = *s;
  171.                 (const_cast<Sequence*>(*s))->molecule = m->second;
  172.                 break;
  173.             }
  174.         }
  175.         if (s == se)
  176.             ERRORMSG("LoadSequencesForSingleStructure() - can't find sequence for molecule "
  177.                 << m->second->identifier->ToString());
  178.     }
  179. }
  181. bool StructureSet::LoadMaster(int masterMMDBID)
  182. {
  183.     if (objects.size() > 0) return false;
  184.     TRACEMSG("loading master " << masterMMDBID);
  186.     if (dataManager->GetMasterStructure()) {
  187.         objects.push_back(new StructureObject(this, *(dataManager->GetMasterStructure()), true));
  188.         if (masterMMDBID != MoleculeIdentifier::VALUE_NOT_SET && objects.front()->mmdbID != masterMMDBID)
  189.             ERRORMSG("StructureSet::LoadMaster() - mismatched master MMDB ID");
  190.     } else if (masterMMDBID != MoleculeIdentifier::VALUE_NOT_SET && dataManager->GetStructureList()) {
  191.         ASNDataManager::BiostrucList::const_iterator b, be = dataManager->GetStructureList()->end();
  192.         for (b=dataManager->GetStructureList()->begin(); b!=be; ++b) {
  193.             if ((*b)->GetId().front()->IsMmdb_id() &&
  194.                 (*b)->GetId().front()->GetMmdb_id().Get() == masterMMDBID) {
  195.                 objects.push_back(new StructureObject(this, **b, true));
  196.                 usedStructures[b->GetPointer()] = true;
  197.                 break;
  198.             }
  199.         }
  200.     }
  201.     if (masterMMDBID != MoleculeIdentifier::VALUE_NOT_SET && objects.size() == 0) {
  202.         CRef < CBiostruc > biostruc;
  203.         wxString id;
  204.         id.Printf("%i", masterMMDBID);
  205.         if (LoadStructureViaCache(id.c_str(), dataManager->GetBiostrucModelType(), biostruc, NULL))
  206.             objects.push_back(new StructureObject(this, *biostruc, true));
  207.     }
  208.     return (objects.size() > 0);
  209. }
  211. bool StructureSet::MatchSequenceToMoleculeInObject(const Sequence *seq,
  212.     const StructureObject *obj, const Sequence **seqHandle)
  213. {
  214.     ChemicalGraph::MoleculeMap::const_iterator m, me = obj->graph->molecules.end();
  215.     for (m=obj->graph->molecules.begin(); m!=me; ++m) {
  216.         if (!(m->second->IsProtein() || m->second->IsNucleotide())) continue;
  218.         if (m->second->identifier == seq->identifier) {
  220.             // verify length
  221.             if (m->second->residues.size() != seq->Length()) {
  222.                 ERRORMSG(
  223.                     "MatchSequenceToMoleculeInObject() - length mismatch between sequence gi "
  224.                     << "and matching molecule " << m->second->identifier->ToString());
  225.                 continue;
  226.             }
  227.             TRACEMSG("matched sequence " << " gi " << seq->identifier->gi << " with object "
  228.                 << obj->pdbID << " moleculeID " << m->second->id);
  230.             // sanity check
  231.             if (m->second->sequence) {
  232.                 ERRORMSG("Molecule " << m->second->identifier->ToString()
  233.                     << " already has an associated sequence");
  234.                 continue;
  235.             }
  237.             // automatically duplicate Sequence if it's already associated with a molecule
  238.             if (seq->molecule) {
  239.                 TRACEMSG("duplicating sequence " << seq->identifier->ToString());
  240. SequenceSet *seqSetMod = const_cast<SequenceSet*>(sequenceSet);
  241.                 CBioseq& bioseqMod = const_cast<CBioseq&>(seq->bioseqASN.GetObject());
  242.                 seq = new Sequence(seqSetMod, bioseqMod);
  243.                 seqSetMod->sequences.push_back(seq);
  244.                 // update Sequence handle, which should be a handle to a MasterSlaveAlignment slave,
  245.                 // so that this new Sequence* is correctly loaded into the BlockMultipleAlignment
  246.                 if (seqHandle) *seqHandle = seq;
  247.             }
  249.             // do the cross-match
  250.             (const_cast<Molecule*>(m->second))->sequence = seq;
  251.             (const_cast<Sequence*>(seq))->molecule = m->second;
  252.             break;
  253.         }
  254.     }
  255.     return (m != me);
  256. }
  258. static void SetStructureRowFlags(const AlignmentSet *alignmentSet, int *structureLimit,
  259.     vector < bool > *dontLoadRowStructure)
  260. {
  261.     vector < string > titles;
  262.     vector < int > rows;
  264.     // find slave rows with associated structure
  265.     AlignmentSet::AlignmentList::const_iterator l, le = alignmentSet->alignments.end();
  266.     int row;
  267.     for (l=alignmentSet->alignments.begin(), row=0; l!=le; ++l, ++row) {
  268.         if ((*l)->slave->identifier->mmdbID != MoleculeIdentifier::VALUE_NOT_SET) {
  269.             titles.push_back((*l)->slave->identifier->ToString());
  270.             rows.push_back(row);
  271.         }
  272.     }
  274.     if (*structureLimit - 1 >= titles.size()) return;
  276.     // let user select which slaves to load
  277.     wxString *items = new wxString[titles.size()];
  278.     vector < bool > itemsOn(titles.size(), false);
  279.     for (row=0; row<titles.size(); ++row) {
  280.         items[row] = titles[row].c_str();
  281.         if (row < *structureLimit - 1)      // by default, first N-1 are selected
  282.             itemsOn[row] = true;
  283.     }
  284.     ShowHideDialog dialog(items, &itemsOn, NULL, true, NULL, -1, "Choose structures:");
  286.     if (dialog.ShowModal() == wxOK) {
  287.         // figure out which rows the user selected, and adjust structureLimit accordingly
  288.         *structureLimit = 1;     // master always visible
  289.         for (row=0; row<itemsOn.size(); ++row) {
  290.             if (itemsOn[row])
  291.                 (*structureLimit)++;                        // structure should be loaded
  292.             else
  293.                 (*dontLoadRowStructure)[rows[row]] = true;  // structure should not be loaded
  294.         }
  295.     }
  297.     delete[] items;
  298. }
  300. void StructureSet::LoadAlignmentsAndStructures(int structureLimit)
  301. {
  302.     // try to determine the master structure
  303.     int masterMMDBID = MoleculeIdentifier::VALUE_NOT_SET;
  304.     // explicit master structure
  305.     if (dataManager->GetMasterStructure() &&
  306.         dataManager->GetMasterStructure()->GetId().front()->IsMmdb_id())
  307.         masterMMDBID = dataManager->GetMasterStructure()->GetId().front()->GetMmdb_id().Get();
  308.     // master of structure alignments
  309.     else if (dataManager->GetStructureAlignments() &&
  310.                 dataManager->GetStructureAlignments()->IsSetId() &&
  311.                 dataManager->GetStructureAlignments()->GetId().front()->IsMmdb_id())
  312.         masterMMDBID = dataManager->GetStructureAlignments()->GetId().front()->GetMmdb_id().Get();
  313.     // SPECIAL CASE: if there's no explicit master structure, but there is a structure
  314.     // list that contains a single structure, then assume that it is the master structure
  315.     else if (dataManager->GetStructureList() &&
  316.                 dataManager->GetStructureList()->size() == 1 &&
  317.                 dataManager->GetStructureList()->front()->GetId().front()->IsMmdb_id())
  318.         masterMMDBID = dataManager->GetStructureList()->front()->GetId().front()->GetMmdb_id().Get();
  320.     // assume data manager has already screened the alignment list
  321.     ASNDataManager::SeqAnnotList *alignments = dataManager->GetSequenceAlignments();
  322.     typedef list < CRef < CSeq_align > > SeqAlignList;
  323.     const SeqAlignList& seqAligns = alignments->front()->GetData().GetAlign();
  325.     // we need to determine the identity of the master sequence; most rigorous way is to look
  326.     // for a Seq-id that is present in all pairwise alignments
  327.     const Sequence *seq1 = NULL, *seq2 = NULL, *master = NULL;
  328.     bool seq1PresentInAll = true, seq2PresentInAll = true;
  330.     // first, find sequences for first pairwise alignment
  331.     const CSeq_id& frontSid = seqAligns.front()->GetSegs().IsDendiag() ?
  332.         seqAligns.front()->GetSegs().GetDendiag().front()->GetIds().front().GetObject() :
  333.         seqAligns.front()->GetSegs().GetDenseg().GetIds().front().GetObject();
  334.     const CSeq_id& backSid = seqAligns.front()->GetSegs().IsDendiag() ?
  335.         seqAligns.front()->GetSegs().GetDendiag().front()->GetIds().back().GetObject() :
  336.         seqAligns.front()->GetSegs().GetDenseg().GetIds().back().GetObject();
  337.     SequenceSet::SequenceList::const_iterator s, se = sequenceSet->sequences.end();
  338.     for (s=sequenceSet->sequences.begin(); s!=se; ++s) {
  339.         if ((*s)->identifier->MatchesSeqId(frontSid)) seq1 = *s;
  340.         if ((*s)->identifier->MatchesSeqId(backSid)) seq2 = *s;
  341.         if (seq1 && seq2) break;
  342.     }
  343.     if (!(seq1 && seq2)) {
  344.         ERRORMSG("Can't match first pair of Seq-ids to Sequences");
  345.         return;
  346.     }
  348.     // now, make sure one of these sequences is present in all the other pairwise alignments
  349.     SeqAlignList::const_iterator a = seqAligns.begin(), ae = seqAligns.end();
  350.     for (++a; a!=ae; ++a) {
  351.         const CSeq_id& frontSid2 = (*a)->GetSegs().IsDendiag() ?
  352.             (*a)->GetSegs().GetDendiag().front()->GetIds().front().GetObject() :
  353.             (*a)->GetSegs().GetDenseg().GetIds().front().GetObject();
  354.         const CSeq_id& backSid2 = (*a)->GetSegs().IsDendiag() ?
  355.             (*a)->GetSegs().GetDendiag().front()->GetIds().back().GetObject() :
  356.             (*a)->GetSegs().GetDenseg().GetIds().back().GetObject();
  357.         if (!seq1->identifier->MatchesSeqId(frontSid2) && !seq1->identifier->MatchesSeqId(backSid2))
  358.             seq1PresentInAll = false;
  359.         if (!seq2->identifier->MatchesSeqId(frontSid2) && !seq2->identifier->MatchesSeqId(backSid2))
  360.             seq2PresentInAll = false;
  361.     }
  362.     if (!seq1PresentInAll && !seq2PresentInAll) {
  363.         ERRORMSG("All pairwise sequence alignments must have a common master sequence");
  364.         return;
  365.     } else if (seq1PresentInAll && !seq2PresentInAll)
  366.         master = seq1;
  367.     else if (seq2PresentInAll && !seq1PresentInAll)
  368.         master = seq2;
  369.     else if (seq1PresentInAll && seq2PresentInAll && seq1 == seq2)
  370.         master = seq1;
  372.     // if still ambiguous, see if master3d is set in CDD data
  373.     if (!master && dataManager->GetCDDMaster3d()) {
  374.         if (seq1->identifier->MatchesSeqId(*(dataManager->GetCDDMaster3d())))
  375.             master = seq1;
  376.         else if (seq2->identifier->MatchesSeqId(*(dataManager->GetCDDMaster3d())))
  377.             master = seq2;
  378.         else
  379.             ERRORMSG("Unable to match CDD's master3d with either sequence in first pairwise alignment");
  380.     }
  382.     // if still ambiguous, try to use master structure info to find master sequence
  383.     if (!master && masterMMDBID != MoleculeIdentifier::VALUE_NOT_SET) {
  384.         // load master - has side affect of matching gi's with PDB/molecule ID during graph evaluation
  385.         if (structureLimit > 0)
  386.             LoadMaster(masterMMDBID);
  388.         // see if there's a sequence in the master structure that matches
  389.         if (seq1->identifier->mmdbID != MoleculeIdentifier::VALUE_NOT_SET &&
  390.             seq1->identifier->mmdbID != seq2->identifier->mmdbID) {
  391.             if (masterMMDBID == seq1->identifier->mmdbID)
  392.                 master = seq1;
  393.             else if (masterMMDBID == seq2->identifier->mmdbID)
  394.                 master = seq2;
  395.             else {
  396.                 ERRORMSG("Structure master does not contain either sequence in first pairwise alignment");
  397.                 return;
  398.             }
  399.         }
  400.     }
  402.     // if still ambiguous, just use the first one
  403.     if (!master) {
  404.         WARNINGMSG("Ambiguous master; using " << seq1->identifier->ToString());
  405.         master = seq1;
  406.     }
  408.     TRACEMSG("determined that master sequence is " << master->identifier->ToString());
  410.     // load alignments now that we know the identity of the master
  411.     alignmentSet = new AlignmentSet(this, master, *(dataManager->GetSequenceAlignments()));
  413.     // check mmdb id's and load master if not already present (and if master has structure)
  414.     if (masterMMDBID == MoleculeIdentifier::VALUE_NOT_SET) {
  415.         masterMMDBID = master->identifier->mmdbID;
  416.     } else if (master->identifier->mmdbID != MoleculeIdentifier::VALUE_NOT_SET &&
  417.                master->identifier->mmdbID != masterMMDBID) {
  418.         ERRORMSG("master structure (" << masterMMDBID <<
  419.             ") disagrees with master sequence (apparently from " << master->identifier->mmdbID << ')');
  420.         return;
  421.     }
  422.     if (objects.size() == 0 && structureLimit > 0 && masterMMDBID != MoleculeIdentifier::VALUE_NOT_SET &&
  423.             master->identifier->accession != "consensus")   // special case for looking at "raw" CD's
  424.         LoadMaster(masterMMDBID);
  426.     // cross-match master sequence and structure
  427.     if (objects.size() == 1 && !MatchSequenceToMoleculeInObject(master, objects.front())) {
  428.         ERRORMSG("MatchSequenceToMoleculeInObject() - can't find molecule in object "
  429.             << objects.front()->pdbID << " to match master sequence "
  430.             << master->identifier->ToString());
  431.         return;
  432.     }
  434.     // IFF there's a master structure, then also load slave structures and cross-match sequences
  435.     if (objects.size() == 1 && structureLimit > 1) {
  436.         ASNDataManager::BiostrucList::const_iterator b, be;
  437.         if (dataManager->GetStructureList()) be = dataManager->GetStructureList()->end();
  438.         int row;
  439.         vector < bool > loadedStructureForSlaveRow(alignmentSet->alignments.size(), false);
  441.         // first, load each remaining slave structure, and for each one, find the first slave
  442.         // sequence that matches it (and that doesn't already have structure)
  443.         AlignmentSet::AlignmentList::const_iterator l, le = alignmentSet->alignments.end();
  444.         if (dataManager->GetStructureList()) {
  445.             for (b=dataManager->GetStructureList()->begin(); b!=be && objects.size()<structureLimit; ++b) {
  447.                 // load structure
  448.                 if (usedStructures.find(b->GetPointer()) != usedStructures.end()) continue;
  449.                 StructureObject *object = new StructureObject(this, **b, false);
  450.                 objects.push_back(object);
  451.                 if (dataManager->GetStructureAlignments())
  452.                     object->SetTransformToMaster(
  453.                         *(dataManager->GetStructureAlignments()), master->identifier->mmdbID);
  454.                 usedStructures[b->GetPointer()] = true;
  456.                 // find matching unstructured slave sequence
  457.                 for (l=alignmentSet->alignments.begin(), row=0; l!=le; ++l, ++row) {
  458.                     if (loadedStructureForSlaveRow[row]) continue;
  459.                     if (MatchSequenceToMoleculeInObject((*l)->slave, object,
  460.                             &((const_cast<MasterSlaveAlignment*>(*l))->slave))) {
  461.                         loadedStructureForSlaveRow[row] = true;
  462.                         break;
  463.                     }
  464.                 }
  465.                 if (l == le)
  466.                     ERRORMSG("Warning: Structure " << object->pdbID
  467.                         << " doesn't have a matching slave sequence in the multiple alignment");
  468.             }
  469.         }
  471.         // now loop through slave rows of the alignment; if the slave
  472.         // sequence has an MMDB ID but no structure yet, then load it.
  473.         if (objects.size() < structureLimit && (dataManager->IsCDD() || dataManager->IsGeneralMime())) {
  475.             // for CDD's, ask user which structures to load if structureLimit is low
  476.             if (dataManager->IsCDD())
  477.                 SetStructureRowFlags(alignmentSet, &structureLimit, &loadedStructureForSlaveRow);
  479.             for (l=alignmentSet->alignments.begin(), row=0; l!=le && objects.size()<structureLimit; ++l, ++row) {
  481.                 if ((*l)->slave->identifier->mmdbID != MoleculeIdentifier::VALUE_NOT_SET &&
  482.                     !loadedStructureForSlaveRow[row]) {
  484.                     // first check the biostruc list to see if this structure is present already
  485.                     CRef < CBiostruc > biostruc;
  486.                     if (dataManager->GetStructureList()) {
  487.                         for (b=dataManager->GetStructureList()->begin(); b!=be ; ++b) {
  488.                             if ((*b)->GetId().front()->IsMmdb_id() &&
  489.                                 (*b)->GetId().front()->GetMmdb_id().Get() == (*l)->slave->identifier->mmdbID) {
  490.                                 biostruc = *b;
  491.                                 break;
  492.                             }
  493.                         }
  494.                     }
  496.                     // if not in list, load Biostruc via HTTP/cache
  497.                     if (biostruc.Empty()) {
  498.                         wxString id;
  499.                         id.Printf("%i", (*l)->slave->identifier->mmdbID);
  500.                         if (!LoadStructureViaCache(id.c_str(),
  501.                                 dataManager->GetBiostrucModelType(), biostruc, NULL)) {
  502.                             ERRORMSG("Failed to load MMDB #" << (*l)->slave->identifier->mmdbID);
  503.                             continue;
  504.                         }
  505.                     }
  507.                     // create StructureObject and cross-match
  508.                     StructureObject *object = new StructureObject(this, *biostruc, false);
  509.                     objects.push_back(object);
  510.                     if (dataManager->GetStructureAlignments())
  511.                         object->SetTransformToMaster(
  512.                             *(dataManager->GetStructureAlignments()), master->identifier->mmdbID);
  513.                     if (!MatchSequenceToMoleculeInObject((*l)->slave, object,
  514.                             &((const_cast<MasterSlaveAlignment*>(*l))->slave)))
  515.                         ERRORMSG("Failed to match any molecule in structure " << object->pdbID
  516.                             << " with sequence " << (*l)->slave->identifier->ToString());
  517.                     loadedStructureForSlaveRow[row] = true;
  518.                 }
  519.             }
  520.         }
  521.     }
  522. }
  524. void StructureSet::Load(int structureLimit)
  525. {
  526.     // member data initialization
  527.     lastAtomName = OpenGLRenderer::NO_NAME;
  528.     lastDisplayList = OpenGLRenderer::NO_LIST;
  529.     sequenceSet = NULL;
  530.     alignmentSet = NULL;
  531.     alignmentManager = NULL;
  532.     nDomains = 0;
  533.     isAlphaOnly = false;
  534.     parentSet = this;
  535.     showHideManager = new ShowHideManager();
  536.     styleManager = new StyleManager(this);
  537.     havePrevPickedAtomCoord = false;
  538.     hasUserStyle = false;
  540.     // if this is a single structure, then there should be one sequence per biopolymer
  541.     if (dataManager->IsSingleStructure()) {
  542.         const CBiostruc *masterBiostruc = dataManager->GetMasterStructure();
  543.         if (!masterBiostruc && dataManager->GetStructureList() && dataManager->GetStructureList()->size() == 1)
  544.             masterBiostruc = dataManager->GetStructureList()->front().GetPointer();
  545.         if (masterBiostruc)
  546.             objects.push_back(new StructureObject(this, *masterBiostruc, true));
  547.         if (dataManager->GetSequences())
  548.             LoadSequencesForSingleStructure();
  549.     }
  551.     // multiple structure: should have exactly one sequence per structure (plus unstructured sequences)
  552.     else {
  553.         if (!dataManager->GetSequences() || !dataManager->GetSequenceAlignments()) {
  554.             ERRORMSG("Data interpreted as multiple alignment, "
  555.                 "but missing sequences and/or sequence alignments");
  556.             return;
  557.         }
  558.         sequenceSet = new SequenceSet(this, *(dataManager->GetSequences()));
  559.         LoadAlignmentsAndStructures(structureLimit);
  560.     }
  562.     // find center of coordinates
  563.     SetCenter();
  565.     // create alignment manager
  566.     if (sequenceSet) {
  567.         if (dataManager->GetUpdates())
  568.             // if updates present, alignment manager will load those into update viewer
  569.             alignmentManager = new AlignmentManager(sequenceSet, alignmentSet, *(dataManager->GetUpdates()));
  570.         else
  571.             alignmentManager = new AlignmentManager(sequenceSet, alignmentSet);
  572.     }
  574. #ifdef _DEBUG
  575.     VerifyFrameMap();
  576. #endif
  578.     // setup show/hide items
  579.     showHideManager->ConstructShowHideArray(this);
  580.     // alignments always start with aligned domains only
  581.     if (alignmentSet)
  582.         showHideManager->ShowAlignedDomains(this);
  584.     // load style dictionary and user annotations
  585.     const CCn3d_style_dictionary *styles = dataManager->GetStyleDictionary();
  586.     if (styles) {
  587.         if (!styleManager->LoadFromASNStyleDictionary(*styles) ||
  588.             !styleManager->CheckGlobalStyleSettings())
  589.             ERRORMSG("Error loading style dictionary");
  590.         dataManager->RemoveStyleDictionary();   // remove now; recreated with current settings upon save
  591.         hasUserStyle = true;
  592.     }
  594.     const CCn3d_user_annotations *annots = dataManager->GetUserAnnotations();
  595.     if (annots) {
  596.         if (!styleManager->LoadFromASNUserAnnotations(*annots) ||
  597.             !renderer->LoadFromASNViewSettings(*annots))
  598.             ERRORMSG("Error loading user annotations or camera settings");
  599.         dataManager->RemoveUserAnnotations();   // remove now; recreated with current settings upon save
  600.     }
  602.     dataManager->SetDataUnchanged();
  603. }
  605. StructureSet::~StructureSet(void)
  606. {
  607.     delete dataManager;
  608.     delete showHideManager;
  609.     delete styleManager;
  610.     if (alignmentManager) delete alignmentManager;
  612.     GlobalMessenger()->RemoveAllHighlights(false);
  613.     MoleculeIdentifier::ClearIdentifiers();
  615.     BioseqMap::iterator i, ie = bioseqs.end();
  616.     for (i=bioseqs.begin(); i!=ie; ++i) BioseqFree(i->second);
  617. }
  619. BioseqPtr StructureSet::GetOrCreateBioseq(const Sequence *sequence)
  620. {
  621.     if (!sequence || !sequence->isProtein) {
  622.         ERRORMSG("StructureSet::GetOrCreateBioseq() - got non-protein or NULL Sequence");
  623.         return NULL;
  624.     }
  626.     // if already done
  627.     BioseqMap::const_iterator b = bioseqs.find(sequence);
  628.     if (b != bioseqs.end()) return b->second;
  630.     // create new Bioseq and fill it in from Sequence data
  631.     BioseqPtr bioseq = BioseqNew();
  632.     bioseq->mol = Seq_mol_aa;
  633.     bioseq->seq_data_type = Seq_code_ncbieaa;
  634.     bioseq->repr = Seq_repr_raw;
  635.     bioseq->length = sequence->Length();
  636.     bioseq->seq_data = BSNew(bioseq->length);
  637.     BSWrite(bioseq->seq_data, const_cast<char*>(sequence->sequenceString.c_str()), bioseq->length);
  639.     // create Seq-id
  640.     sequence->AddCSeqId(&(bioseq->id), true);
  642.     // store Bioseq
  643.     bioseqs[sequence] = bioseq;
  645.     return bioseq;
  646. }
  648. void StructureSet::CreateAllBioseqs(const BlockMultipleAlignment *multiple)
  649. {
  650.     for (int row=0; row<multiple->NRows(); ++row)
  651.         GetOrCreateBioseq(multiple->GetSequenceOfRow(row));
  652. }
  654. bool StructureSet::AddBiostrucToASN(ncbi::objects::CBiostruc *biostruc)
  655. {
  656.     bool added = dataManager->AddBiostrucToASN(biostruc);
  657.     if (added && objects.size() == 1)
  658.         InitStructureAlignments(objects.front()->mmdbID);
  659.     return added;
  660. }
  662. static const int NO_DOMAIN = -1, MULTI_DOMAIN = 0;
  664. void StructureSet::InitStructureAlignments(int masterMMDBID)
  665. {
  666.     // create or empty the Biostruc-annot-set that will contain these alignments
  667.     // in the asn data, erasing any structure alignments currently stored there
  668.     CBiostruc_annot_set *structureAlignments = dataManager->GetStructureAlignments();
  669.     if (structureAlignments) {
  670.         structureAlignments->SetId().clear();
  671.         structureAlignments->SetDescr().clear();
  672.         structureAlignments->SetFeatures().clear();
  673.     } else {
  674.         structureAlignments = new CBiostruc_annot_set();
  675.         dataManager->SetStructureAlignments(structureAlignments);
  676.     }
  678.     // set up the skeleton of the new Biostruc-annot-set
  679.     // new Mmdb-id
  680.     structureAlignments->SetId().resize(1);
  681.     structureAlignments->SetId().front().Reset(new CBiostruc_id());
  682.     CMmdb_id *mid = new CMmdb_id(masterMMDBID);
  683.     structureAlignments->SetId().front().GetObject().SetMmdb_id(*mid);
  684.     // new Biostruc-feature-set
  685.     CRef<CBiostruc_feature_set> featSet(new CBiostruc_feature_set());
  686.     featSet->SetId().Set(NO_DOMAIN);
  687.     featSet->SetFeatures();    // just create an empty list
  688.     structureAlignments->SetFeatures().resize(1, featSet);
  690.     // flag a change in data
  691.     SetDataChanged(eStructureAlignmentData);
  692. }
  694. void StructureSet::AddStructureAlignment(CBiostruc_feature *feature,
  695.     int masterDomainID, int slaveDomainID)
  696. {
  697.     CBiostruc_annot_set *structureAlignments = dataManager->GetStructureAlignments();
  698.     if (!structureAlignments) {
  699.         WARNINGMSG("StructureSet::AddStructureAlignment() - creating new structure alignment list");
  700.         InitStructureAlignments(objects.front()->mmdbID);
  701.         structureAlignments = dataManager->GetStructureAlignments();
  702.     }
  704.     // check master domain ID, to see if alignments have crossed master's domain boundaries
  705.     int *currentMasterDomainID = &(structureAlignments->SetFeatures().front().GetObject().SetId().Set());
  706.     if (*currentMasterDomainID == NO_DOMAIN)
  707.         *currentMasterDomainID = masterDomainID;
  708.     else if ((*currentMasterDomainID % 100) != (masterDomainID % 100))
  709.         *currentMasterDomainID = (*currentMasterDomainID / 100) * 100;
  711.     // check to see if this slave domain already has an alignment; if so, increment alignment #
  712.     CBiostruc_feature_set::TFeatures::const_iterator
  713.         f, fe = structureAlignments->GetFeatures().front().GetObject().GetFeatures().end();
  714.     for (f=structureAlignments->GetFeatures().front().GetObject().GetFeatures().begin(); f!=fe; ++f) {
  715.         if ((f->GetObject().GetId().Get() / 10) == (slaveDomainID / 10))
  716.             ++slaveDomainID;
  717.     }
  718.     CBiostruc_feature_id id(slaveDomainID);
  719.     feature->SetId(id);
  721.     CRef<CBiostruc_feature> featureRef(feature);
  722.     structureAlignments->SetFeatures().front().GetObject().SetFeatures().resize(
  723.         structureAlignments->GetFeatures().front().GetObject().GetFeatures().size() + 1, featureRef);
  725.     // flag a change in data
  726.     SetDataChanged(eStructureAlignmentData);
  727. }
  729. void StructureSet::RemoveStructureAlignments(void)
  730. {
  731.     dataManager->SetStructureAlignments(NULL);
  732.     // flag a change in data
  733.     SetDataChanged(eStructureAlignmentData);
  734. }
  736. void StructureSet::ReplaceAlignmentSet(AlignmentSet *newAlignmentSet)
  737. {
  738.     ASNDataManager::SeqAnnotList *seqAnnots = dataManager->GetOrCreateSequenceAlignments();
  739.     if (!seqAnnots) {
  740.         ERRORMSG("StructureSet::ReplaceAlignmentSet() - "
  741.             << "can't figure out where in the asn the alignments are to go");
  742.         return;
  743.     }
  745.     // update the AlignmentSet
  746.     if (alignmentSet) {
  747.         _RemoveChild(alignmentSet);
  748.         delete alignmentSet;
  749.     }
  750.     alignmentSet = newAlignmentSet;
  752.     // update the asn alignments
  753.     seqAnnots->resize(alignmentSet->newAsnAlignmentData->size());
  754.     ASNDataManager::SeqAnnotList::iterator o = seqAnnots->begin();
  755.     ASNDataManager::SeqAnnotList::iterator n, ne = alignmentSet->newAsnAlignmentData->end();
  756.     for (n=alignmentSet->newAsnAlignmentData->begin(); n!=ne; ++n, ++o)
  757.         o->Reset(n->GetPointer());   // copy each Seq-annot CRef
  759.     // don't set data PSSM/row order flags here; done by AlignmentManager::SavePairwiseFromMultiple()
  760.     SetDataChanged(eAnyAlignmentData);
  761. }
  763. void StructureSet::ReplaceUpdates(ncbi::objects::CCdd::TPending& newUpdates)
  764. {
  765.     dataManager->ReplaceUpdates(newUpdates);
  766. }
  768. void StructureSet::RemoveUnusedSequences(void)
  769. {
  770. ASNDataManager::SequenceList updateSequences;
  771.     if (alignmentManager) alignmentManager->GetUpdateSequences(&updateSequences);
  772.     dataManager->RemoveUnusedSequences(alignmentSet, updateSequences);
  773. }
  775. bool StructureSet::SaveASNData(const char *filename, bool doBinary, unsigned int *changeFlags)
  776. {
  777.     // force a save of any edits to alignment and updates first (it's okay if this has already been done)
  778.     GlobalMessenger()->SequenceWindowsSave(true);
  779.     /*if (dataManager->HasDataChanged())*/ RemoveUnusedSequences();
  781.     // create and temporarily attach a style dictionary, and annotation set + camera info
  782.     // to the data (and then remove it again, so it's never out of date)
  783.     CRef < CCn3d_style_dictionary > styleDictionary(styleManager->CreateASNStyleDictionary());
  784.     dataManager->SetStyleDictionary(*styleDictionary);
  785.     CRef < CCn3d_user_annotations > userAnnotations(new CCn3d_user_annotations());
  786.     if (!styleManager->SaveToASNUserAnnotations(userAnnotations.GetPointer()) ||
  787.         (objects.size() >= 1 && !renderer->SaveToASNViewSettings(userAnnotations.GetPointer()))) {
  788.         ERRORMSG("StructureSet::SaveASNData() - error creating user annotations blob");
  789.         return false;
  790.     }
  791.     if (userAnnotations->IsSetAnnotations() || userAnnotations->IsSetView())
  792.         dataManager->SetUserAnnotations(*userAnnotations);
  794.     string err;
  795.     bool writeOK = dataManager->WriteDataToFile(filename, doBinary, &err, eFNP_Replace);
  797.     // remove style dictionary and annotations from asn
  798.     dataManager->RemoveStyleDictionary();
  799.     dataManager->RemoveUserAnnotations();
  801.     if (writeOK) {
  802.         *changeFlags = dataManager->GetDataChanged();
  803.         dataManager->SetDataUnchanged();
  804.     } else {
  805.         ERRORMSG("Write failed: " << err);
  806.     }
  807.     return writeOK;
  808. }
  810. // because the frame map (for each frame, a list of diplay lists) is complicated
  811. // to create, this just verifies that all display lists occur exactly once
  812. // in the map. Also, make sure that total # display lists in all frames adds up.
  813. void StructureSet::VerifyFrameMap(void) const
  814. {
  815.     for (unsigned int l=OpenGLRenderer::FIRST_LIST; l<=lastDisplayList; ++l) {
  816.         bool found = false;
  817.         for (unsigned int f=0; f<frameMap.size(); ++f) {
  818.             DisplayLists::const_iterator d, de=frameMap[f].end();
  819.             for (d=frameMap[f].begin(); d!=de; ++d) {
  820.                 if (*d == l) {
  821.                     if (!found)
  822.                         found = true;
  823.                     else
  824.                         ERRORMSG("frameMap: repeated display list " << l);
  825.                 }
  826.             }
  827.         }
  828.         if (!found)
  829.             ERRORMSG("display list " << l << " not in frameMap");
  830.     }
  832.     unsigned int nLists = 0;
  833.     for (unsigned int f=0; f<frameMap.size(); ++f) {
  834.         DisplayLists::const_iterator d, de=frameMap[f].end();
  835.         for (d=frameMap[f].begin(); d!=de; ++d) ++nLists;
  836.     }
  837.     if (nLists != lastDisplayList)
  838.         ERRORMSG("frameMap has too many display lists");
  839. }
  841. void StructureSet::SetCenter(const Vector *given)
  842. {
  843.     Vector siteSum;
  844.     int nAtoms = 0;
  845.     double dist;
  846.     maxDistFromCenter = 0.0;
  848.     // set new center if given one
  849.     if (given) center = *given;
  851.     // loop trough all atoms twice - once to get average center, then once to
  852.     // find max distance from this center
  853.     for (int i=0; i<2; ++i) {
  854.         if (given && i==0) continue; // skip center calculation if given one
  855.         ObjectList::const_iterator o, oe=objects.end();
  856.         for (o=objects.begin(); o!=oe; ++o) {
  857.             StructureObject::CoordSetList::const_iterator c, ce=(*o)->coordSets.end();
  858.             for (c=(*o)->coordSets.begin(); c!=ce; ++c) {
  859.                 AtomSet::AtomMap::const_iterator a, ae=(*c)->atomSet->atomMap.end();
  860.                 for (a=(*c)->atomSet->atomMap.begin(); a!=ae; ++a) {
  861.                     Vector site(a->second.front()->site);
  862.                     if ((*o)->IsSlave() && (*o)->transformToMaster)
  863.                         ApplyTransformation(&site, *((*o)->transformToMaster));
  864.                     if (i==0) {
  865.                         siteSum += site;
  866.                         ++nAtoms;
  867.                     } else {
  868.                         dist = (site - center).length();
  869.                         if (dist > maxDistFromCenter)
  870.                             maxDistFromCenter = dist;
  871.                     }
  872.                 }
  873.             }
  874.         }
  875.         if (i==0) {
  876.             if (nAtoms == 0) {
  877.                 center.Set(0.0, 0.0, 0.0);
  878.                 break;
  879.             }
  880.             center = siteSum / nAtoms;
  881.         }
  882.     }
  883.     TRACEMSG("center: " << center << ", maxDistFromCenter " << maxDistFromCenter);
  884.     rotationCenter = center;
  885. }
  887. void StructureSet::CenterViewOnAlignedResidues(void)
  888. {
  889.     const BlockMultipleAlignment *alignment = alignmentManager->GetCurrentMultipleAlignment();
  890.     if (!alignment || !alignment->GetSequenceOfRow(0) || !alignment->GetSequenceOfRow(0))
  891.         return;                     // no alignment
  892.     const Molecule *masterMolecule = alignment->GetSequenceOfRow(0)->molecule;
  893.     if (!masterMolecule) return;    // no structured master
  894.     const StructureObject *masterObject;
  895.     if (!masterMolecule->GetParentOfType(&masterObject)) return;
  897.     // get coords of all aligned c-alphas
  898.     deque < Vector > coords;
  899.     Molecule::ResidueMap::const_iterator r, re = masterMolecule->residues.end();
  900.     for (r=masterMolecule->residues.begin(); r!=re; ++r) {
  901.         if (!alignment->IsAligned(0, r->first - 1)) continue;
  902.         if (r->second->alphaID == Residue::NO_ALPHA_ID) continue;
  903.         AtomPntr ap(masterMolecule->id, r->first, r->second->alphaID);
  904.         const AtomCoord* atom = masterObject->coordSets.front()->atomSet->GetAtom(ap, true, true);
  905.         if (atom) coords.push_back(atom->site);
  906.     }
  908.     // calculate center
  909.     int i;
  910.     Vector alignedCenter;
  911.     for (i=0; i<coords.size(); ++i) alignedCenter += coords[i];
  912.     alignedCenter /= coords.size();
  914.     // find radius
  915.     double radius = 0.0, d;
  916.     for (i=0; i<coords.size(); ++i) {
  917.         d = (coords[i] - alignedCenter).length();
  918.         if (d > radius) radius = d;
  919.     }
  921.     // set view
  922.     renderer->CenterView(alignedCenter, radius);
  923. }
  925. bool StructureSet::Draw(const AtomSet *atomSet) const
  926. {
  927.     TRACEMSG("drawing StructureSet");
  928.     if (!styleManager->CheckGlobalStyleSettings()) return false;
  929.     return true;
  930. }
  932. unsigned int StructureSet::CreateName(const Residue *residue, int atomID)
  933. {
  934.     ++lastAtomName;
  935.     nameMap[lastAtomName] = make_pair(residue, atomID);
  936.     return lastAtomName;
  937. }
  939. bool StructureSet::GetAtomFromName(unsigned int name, const Residue **residue, int *atomID) const
  940. {
  941.     NameMap::const_iterator i = nameMap.find(name);
  942.     if (i == nameMap.end()) return false;
  943.     *residue = i->second.first;
  944.     *atomID = i->second.second;
  945. return true;
  946. }
  948. void StructureSet::SelectedAtom(unsigned int name, bool setCenter)
  949. {
  950.     const Residue *residue;
  951.     int atomID;
  953.     if (name == OpenGLRenderer::NO_NAME || !GetAtomFromName(name, &residue, &atomID)) {
  954.         INFOMSG("nothing selected");
  955.         return;
  956.     }
  958.     // add highlight
  959.     const Molecule *molecule;
  960.     if (!residue->GetParentOfType(&molecule)) return;
  961.     GlobalMessenger()->ToggleHighlight(molecule, residue->id, true);
  962.     wxString molresid;
  963.     if (molecule->IsHeterogen() || molecule->IsSolvent()) {
  964.         const StructureObject *object;
  965.         if (molecule->GetParentOfType(&object)) {
  966.             // assume hets/solvents are single residue
  967.             if (object->pdbID.size() > 0)
  968.                 molresid.Printf("%s heterogen/solvent molecule %i", object->pdbID.c_str(), molecule->id);
  969.             else
  970.                 molresid = molecule->identifier->ToString().c_str();
  971.         }
  972.     } else
  973.         molresid.Printf("chain %s residue %i", molecule->identifier->ToString().c_str(), residue->id);
  974.     INFOMSG("selected " << molresid.c_str() << " (PDB: " << residue->nameGraph << ' ' << residue->namePDB
  975.         << ") atom " << atomID << " (PDB: " << residue->GetAtomInfo(atomID)->name << ')');
  977.     // get coordinate of picked atom, in coordinates of master frame
  978.     const StructureObject *object;
  979.     if (!molecule->GetParentOfType(&object)) return;
  980.     object->coordSets.front()->atomSet->SetActiveEnsemble(NULL);    // don't actually know which alternate...
  981.     Vector pickedAtomCoord = object->coordSets.front()->atomSet->
  982.         GetAtom(AtomPntr(molecule->id, residue->id, atomID)) ->site;
  983.     if (object->IsSlave() && object->transformToMaster)
  984.         ApplyTransformation(&pickedAtomCoord, *(object->transformToMaster));
  986.     // print out distance to previous picked atom
  987.     if (havePrevPickedAtomCoord)
  988.         INFOMSG("distance to previously selected atom: " << setprecision(3) <<
  989.             (pickedAtomCoord - prevPickedAtomCoord).length() << setprecision(6) << " A");
  990.     prevPickedAtomCoord = pickedAtomCoord;
  991.     havePrevPickedAtomCoord = true;
  993.     // if indicated, use atom site as rotation center; use coordinate from first CoordSet, default altConf
  994.     if (setCenter) {
  995.         INFOMSG("rotating about " << object->pdbID
  996.             << " molecule " << molecule->id << " residue " << residue->id << ", atom " << atomID);
  997.         rotationCenter = pickedAtomCoord;
  998.     }
  999. }
  1001. void StructureSet::SelectByDistance(double cutoff, bool biopolymersOnly, bool otherMoleculesOnly) const
  1002. {
  1003.     StructureObject::ResidueMap residuesToHighlight;
  1005.     // add residues to highlight to master list, based on proximities within objects
  1006.     ObjectList::const_iterator o, oe = objects.end();
  1007.     for (o=objects.begin(); o!=oe; ++o)
  1008.         (*o)->SelectByDistance(cutoff, biopolymersOnly, otherMoleculesOnly, &residuesToHighlight);
  1010.     // now actually add highlights for new selected residues
  1011.     StructureObject::ResidueMap::const_iterator r, re = residuesToHighlight.end();
  1012.     for (r=residuesToHighlight.begin(); r!=re; ++r)
  1013.         if (!GlobalMessenger()->IsHighlighted(r->second, r->first->id))
  1014.             GlobalMessenger()->ToggleHighlight(r->second, r->first->id, false);
  1015. }
  1017. const Sequence * StructureSet::CreateNewSequence(ncbi::objects::CBioseq& bioseq)
  1018. {
  1019.     // add Sequence to SequenceSet
  1020.     SequenceSet *modifiableSet = const_cast<SequenceSet*>(sequenceSet);
  1021.     const Sequence *newSeq = new Sequence(modifiableSet, bioseq);
  1022.     if (!newSeq->identifier) {
  1023.         ERRORMSG("StructureSet::CreateNewSequence() - identifier conflict, no new sequence created");
  1024.         delete newSeq;
  1025.         return NULL;
  1026.     }
  1027.     modifiableSet->sequences.push_back(newSeq);
  1029.     // add asn sequence to asn data
  1030.     if (dataManager->GetSequences()) {
  1031.         CSeq_entry *se = new CSeq_entry();
  1032.         se->SetSeq(bioseq);
  1033.         dataManager->GetSequences()->push_back(CRef<CSeq_entry>(se));
  1034.     } else
  1035.         ERRORMSG("StructureSet::CreateNewSequence() - no sequence list in asn data");
  1036.     SetDataChanged(eSequenceData);
  1038.     return newSeq;
  1039. }
  1041. void StructureSet::RejectAndPurgeSequence(const Sequence *reject, string reason, bool purge)
  1042. {
  1043.     if (!dataManager->IsCDD() || !reject || reason.size() == 0) return;
  1045.     CReject_id *rejectID = new CReject_id();
  1046.     rejectID->SetIds() = reject->bioseqASN->GetId();    // copy Seq-id lists
  1048.     CUpdate_comment *comment = new CUpdate_comment();
  1049.     comment->SetComment(reason);
  1050.     rejectID->SetDescription().push_back(CRef < CUpdate_comment > (comment));
  1052.     dataManager->AddReject(rejectID);
  1054.     if (purge)
  1055.         alignmentManager->PurgeSequence(reject->identifier);
  1056. }
  1058. const StructureSet::RejectList * StructureSet::GetRejects(void) const
  1059. {
  1060.     return dataManager->GetRejects();
  1061. }
  1063. void StructureSet::ShowRejects(void) const
  1064. {
  1065.     const RejectList *rejects = GetRejects();
  1066.     if (!rejects) {
  1067.         INFOMSG("No rejects in this CD");
  1068.         return;
  1069.     }
  1071.     INFOMSG("Rejects:");
  1072.     RejectList::const_iterator r, re = rejects->end();
  1073.     for (r=rejects->begin(); r!=re; ++r) {
  1074.         string idstr;
  1075.         CReject_id::TIds::const_iterator i, ie = (*r)->GetIds().end();
  1076.         for (i=(*r)->GetIds().begin(); i!=ie; ++i)
  1077.             idstr += (*i)->AsFastaString() + ", ";
  1078.         INFOMSG(idstr << "Reason: " <<
  1079.             (((*r)->IsSetDescription() && (*r)->GetDescription().front()->IsComment()) ?
  1080.                 (*r)->GetDescription().front()->GetComment() : string("none given")));
  1081.     }
  1082. }
  1084. bool StructureSet::ConvertMimeDataToCDD(const std::string& cddName)
  1085. {
  1086.     if (!dataManager->ConvertMimeDataToCDD(cddName))
  1087.         return false;
  1089.     // make sure all structured sequences have MMDB annot tags
  1090.     SequenceSet::SequenceList::const_iterator s, se = sequenceSet->sequences.end();
  1091.     for (s=sequenceSet->sequences.begin(); s!=se; ++s) {
  1092.         if ((*s)->molecule) {
  1093.             if ((*s)->identifier->mmdbID == MoleculeIdentifier::VALUE_NOT_SET) {
  1094.                 ERRORMSG("sequence " << (*s)->identifier->ToString()
  1095.                     << " has associated molecule but no MMDB id");
  1096.                 return false;
  1097.             }
  1098.             (*s)->AddMMDBAnnotTag((*s)->identifier->mmdbID);
  1099.         }
  1100.     }
  1102.     return true;
  1103. }
  1105. // trivial methods...
  1106. bool StructureSet::IsMultiStructure(void) const { return !dataManager->IsSingleStructure(); }
  1107. bool StructureSet::HasDataChanged(void) const { return dataManager->HasDataChanged(); }
  1108. void StructureSet::SetDataChanged(unsigned int what) const { dataManager->SetDataChanged(what); }
  1109. bool StructureSet::IsCDD(void) const { return dataManager->IsCDD(); }
  1110. bool StructureSet::IsCDDInMime(void) const { return dataManager->IsCDDInMime(); }
  1111. const string& StructureSet::GetCDDName(void) const { return dataManager->GetCDDName(); }
  1112. bool StructureSet::SetCDDName(const string& name) { return dataManager->SetCDDName(name); }
  1113. const string& StructureSet::GetCDDDescription(void) const { return dataManager->GetCDDDescription(); }
  1114. bool StructureSet::SetCDDDescription(const string& descr) { return dataManager->SetCDDDescription(descr); }
  1115. bool StructureSet::GetCDDNotes(StructureSet::TextLines *lines) const { return dataManager->GetCDDNotes(lines); }
  1116. bool StructureSet::SetCDDNotes(const StructureSet::TextLines& lines) { return dataManager->SetCDDNotes(lines); }
  1117. ncbi::objects::CCdd_descr_set * StructureSet::GetCDDDescrSet(void) { return dataManager->GetCDDDescrSet(); }
  1118. ncbi::objects::CAlign_annot_set * StructureSet::GetCDDAnnotSet(void) { return dataManager->GetCDDAnnotSet(); }
  1122. ///// StructureObject stuff /////
  1124. const int StructureObject::NO_MMDB_ID = -1;
  1125. const double StructureObject::NO_TEMPERATURE = kMin_Double;
  1127. StructureObject::StructureObject(StructureBase *parent, const CBiostruc& biostruc, bool master) :
  1128.     StructureBase(parent), isMaster(master), mmdbID(NO_MMDB_ID), transformToMaster(NULL),
  1129.     minTemperature(NO_TEMPERATURE), maxTemperature(NO_TEMPERATURE)
  1130. {
  1131.     // set numerical id simply based on # objects in parentSet
  1132.     id = parentSet->objects.size() + 1;
  1134.     // get MMDB id
  1135.     CBiostruc::TId::const_iterator j, je=biostruc.GetId().end();
  1136.     for (j=biostruc.GetId().begin(); j!=je; ++j) {
  1137.         if (j->GetObject().IsMmdb_id()) {
  1138.             mmdbID = j->GetObject().GetMmdb_id().Get();
  1139.             break;
  1140.         }
  1141.     }
  1142.     TRACEMSG("MMDB id " << mmdbID);
  1144.     // get PDB id
  1145.     if (biostruc.IsSetDescr()) {
  1146.         CBiostruc::TDescr::const_iterator k, ke=biostruc.GetDescr().end();
  1147.         for (k=biostruc.GetDescr().begin(); k!=ke; ++k) {
  1148.             if (k->GetObject().IsName()) {
  1149.                 pdbID = k->GetObject().GetName();
  1150.                 break;
  1151.             }
  1152.         }
  1153.     }
  1154.     TRACEMSG("PDB id " << pdbID);
  1156.     // get atom and feature spatial coordinates
  1157.     if (biostruc.IsSetModel()) {
  1158.         // iterate SEQUENCE OF Biostruc-model
  1159.         CBiostruc::TModel::const_iterator i, ie=biostruc.GetModel().end();
  1160.         for (i=biostruc.GetModel().begin(); i!=ie; ++i) {
  1162.             // don't know how to deal with these...
  1163.             if (i->GetObject().GetType() == eModel_type_ncbi_vector ||
  1164.                 i->GetObject().GetType() == eModel_type_other) continue;
  1166. //            // special case, typically for loading CDD's, when we're only interested in a single model type
  1167. //            if (isRawBiostrucFromMMDB && i->GetObject().GetType() != eModel_type_ncbi_all_atom) continue;
  1169.             // otherwise, assume all models in this set are of same type
  1170.             if (i->GetObject().GetType() == eModel_type_ncbi_backbone)
  1171.                 parentSet->isAlphaOnly = true;
  1172.             else
  1173.                 parentSet->isAlphaOnly = false;
  1175.             // load each Biostruc-model into a CoordSet
  1176.             if (i->GetObject().IsSetModel_coordinates()) {
  1177.                 CoordSet *coordSet =
  1178.                     new CoordSet(this, i->GetObject().GetModel_coordinates());
  1179.                 coordSets.push_back(coordSet);
  1180.             }
  1181.         }
  1182.     }
  1184.     TRACEMSG("temperature range: " << minTemperature << " to " << maxTemperature);
  1186.     // get graph - must be done after atom coordinates are loaded, so we can
  1187.     // avoid storing graph nodes for atoms not present in the model
  1188.     graph = new ChemicalGraph(this, biostruc.GetChemical_graph(), biostruc.GetFeatures());
  1189. }
  1191. bool StructureObject::SetTransformToMaster(const CBiostruc_annot_set& annot, int masterMMDBID)
  1192. {
  1193.     CBiostruc_annot_set::TFeatures::const_iterator f1, f1e=annot.GetFeatures().end();
  1194.     for (f1=annot.GetFeatures().begin(); f1!=f1e; ++f1) {
  1195.         CBiostruc_feature_set::TFeatures::const_iterator f2, f2e=f1->GetObject().GetFeatures().end();
  1196.         for (f2=f1->GetObject().GetFeatures().begin(); f2!=f2e; ++f2) {
  1198.             // skip if already used
  1199.             if (f2->GetObject().IsSetId() &&
  1200.                     parentSet->usedFeatures.find(f2->GetObject().GetId().Get()) !=
  1201.                         parentSet->usedFeatures.end())
  1202.                 continue;
  1204.             // look for alignment feature
  1205.             if (f2->GetObject().IsSetType() &&
  1206. f2->GetObject().GetType() == CBiostruc_feature::eType_alignment &&
  1207.                 f2->GetObject().IsSetLocation() &&
  1208.                 f2->GetObject().GetLocation().IsAlignment()) {
  1209.                 const CChem_graph_alignment& graphAlign =
  1210. f2->GetObject().GetLocation().GetAlignment();
  1212.                 // find transform alignment of this object with master
  1213.                 if (graphAlign.GetDimension() == 2 &&
  1214.                     graphAlign.GetBiostruc_ids().size() == 2 &&
  1215.                     graphAlign.IsSetTransform() &&
  1216.                     graphAlign.GetTransform().size() == 1 &&
  1217.                     graphAlign.GetBiostruc_ids().front().GetObject().IsMmdb_id() &&
  1218.                     graphAlign.GetBiostruc_ids().front().GetObject().GetMmdb_id().Get() == masterMMDBID &&
  1219.                     graphAlign.GetBiostruc_ids().back().GetObject().IsMmdb_id() &&
  1220.                     graphAlign.GetBiostruc_ids().back().GetObject().GetMmdb_id().Get() == mmdbID) {
  1222.                     // mark feature as used
  1223.                     if (f2->GetObject().IsSetId())
  1224.                         parentSet->usedFeatures[f2->GetObject().GetId().Get()] = true;
  1225.                     TRACEMSG("got transform for " << pdbID << "->master");
  1227.                     // unpack transform into matrix, moves in reverse order;
  1228.                     Matrix xform;
  1229.                     transformToMaster = new Matrix();
  1230.                     CTransform::TMoves::const_iterator
  1231.                         m, me=graphAlign.GetTransform().front().GetObject().GetMoves().end();
  1232.                     for (m=graphAlign.GetTransform().front().GetObject().GetMoves().begin(); m!=me; ++m) {
  1233.                         Matrix xmat;
  1234.                         double scale;
  1235.                         if (m->GetObject().IsTranslate()) {
  1236.                             const CTrans_matrix& trans = m->GetObject().GetTranslate();
  1237.                             scale = 1.0 / trans.GetScale_factor();
  1238.                             SetTranslationMatrix(&xmat,
  1239.                                 Vector(scale * trans.GetTran_1(),
  1240.                                        scale * trans.GetTran_2(),
  1241.                                        scale * trans.GetTran_3()));
  1242.                         } else { // rotate
  1243.                             const CRot_matrix& rot = m->GetObject().GetRotate();
  1244.                             scale = 1.0 / rot.GetScale_factor();
  1245.                             xmat.m[0]=scale*rot.GetRot_11(); xmat.m[1]=scale*rot.GetRot_21(); xmat.m[2]= scale*rot.GetRot_31(); xmat.m[3]=0;
  1246.                             xmat.m[4]=scale*rot.GetRot_12(); xmat.m[5]=scale*rot.GetRot_22(); xmat.m[6]= scale*rot.GetRot_32(); xmat.m[7]=0;
  1247.                             xmat.m[8]=scale*rot.GetRot_13(); xmat.m[9]=scale*rot.GetRot_23(); xmat.m[10]=scale*rot.GetRot_33(); xmat.m[11]=0;
  1248.                             xmat.m[12]=0;                    xmat.m[13]=0;                    xmat.m[14]=0;                     xmat.m[15]=1;
  1249.                         }
  1250.                         ComposeInto(transformToMaster, xmat, xform);
  1251.                         xform = *transformToMaster;
  1252.                     }
  1253.                     return true;
  1254.                 }
  1255.             }
  1256.         }
  1257.     }
  1259.     WARNINGMSG("Can't get structure alignment for slave " << pdbID
  1260.         << " with master " << masterMMDBID << ";nwill likely require manual realignment");
  1261.     return false;
  1262. }
  1264. static void AddDomain(int *domain, const Molecule *molecule, const Block::Range *range)
  1265. {
  1266.     const StructureObject *object;
  1267.     if (!molecule->GetParentOfType(&object)) return;
  1268.     for (int l=range->from; l<=range->to; ++l) {
  1269.         if (molecule->residueDomains[l] != Molecule::NO_DOMAIN_SET) {
  1270.             if (*domain == NO_DOMAIN)
  1271.                 *domain = object->domainID2MMDB.find(molecule->residueDomains[l])->second;
  1272.             else if (*domain != MULTI_DOMAIN &&
  1273.                      *domain != object->domainID2MMDB.find(molecule->residueDomains[l])->second)
  1274.                 *domain = MULTI_DOMAIN;
  1275.         }
  1276.     }
  1277. }
  1279. void StructureObject::RealignStructure(int nCoords,
  1280.     const Vector * const *masterCoords, const Vector * const *slaveCoords,
  1281.     const double *weights, int slaveRow)
  1282. {
  1283.     Vector masterCOM, slaveCOM; // centers of mass for master, slave
  1284.     Matrix slaveRotation;       // rotation to align slave with master
  1285.     if (!transformToMaster) transformToMaster = new Matrix();
  1287.     // do the fit
  1288.     RigidBodyFit(nCoords, masterCoords, slaveCoords, weights, masterCOM, slaveCOM, slaveRotation);
  1290.     // apply the resulting transform elements from the fit to this object's transform Matrix
  1291.     Matrix single, combined;
  1292.     SetTranslationMatrix(&single, -slaveCOM);
  1293.     ComposeInto(transformToMaster, slaveRotation, single);
  1294.     combined = *transformToMaster;
  1295.     SetTranslationMatrix(&single, masterCOM);
  1296.     ComposeInto(transformToMaster, single, combined);
  1298.     // print out RMSD
  1299.     Vector x;
  1300.     double rmsd = 0.0, d;
  1301.     int n = 0;
  1302.     for (int c=0; c<nCoords; ++c) {
  1303.         if (!slaveCoords[c] || !masterCoords[c]) continue;
  1304.         x = *(slaveCoords[c]);
  1305.         ApplyTransformation(&x, *transformToMaster);
  1306.         d = (*(masterCoords[c]) - x).length();
  1307.         rmsd += d * d;
  1308.         ++n;
  1309.     }
  1310.     rmsd = sqrt(rmsd / n);
  1311.     INFOMSG("RMSD of aligned alpha coordinates between master structure and " << pdbID << ": "
  1312.         << setprecision(3) << rmsd << setprecision(6) << " A");
  1314.     // create a new Biostruc-feature that contains this alignment
  1315.     CBiostruc_feature *feature = new CBiostruc_feature();
  1316.     feature->SetType(CBiostruc_feature::eType_alignment);
  1317.     CRef<CBiostruc_feature::C_Location> location(new CBiostruc_feature::C_Location());
  1318.     feature->SetLocation(*location);
  1319.     CChem_graph_alignment *graphAlignment = new CChem_graph_alignment();
  1320.     location.GetObject().SetAlignment(*graphAlignment);
  1322.     // fill out the Chem-graph-alignment
  1323.     graphAlignment->SetDimension(2);
  1324.     CMmdb_id
  1325.         *masterMID = new CMmdb_id(parentSet->objects.front()->mmdbID),
  1326.         *slaveMID = new CMmdb_id(mmdbID);
  1327.     CBiostruc_id
  1328.         *masterBID = new CBiostruc_id(),
  1329.         *slaveBID = new CBiostruc_id();
  1330.     masterBID->SetMmdb_id(*masterMID);
  1331.     slaveBID->SetMmdb_id(*slaveMID);
  1332.     graphAlignment->SetBiostruc_ids().resize(2);
  1333.     graphAlignment->SetBiostruc_ids().front().Reset(masterBID);
  1334.     graphAlignment->SetBiostruc_ids().back().Reset(slaveBID);
  1335.     graphAlignment->SetAlignment().resize(2);
  1337.     // fill out sequence alignment intervals, tracking domains in alignment
  1338.     const BlockMultipleAlignment *multiple = parentSet->alignmentManager->GetCurrentMultipleAlignment();
  1339.     int masterDomain = NO_DOMAIN, slaveDomain = NO_DOMAIN;
  1340.     const Molecule
  1341.         *masterMolecule = multiple->GetSequenceOfRow(0)->molecule,
  1342.         *slaveMolecule = multiple->GetSequenceOfRow(slaveRow)->molecule;
  1343.     BlockMultipleAlignment::UngappedAlignedBlockList blocks;
  1344.     multiple->GetUngappedAlignedBlocks(&blocks);
  1345.     if (blocks.size() > 0) {
  1346.         CChem_graph_pntrs
  1347.             *masterCGPs = new CChem_graph_pntrs(),
  1348.             *slaveCGPs = new CChem_graph_pntrs();
  1349.         graphAlignment->SetAlignment().front().Reset(masterCGPs);
  1350.         graphAlignment->SetAlignment().back().Reset(slaveCGPs);
  1351.         CResidue_pntrs
  1352.             *masterRPs = new CResidue_pntrs(),
  1353.             *slaveRPs = new CResidue_pntrs();
  1354.         masterCGPs->SetResidues(*masterRPs);
  1355.         slaveCGPs->SetResidues(*slaveRPs);
  1357.         masterRPs->SetInterval().resize(blocks.size());
  1358.         slaveRPs->SetInterval().resize(blocks.size());
  1359.         BlockMultipleAlignment::UngappedAlignedBlockList::const_iterator b, be = blocks.end();
  1360.         CResidue_pntrs::TInterval::iterator
  1361.             mi = masterRPs->SetInterval().begin(),
  1362.             si = slaveRPs->SetInterval().begin();
  1363.         for (b=blocks.begin(); b!=be; ++b, ++mi, ++si) {
  1364.             CResidue_interval_pntr
  1365.                 *masterRIP = new CResidue_interval_pntr(),
  1366.                 *slaveRIP = new CResidue_interval_pntr();
  1367.             mi->Reset(masterRIP);
  1368.             si->Reset(slaveRIP);
  1370.             masterRIP->SetMolecule_id().Set(masterMolecule->id);
  1371.             slaveRIP->SetMolecule_id().Set(slaveMolecule->id);
  1373.             const Block::Range *range = (*b)->GetRangeOfRow(0);
  1374.             masterRIP->SetFrom().Set(range->from + 1); // +1 to convert seqLoc to residueID
  1375.             masterRIP->SetTo().Set(range->to + 1);
  1376.             AddDomain(&masterDomain, masterMolecule, range);
  1378.             range = (*b)->GetRangeOfRow(slaveRow);
  1379.             slaveRIP->SetFrom().Set(range->from + 1);
  1380.             slaveRIP->SetTo().Set(range->to + 1);
  1381.             AddDomain(&slaveDomain, slaveMolecule, range);
  1382.         }
  1383.     }
  1385.     // fill out structure alignment transform
  1386.     CTransform *xform = new CTransform();
  1387.     graphAlignment->SetTransform().resize(1);
  1388.     graphAlignment->SetTransform().front().Reset(xform);
  1389.     xform->SetId(1);
  1390.     xform->SetMoves().resize(3);
  1391.     CTransform::TMoves::iterator m = xform->SetMoves().begin();
  1392.     for (int i=0; i<3; ++i, ++m) {
  1393.         CMove *move = new CMove();
  1394.         m->Reset(move);
  1395.         static const int scaleFactor = 100000;
  1396.         if (i == 0) {   // translate slave so its COM is at origin
  1397.             CTrans_matrix *trans = new CTrans_matrix();
  1398.             move->SetTranslate(*trans);
  1399.             trans->SetScale_factor(scaleFactor);
  1400.             trans->SetTran_1((int)(-(slaveCOM.x * scaleFactor)));
  1401.             trans->SetTran_2((int)(-(slaveCOM.y * scaleFactor)));
  1402.             trans->SetTran_3((int)(-(slaveCOM.z * scaleFactor)));
  1403.         } else if (i == 1) {
  1404.             CRot_matrix *rot = new CRot_matrix();
  1405.             move->SetRotate(*rot);
  1406.             rot->SetScale_factor(scaleFactor);
  1407.             rot->SetRot_11((int)(slaveRotation[0] * scaleFactor));
  1408.             rot->SetRot_12((int)(slaveRotation[4] * scaleFactor));
  1409.             rot->SetRot_13((int)(slaveRotation[8] * scaleFactor));
  1410.             rot->SetRot_21((int)(slaveRotation[1] * scaleFactor));
  1411.             rot->SetRot_22((int)(slaveRotation[5] * scaleFactor));
  1412.             rot->SetRot_23((int)(slaveRotation[9] * scaleFactor));
  1413.             rot->SetRot_31((int)(slaveRotation[2] * scaleFactor));
  1414.             rot->SetRot_32((int)(slaveRotation[6] * scaleFactor));
  1415.             rot->SetRot_33((int)(slaveRotation[10] * scaleFactor));
  1416.         } else if (i == 2) {    // translate slave so its COM is at COM of master
  1417.             CTrans_matrix *trans = new CTrans_matrix();
  1418.             move->SetTranslate(*trans);
  1419.             trans->SetScale_factor(scaleFactor);
  1420.             trans->SetTran_1((int)(masterCOM.x * scaleFactor));
  1421.             trans->SetTran_2((int)(masterCOM.y * scaleFactor));
  1422.             trans->SetTran_3((int)(masterCOM.z * scaleFactor));
  1423.         }
  1424.     }
  1426.     // store the new alignment in the Biostruc-annot-set,
  1427.     // setting the feature id depending on the aligned domain(s)
  1428.     if (masterDomain == NO_DOMAIN) masterDomain = 0;    // can happen if single domain chain
  1429.     if (slaveDomain == NO_DOMAIN) slaveDomain = 0;
  1430.     const StructureObject *masterObject;
  1431.     if (!masterMolecule->GetParentOfType(&masterObject)) return;
  1432.     int
  1433.         masterDomainID = masterObject->mmdbID*10000 + masterMolecule->id*100 + masterDomain,
  1434.         slaveDomainID = mmdbID*100000 + slaveMolecule->id*1000 + slaveDomain*10 + 1;
  1435.     parentSet->AddStructureAlignment(feature, masterDomainID, slaveDomainID);
  1437.     // for backward-compatibility with Cn3D 3.5, need name to encode chain/domain
  1438.     CNcbiOstrstream oss;
  1439.     oss << masterMolecule->identifier->pdbID << ((char) masterMolecule->identifier->pdbChain) << masterDomain << ' '
  1440.         << slaveMolecule->identifier->pdbID << ((char) slaveMolecule->identifier->pdbChain) << slaveDomain << ' '
  1441.         << "Structure alignment of slave " << multiple->GetSequenceOfRow(slaveRow)->identifier->ToString()
  1442.         << " with master " << multiple->GetSequenceOfRow(0)->identifier->ToString()
  1443.         << ", as computed by Cn3D" << '';
  1444.     feature->SetName(string(oss.str()));
  1445.     delete oss.str();
  1446. }
  1448. void StructureObject::SelectByDistance(double cutoff, bool biopolymersOnly, bool otherMoleculesOnly,
  1449.     ResidueMap *selectedResidues) const
  1450. {
  1451.     // first make a list of coordinates of atoms in selected residues,
  1452.     typedef vector < const AtomCoord * > CoordList;
  1453.     CoordList highlightedAtoms;
  1454.     typedef map < const Molecule * , bool > MoleculeList;
  1455.     MoleculeList moleculesWithHighlights;
  1457.     ChemicalGraph::MoleculeMap::const_iterator m, me = graph->molecules.end();
  1458.     for (m=graph->molecules.begin(); m!=me; ++m) {
  1459.         Molecule::ResidueMap::const_iterator r, re = m->second->residues.end();
  1460.         for (r=m->second->residues.begin(); r!=re; ++r) {
  1462.             if (GlobalMessenger()->IsHighlighted(m->second, r->second->id)) {
  1463.                 const Residue::AtomInfoMap& atomInfos = r->second->GetAtomInfos();
  1464.                 Residue::AtomInfoMap::const_iterator a, ae = atomInfos.end();
  1465.                 for (a=atomInfos.begin(); a!=ae; ++a) {
  1466.                     const AtomCoord *atomCoord = coordSets.front()->atomSet->
  1467.                         GetAtom(AtomPntr(m->second->id, r->second->id, a->first), true, true);
  1468.                     if (atomCoord) highlightedAtoms.push_back(atomCoord);
  1469.                 }
  1470.                 moleculesWithHighlights[m->second] = true;
  1471.             }
  1472.         }
  1473.     }
  1474.     if (highlightedAtoms.size() == 0) return;
  1476.     // now make a list of unselected residues to check for proximity
  1477.     typedef vector < const Residue * > ResidueList;
  1478.     ResidueList unhighlightedResidues;
  1480.     for (m=graph->molecules.begin(); m!=me; ++m) {
  1481.         Molecule::ResidueMap::const_iterator r, re = m->second->residues.end();
  1483.         if (otherMoleculesOnly &&
  1484.                 moleculesWithHighlights.find(m->second) != moleculesWithHighlights.end())
  1485.             continue;
  1487.         for (r=m->second->residues.begin(); r!=re; ++r) {
  1489.             if (!GlobalMessenger()->IsHighlighted(m->second, r->second->id) &&
  1490.                     (!biopolymersOnly || m->second->IsProtein() || m->second->IsNucleotide()))
  1491.                 unhighlightedResidues.push_back(r->second);
  1492.         }
  1493.     }
  1494.     if (unhighlightedResidues.size() == 0) return;
  1496.     // now check all unhighlighted residues, to see if any atoms are within cutoff distance
  1497.     // of any highlighted atoms; if so, add to residue selection list
  1498.     CoordList::const_iterator h, he = highlightedAtoms.end();
  1499.     ResidueList::const_iterator u, ue = unhighlightedResidues.end();
  1500.     for (u=unhighlightedResidues.begin(); u!=ue; ++u) {
  1501.         const Molecule *molecule;
  1502.         if (!(*u)->GetParentOfType(&molecule)) continue;
  1504.         const Residue::AtomInfoMap& atomInfos = (*u)->GetAtomInfos();
  1505.         Residue::AtomInfoMap::const_iterator a, ae = atomInfos.end();
  1506.         for (a=atomInfos.begin(); a!=ae; ++a) {
  1507.             const AtomCoord *uAtomCoord = coordSets.front()->atomSet->
  1508.                 GetAtom(AtomPntr(molecule->id, (*u)->id, a->first), true, true);
  1509.             if (!uAtomCoord) continue;
  1511.             // for each unhighlighted atom, check for proximity to a highlighted atom
  1512.             for (h=highlightedAtoms.begin(); h!=he; ++h) {
  1513.                 if ((uAtomCoord->site - (*h)->site).length() <= cutoff)
  1514.                     break;
  1515.             }
  1516.             if (h != he)
  1517.                 break;
  1518.         }
  1520.         // if any atom of this residue is near a highlighted atom, add to selection list
  1521.         if (a != ae)
  1522.             (*selectedResidues)[*u] = molecule;
  1523.     }
  1524. }
  1526. END_SCOPE(Cn3D)
  1528. /*
  1529. * ---------------------------------------------------------------------------
  1530. * $Log: structure_set.cpp,v $
  1531. * Revision 1000.2  2004/06/01 18:29:28  gouriano
  1532. * PRODUCTION: UPGRADED [GCC34_MSVC7] Dev-tree R1.141
  1533. *
  1534. * Revision 1.141  2004/05/21 21:41:40  gorelenk
  1535. * Added PCH ncbi_pch.hpp
  1536. *
  1537. * Revision 1.140  2004/05/21 17:29:51  thiessen
  1538. * allow conversion of mime to cdd data
  1539. *
  1540. * Revision 1.139  2004/03/15 18:32:04  thiessen
  1541. * prefer prefix vs. postfix ++/-- operators
  1542. *
  1543. * Revision 1.138  2004/02/19 17:05:14  thiessen
  1544. * remove cn3d/ from include paths; add pragma to disable annoying msvc warning
  1545. *
  1546. * Revision 1.137  2004/02/05 19:00:08  thiessen
  1547. * change mmdb->pdb upon het selection
  1548. *
  1549. * Revision 1.136  2004/01/05 17:09:16  thiessen
  1550. * abort import and warn if same accession different gi
  1551. *
  1552. * Revision 1.135  2003/10/21 13:48:48  grichenk
  1553. * Redesigned type aliases in serialization library.
  1554. * Fixed the code (removed CRef-s, added explicit
  1555. * initializers etc.)
  1556. *
  1557. * Revision 1.134  2003/07/17 20:13:51  thiessen
  1558. * don't ask for filename on save/termination with -f; add eAnyAlignmentData flag
  1559. *
  1560. * Revision 1.133  2003/07/17 16:52:34  thiessen
  1561. * add FileSaved message with edit typing
  1562. *
  1563. * Revision 1.132  2003/07/14 18:37:08  thiessen
  1564. * change GetUngappedAlignedBlocks() param types; other syntax changes
  1565. *
  1566. * Revision 1.131  2003/06/12 14:38:46  thiessen
  1567. * fix empty feature list bug in blank structure alignment data
  1568. *
  1569. * Revision 1.130  2003/04/02 17:49:18  thiessen
  1570. * allow pdb id's in structure import dialog
  1571. *
  1572. * Revision 1.129  2003/02/05 14:55:22  thiessen
  1573. * always load single structure even if structureLimit == 0
  1574. *
  1575. * Revision 1.128  2003/02/03 19:20:07  thiessen
  1576. * format changes: move CVS Log to bottom of file, remove std:: from .cpp files, and use new diagnostic macros
  1577. *
  1578. * Revision 1.127  2003/01/27 15:52:22  thiessen
  1579. * merge after highlighted row; show rejects; trim rejects from new reject list
  1580. *
  1581. * Revision 1.126  2003/01/10 20:10:53  thiessen
  1582. * undo previous change - don't want full Bioseqs (with multiple Seq-ids) in converted C-Bioseqs
  1583. *
  1584. * Revision 1.125  2003/01/09 16:15:42  thiessen
  1585. * use ConvertAsnFromCPPToC to convert full Bioseqs
  1586. *
  1587. * Revision 1.124  2002/12/19 18:52:28  thiessen
  1588. * rms -> rmsd
  1589. *
  1590. * Revision 1.123  2002/11/21 17:48:12  thiessen
  1591. * fix yet another user style bug
  1592. *
  1593. * Revision 1.122  2002/11/21 15:26:33  thiessen
  1594. * fix style dictionary loading bug
  1595. *
  1596. * Revision 1.121  2002/11/19 21:19:44  thiessen
  1597. * more const changes for objects; fix user vs default style bug
  1598. *
  1599. * Revision 1.120  2002/11/18 19:48:40  grichenk
  1600. * Removed "const" from datatool-generated setters
  1601. *
  1602. * Revision 1.119  2002/11/18 15:02:40  thiessen
  1603. * set flags on style shortcut menus
  1604. *
  1605. * Revision 1.118  2002/11/06 00:18:10  thiessen
  1606. * fixes for new CRef/const rules in objects
  1607. *
  1608. * Revision 1.117  2002/10/08 12:35:43  thiessen
  1609. * use delete[] for arrays
  1610. *
  1611. * Revision 1.116  2002/09/30 17:13:02  thiessen
  1612. * change structure import to do sequences as well; change cache to hold mimes; change block aligner vocabulary; fix block aligner dialog bugs
  1613. *
  1614. * Revision 1.115  2002/09/26 18:30:58  thiessen
  1615. * fix RMS bug when coords missing
  1616. *
  1617. * Revision 1.114  2002/09/26 17:32:13  thiessen
  1618. * show distance between picked atoms; show RMS for structure alignments
  1619. *
  1620. * Revision 1.113  2002/09/19 12:51:08  thiessen
  1621. * fix block aligner / update bug; add distance select for other molecules only
  1622. *
  1623. * Revision 1.112  2002/09/16 21:24:58  thiessen
  1624. * add block freezing to block aligner
  1625. *
  1626. * Revision 1.111  2002/09/14 19:18:32  thiessen
  1627. * fix center-on-aligned bug when no structured master
  1628. *
  1629. * Revision 1.110  2002/09/14 18:41:05  thiessen
  1630. * fix center-on-aligned omission
  1631. *
  1632. * Revision 1.109  2002/09/14 17:03:07  thiessen
  1633. * center initial view on aligned residues
  1634. *
  1635. * Revision 1.108  2002/09/09 13:38:23  thiessen
  1636. * separate save and save-as
  1637. *
  1638. * Revision 1.107  2002/08/15 22:13:17  thiessen
  1639. * update for wx2.3.2+ only; add structure pick dialog; fix MultitextDialog bug
  1640. *
  1641. * Revision 1.106  2002/07/12 13:24:10  thiessen
  1642. * fixes for PSSM creation to agree with cddumper/RPSBLAST
  1643. *
  1644. * Revision 1.105  2002/07/03 13:39:40  thiessen
  1645. * fix for redundant sequence removal
  1646. *
  1647. * Revision 1.104  2002/07/01 15:30:21  thiessen
  1648. * fix for container type switch in Dense-seg
  1649. *
  1650. * Revision 1.103  2002/06/07 12:41:04  thiessen
  1651. * change ambiguous master error to warning
  1652. *
  1653. * Revision 1.102  2002/06/05 17:50:08  thiessen
  1654. * title tweaks
  1655. *
  1656. * Revision 1.101  2002/04/27 16:32:14  thiessen
  1657. * fix small leaks/bugs found by BoundsChecker
  1658. *
  1659. * Revision 1.100  2002/04/11 16:39:54  thiessen
  1660. * fix style manager bug
  1661. *
  1662. * Revision 1.99  2002/04/10 13:16:28  thiessen
  1663. * new selection by distance algorithm
  1664. *
  1665. * Revision 1.98  2002/03/28 14:06:02  thiessen
  1666. * preliminary BLAST/PSSM ; new CD startup style
  1667. *
  1668. * Revision 1.97  2002/02/27 16:29:41  thiessen
  1669. * add model type flag to general mime type
  1670. *
  1671. * Revision 1.96  2002/02/19 14:59:40  thiessen
  1672. * add CDD reject and purge sequence
  1673. *
  1674. * Revision 1.95  2002/02/12 17:19:22  thiessen
  1675. * first working structure import
  1676. *
  1677. * Revision 1.94  2002/02/05 18:53:26  thiessen
  1678. * scroll to residue in sequence windows when selected in structure window
  1679. *
  1680. * Revision 1.93  2002/01/19 02:34:50  thiessen
  1681. * fixes for changes in asn serialization API
  1682. *
  1683. * Revision 1.92  2002/01/03 16:18:40  thiessen
  1684. * add distance selection
  1685. *
  1686. * Revision 1.91  2001/12/21 23:08:48  thiessen
  1687. * add residue info when structure selected
  1688. *
  1689. * Revision 1.90  2001/12/15 03:15:59  thiessen
  1690. * adjustments for slightly changed object loader Set...() API
  1691. *
  1692. * Revision 1.89  2001/12/12 14:04:15  thiessen
  1693. * add missing object headers after object loader change
  1694. *
  1695. * Revision 1.88  2001/12/06 23:13:46  thiessen
  1696. * finish import/align new sequences into single-structure data; many small tweaks
  1697. *
  1698. * Revision 1.87  2001/11/30 14:02:05  thiessen
  1699. * progress on sequence imports to single structures
  1700. *
  1701. * Revision 1.86  2001/11/27 16:26:09  thiessen
  1702. * major update to data management system
  1703. *
  1704. * Revision 1.85  2001/10/30 02:54:13  thiessen
  1705. * add Biostruc cache
  1706. *
  1707. * Revision 1.84  2001/10/17 17:46:22  thiessen
  1708. * save camera setup and rotation center in files
  1709. *
  1710. * Revision 1.83  2001/10/09 18:57:05  thiessen
  1711. * add CDD references editing dialog
  1712. *
  1713. * Revision 1.82  2001/10/08 00:00:09  thiessen
  1714. * estimate threader N random starts; edit CDD name
  1715. *
  1716. * Revision 1.81  2001/10/01 16:04:25  thiessen
  1717. * make CDD annotation window non-modal; add SetWindowTitle to viewers
  1718. *
  1719. * Revision 1.80  2001/09/27 20:58:28  thiessen
  1720. * add VisibleString filter option
  1721. *
  1722. * Revision 1.79  2001/09/27 15:37:59  thiessen
  1723. * decouple sequence import and BLAST
  1724. *
  1725. * Revision 1.78  2001/09/19 22:55:39  thiessen
  1726. * add preliminary net import and BLAST
  1727. *
  1728. * Revision 1.77  2001/09/18 03:10:45  thiessen
  1729. * add preliminary sequence import pipeline
  1730. *
  1731. * Revision 1.76  2001/08/27 00:06:23  thiessen
  1732. * add structure evidence to CDD annotation
  1733. *
  1734. * Revision 1.75  2001/08/13 22:30:59  thiessen
  1735. * add structure window mouse drag/zoom; add highlight option to render settings
  1736. *
  1737. * Revision 1.74  2001/08/10 15:01:57  thiessen
  1738. * fill out shortcuts; add update show/hide menu
  1739. *
  1740. * Revision 1.73  2001/08/09 19:07:13  thiessen
  1741. * add temperature and hydrophobicity coloring
  1742. *
  1743. * Revision 1.72  2001/07/19 19:14:38  thiessen
  1744. * working CDD alignment annotator ; misc tweaks
  1745. *
  1746. * Revision 1.71  2001/07/12 17:35:15  thiessen
  1747. * change domain mapping ; add preliminary cdd annotation GUI
  1748. *
  1749. * Revision 1.70  2001/07/10 16:39:55  thiessen
  1750. * change selection control keys; add CDD name/notes dialogs
  1751. *
  1752. * Revision 1.69  2001/07/04 19:39:17  thiessen
  1753. * finish user annotation system
  1754. *
  1755. * Revision 1.68  2001/06/29 18:13:58  thiessen
  1756. * initial (incomplete) user annotation system
  1757. *
  1758. * Revision 1.67  2001/06/21 02:02:34  thiessen
  1759. * major update to molecule identification and highlighting ; add toggle highlight (via alt)
  1760. *
  1761. * Revision 1.66  2001/06/15 14:06:40  thiessen
  1762. * save/load asn styles now complete
  1763. *
  1764. * Revision 1.65  2001/06/14 17:45:10  thiessen
  1765. * progress in styles<->asn ; add structure limits
  1766. *
  1767. * Revision 1.64  2001/06/14 00:34:01  thiessen
  1768. * asn additions
  1769. *
  1770. * Revision 1.63  2001/06/07 19:05:38  thiessen
  1771. * functional (although incomplete) render settings panel ; highlight title - not sequence - upon mouse click
  1772. *
  1773. * Revision 1.62  2001/06/05 13:21:08  thiessen
  1774. * fix structure alignment list problems
  1775. *
  1776. * Revision 1.61  2001/05/31 18:47:09  thiessen
  1777. * add preliminary style dialog; remove LIST_TYPE; add thread single and delete all; misc tweaks
  1778. *
  1779. * Revision 1.60  2001/05/17 18:32:47  thiessen
  1780. * fix sequence/molecule/master match bug
  1781. *
  1782. * Revision 1.59  2001/05/14 16:04:32  thiessen
  1783. * fix minor row reordering bug
  1784. *
  1785. * Revision 1.58  2001/05/11 18:20:11  thiessen
  1786. * fix loading of cdd with repeated master structure
  1787. *
  1788. * Revision 1.57  2001/05/02 13:46:28  thiessen
  1789. * major revision of stuff relating to saving of updates; allow stored null-alignments
  1790. *
  1791. * Revision 1.56  2001/04/17 20:15:39  thiessen
  1792. * load 'pending' Cdd alignments into update window
  1793. *
  1794. * Revision 1.55  2001/03/23 15:14:07  thiessen
  1795. * load sidechains in CDD's
  1796. *
  1797. * Revision 1.54  2001/03/23 04:18:52  thiessen
  1798. * parse and display disulfides
  1799. *
  1800. * Revision 1.53  2001/03/22 00:33:17  thiessen
  1801. * initial threading working (PSSM only); free color storage in undo stack
  1802. *
  1803. * Revision 1.52  2001/03/01 20:15:51  thiessen
  1804. * major rearrangement of sequence viewer code into base and derived classes
  1805. *
  1806. * Revision 1.51  2001/02/22 00:30:07  thiessen
  1807. * make directories global ; allow only one Sequence per StructureObject
  1808. *
  1809. * Revision 1.50  2001/02/16 00:40:01  thiessen
  1810. * remove unused sequences from asn data
  1811. *
  1812. * Revision 1.49  2001/02/13 01:03:57  thiessen
  1813. * backward-compatible domain ID's in output; add ability to delete rows
  1814. *
  1815. * Revision 1.48  2001/02/08 23:01:51  thiessen
  1816. * hook up C-toolkit stuff for threading; working PSSM calculation
  1817. *
  1818. * Revision 1.47  2001/02/02 20:17:33  thiessen
  1819. * can read in CDD with multi-structure but no struct. alignments
  1820. *
  1821. * Revision 1.46  2001/01/30 20:51:19  thiessen
  1822. * minor fixes
  1823. *
  1824. * Revision 1.45  2001/01/25 20:21:18  thiessen
  1825. * fix ostrstream memory leaks
  1826. *
  1827. * Revision 1.44  2001/01/18 19:37:28  thiessen
  1828. * save structure (re)alignments to asn output
  1829. *
  1830. * Revision 1.43  2001/01/12 01:34:18  thiessen
  1831. * network Biostruc load
  1832. *
  1833. * Revision 1.42  2000/12/29 19:23:40  thiessen
  1834. * save row order
  1835. *
  1836. * Revision 1.41  2000/12/22 19:26:40  thiessen
  1837. * write cdd output files
  1838. *
  1839. * Revision 1.40  2000/12/21 23:42:16  thiessen
  1840. * load structures from cdd's
  1841. *
  1842. * Revision 1.39  2000/12/20 23:47:48  thiessen
  1843. * load CDD's
  1844. *
  1845. * Revision 1.38  2000/12/15 15:51:47  thiessen
  1846. * show/hide system installed
  1847. *
  1848. * Revision 1.37  2000/12/01 19:35:57  thiessen
  1849. * better domain assignment; basic show/hide mechanism
  1850. *
  1851. * Revision 1.36  2000/11/30 15:49:40  thiessen
  1852. * add show/hide rows; unpack sec. struc. and domain features
  1853. *
  1854. * Revision 1.35  2000/11/13 18:06:53  thiessen
  1855. * working structure re-superpositioning
  1856. *
  1857. * Revision 1.34  2000/11/12 04:03:00  thiessen
  1858. * working file save including alignment edits
  1859. *
  1860. * Revision 1.33  2000/11/11 21:15:55  thiessen
  1861. * create Seq-annot from BlockMultipleAlignment
  1862. *
  1863. * Revision 1.32  2000/11/02 16:56:03  thiessen
  1864. * working editor undo; dynamic slave transforms
  1865. *
  1866. * Revision 1.31  2000/09/20 22:22:29  thiessen
  1867. * working conservation coloring; split and center unaligned justification
  1868. *
  1869. * Revision 1.30  2000/09/15 19:24:22  thiessen
  1870. * allow repeated structures w/o different local id
  1871. *
  1872. * Revision 1.29  2000/09/11 22:57:33  thiessen
  1873. * working highlighting
  1874. *
  1875. * Revision 1.28  2000/09/11 14:06:29  thiessen
  1876. * working alignment coloring
  1877. *
  1878. * Revision 1.27  2000/09/11 01:46:16  thiessen
  1879. * working messenger for sequence<->structure window communication
  1880. *
  1881. * Revision 1.26  2000/08/30 19:48:43  thiessen
  1882. * working sequence window
  1883. *
  1884. * Revision 1.25  2000/08/29 04:34:27  thiessen
  1885. * working alignment manager, IBM
  1886. *
  1887. * Revision 1.24  2000/08/28 23:47:19  thiessen
  1888. * functional denseg and dendiag alignment parsing
  1889. *
  1890. * Revision 1.23  2000/08/28 18:52:42  thiessen
  1891. * start unpacking alignments
  1892. *
  1893. * Revision 1.22  2000/08/27 18:52:22  thiessen
  1894. * extract sequence information
  1895. *
  1896. * Revision 1.21  2000/08/21 19:31:48  thiessen
  1897. * add style consistency checking
  1898. *
  1899. * Revision 1.20  2000/08/21 17:22:38  thiessen
  1900. * add primitive highlighting for testing
  1901. *
  1902. * Revision 1.19  2000/08/17 14:24:06  thiessen
  1903. * added working StyleManager
  1904. *
  1905. * Revision 1.18  2000/08/13 02:43:01  thiessen
  1906. * added helix and strand objects
  1907. *
  1908. * Revision 1.17  2000/08/07 14:13:16  thiessen
  1909. * added animation frames
  1910. *
  1911. * Revision 1.16  2000/08/07 00:21:18  thiessen
  1912. * add display list mechanism
  1913. *
  1914. * Revision 1.15  2000/08/04 22:49:04  thiessen
  1915. * add backbone atom classification and selection feedback mechanism
  1916. *
  1917. * Revision 1.14  2000/08/03 15:12:23  thiessen
  1918. * add skeleton of style and show/hide managers
  1919. *
  1920. * Revision 1.13  2000/07/27 13:30:51  thiessen
  1921. * remove 'using namespace ...' from all headers
  1922. *
  1923. * Revision 1.12  2000/07/18 16:50:11  thiessen
  1924. * more friendly rotation center setting
  1925. *
  1926. * Revision 1.11  2000/07/18 00:06:00  thiessen
  1927. * allow arbitrary rotation center
  1928. *
  1929. * Revision 1.10  2000/07/17 22:37:18  thiessen
  1930. * fix vector_math typo; correctly set initial view
  1931. *
  1932. * Revision 1.9  2000/07/17 04:20:50  thiessen
  1933. * now does correct structure alignment transformation
  1934. *
  1935. * Revision 1.8  2000/07/16 23:19:11  thiessen
  1936. * redo of drawing system
  1937. *
  1938. * Revision 1.7  2000/07/12 02:00:15  thiessen
  1939. * add basic wxWindows GUI
  1940. *
  1941. * Revision 1.6  2000/07/11 13:45:31  thiessen
  1942. * add modules to parse chemical graph; many improvements
  1943. *
  1944. * Revision 1.5  2000/07/01 15:43:50  thiessen
  1945. * major improvements to StructureBase functionality
  1946. *
  1947. * Revision 1.4  2000/06/29 16:46:07  thiessen
  1948. * use NCBI streams correctly
  1949. *
  1950. * Revision 1.3  2000/06/29 14:35:06  thiessen
  1951. * new atom_set files
  1952. *
  1953. * Revision 1.2  2000/06/28 13:07:55  thiessen
  1954. * store alt conf ensembles
  1955. *
  1956. * Revision 1.1  2000/06/27 20:09:40  thiessen
  1957. * initial checkin
  1958. *
  1959. */