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

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

生物技术

开发平台：

C/C++

sequence_set.cpp：源码内容
  1. /*
  2.  * ===========================================================================
  3.  * PRODUCTION $Log: sequence_set.cpp,v $
  4.  * PRODUCTION Revision 1000.4  2004/06/01 18:29:07  gouriano
  5.  * PRODUCTION PRODUCTION: UPGRADED [GCC34_MSVC7] Dev-tree R1.68
  6.  * PRODUCTION
  7.  * ===========================================================================
  8.  */
  10. /*  $Id: sequence_set.cpp,v 1000.4 2004/06/01 18:29:07 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 sequences
  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/ncbistl.hpp>
  52. #include <objects/seqloc/Seq_id.hpp>
  53. #include <objects/seqloc/PDB_seq_id.hpp>
  54. #include <objects/seqloc/PDB_mol_id.hpp>
  55. #include <objects/general/Object_id.hpp>
  56. #include <objects/seqset/Bioseq_set.hpp>
  57. #include <objects/seq/Seq_inst.hpp>
  58. #include <objects/seq/Seq_data.hpp>
  59. #include <objects/seq/NCBIeaa.hpp>
  60. #include <objects/seq/IUPACaa.hpp>
  61. #include <objects/seq/NCBIstdaa.hpp>
  62. #include <objects/seq/NCBI4na.hpp>
  63. #include <objects/seq/NCBI8na.hpp>
  64. #include <objects/seq/NCBI2na.hpp>
  65. #include <objects/seq/IUPACna.hpp>
  66. #include <objects/seq/Seq_annot.hpp>
  67. #include <objects/general/Dbtag.hpp>
  68. #include <objects/seqloc/Textseq_id.hpp>
  69. #include <objects/seq/Seq_descr.hpp>
  70. #include <objects/seq/Seqdesc.hpp>
  71. #include <objects/seqblock/PDB_block.hpp>
  72. #include <objects/seqfeat/BioSource.hpp>
  73. #include <objects/seqfeat/Org_ref.hpp>
  75. #include <regex.h>  // regex from C-toolkit
  77. #include "sequence_set.hpp"
  78. #include "molecule.hpp"
  79. #include "structure_set.hpp"
  80. #include "cn3d_tools.hpp"
  81. #include "molecule_identifier.hpp"
  82. #include "messenger.hpp"
  84. USING_NCBI_SCOPE;
  85. USING_SCOPE(objects);
  88. BEGIN_SCOPE(Cn3D)
  90. static void UnpackSeqSet(CBioseq_set& bss, SequenceSet *parent, SequenceSet::SequenceList& seqlist)
  91. {
  92.     CBioseq_set::TSeq_set::iterator q, qe = bss.SetSeq_set().end();
  93.     for (q=bss.SetSeq_set().begin(); q!=qe; ++q) {
  94.         if (q->GetObject().IsSeq()) {
  96.             // only store amino acid or nucleotide sequences
  97.             if (q->GetObject().GetSeq().GetInst().GetMol() != CSeq_inst::eMol_aa &&
  98.                 q->GetObject().GetSeq().GetInst().GetMol() != CSeq_inst::eMol_dna &&
  99.                 q->GetObject().GetSeq().GetInst().GetMol() != CSeq_inst::eMol_rna &&
  100.                 q->GetObject().GetSeq().GetInst().GetMol() != CSeq_inst::eMol_na)
  101.                 continue;
  103.             const Sequence *sequence = new Sequence(parent, q->GetObject().SetSeq());
  104.             seqlist.push_back(sequence);
  106.         } else { // Bioseq-set
  107.             UnpackSeqSet(q->GetObject().SetSet(), parent, seqlist);
  108.         }
  109.     }
  110. }
  112. static void UnpackSeqEntry(CSeq_entry& seqEntry, SequenceSet *parent, SequenceSet::SequenceList& seqlist)
  113. {
  114.     if (seqEntry.IsSeq()) {
  115.         const Sequence *sequence = new Sequence(parent, seqEntry.SetSeq());
  116.         seqlist.push_back(sequence);
  117.     } else { // Bioseq-set
  118.         UnpackSeqSet(seqEntry.SetSet(), parent, seqlist);
  119.     }
  120. }
  122. SequenceSet::SequenceSet(StructureBase *parent, SeqEntryList& seqEntries) :
  123.     StructureBase(parent)
  124. {
  125.     SeqEntryList::iterator s, se = seqEntries.end();
  126.     for (s=seqEntries.begin(); s!=se; ++s)
  127.         UnpackSeqEntry(s->GetObject(), this, sequences);
  128.     TRACEMSG("number of sequences: " << sequences.size());
  129. }
  131. #define FIRSTOF2(byte) (((byte) & 0xF0) >> 4)
  132. #define SECONDOF2(byte) ((byte) & 0x0F)
  134. static void StringFrom4na(const vector< char >& vec, string *str, bool isDNA)
  135. {
  136.     if (SECONDOF2(vec.back()) > 0)
  137.         str->resize(vec.size() * 2);
  138.     else
  139.         str->resize(vec.size() * 2 - 1);
  141.     // first, extract 4-bit values
  142.     int i;
  143.     for (i=0; i<vec.size(); ++i) {
  144.         str->at(2*i) = FIRSTOF2(vec[i]);
  145.         if (SECONDOF2(vec[i]) > 0) str->at(2*i + 1) = SECONDOF2(vec[i]);
  146.     }
  148.     // then convert 4-bit values to ascii characters
  149.     for (i=0; i<str->size(); ++i) {
  150.         switch (str->at(i)) {
  151.             case 1: str->at(i) = 'A'; break;
  152.             case 2: str->at(i) = 'C'; break;
  153.             case 4: str->at(i) = 'G'; break;
  154.             case 8: isDNA ? str->at(i) = 'T' : str->at(i) = 'U'; break;
  155.             default:
  156.                 str->at(i) = 'X';
  157.         }
  158.     }
  159. }
  161. #define FIRSTOF4(byte) (((byte) & 0xC0) >> 6)
  162. #define SECONDOF4(byte) (((byte) & 0x30) >> 4)
  163. #define THIRDOF4(byte) (((byte) & 0x0C) >> 2)
  164. #define FOURTHOF4(byte) ((byte) & 0x03)
  166. static void StringFrom2na(const vector< char >& vec, string *str, bool isDNA)
  167. {
  168.     str->resize(vec.size() * 4);
  170.     // first, extract 4-bit values
  171.     int i;
  172.     for (i=0; i<vec.size(); ++i) {
  173.         str->at(4*i) = FIRSTOF4(vec[i]);
  174.         str->at(4*i + 1) = SECONDOF4(vec[i]);
  175.         str->at(4*i + 2) = THIRDOF4(vec[i]);
  176.         str->at(4*i + 3) = FOURTHOF4(vec[i]);
  177.     }
  179.     // then convert 4-bit values to ascii characters
  180.     for (i=0; i<str->size(); ++i) {
  181.         switch (str->at(i)) {
  182.             case 0: str->at(i) = 'A'; break;
  183.             case 1: str->at(i) = 'C'; break;
  184.             case 2: str->at(i) = 'G'; break;
  185.             case 3: isDNA ? str->at(i) = 'T' : str->at(i) = 'U'; break;
  186.         }
  187.     }
  188. }
  190. static void StringFromStdaa(const vector < char >& vec, string *str)
  191. {
  192.     static const char *stdaaMap = "-ABCDEFGHIKLMNPQRSTVWXYZU*";
  194.     str->resize(vec.size());
  195.     for (int i=0; i<vec.size(); ++i)
  196.         str->at(i) = stdaaMap[vec[i]];
  197. }
  199. Sequence::Sequence(SequenceSet *parent, ncbi::objects::CBioseq& bioseq) :
  200.     StructureBase(parent), bioseqASN(&bioseq), molecule(NULL), isProtein(false)
  201. {
  202.     int gi = MoleculeIdentifier::VALUE_NOT_SET, mmdbID = MoleculeIdentifier::VALUE_NOT_SET,
  203.         pdbChain = MoleculeIdentifier::VALUE_NOT_SET;
  204.     string pdbID, accession;
  206.     // get Seq-id info
  207.     CBioseq::TId::const_iterator s, se = bioseq.GetId().end();
  208.     for (s=bioseq.GetId().begin(); s!=se; ++s) {
  209.         if ((*s)->IsGi()) {
  210.             gi = (*s)->GetGi();
  211.         } else if ((*s)->IsPdb()) {
  212.             pdbID = (*s)->GetPdb().GetMol().Get();
  213.             if ((*s)->GetPdb().IsSetChain())
  214.                 pdbChain = (*s)->GetPdb().GetChain();
  215.             else
  216.                 pdbChain = ' ';
  217.         } else if ((*s)->IsLocal() && (*s)->GetLocal().IsStr()) {
  218.             accession = (*s)->GetLocal().GetStr();
  219.             // special case where local accession is actually a PDB chain + extra stuff
  220.             if (pdbID.size() == 0 && accession.size() >= 7 &&
  221.                     accession[4] == ' ' && accession[6] == ' ' && isalpha(accession[5])) {
  222.                 pdbID = accession.substr(0, 4);
  223.                 pdbChain = accession[5];
  224.                 accession.erase();
  225.             }
  226.         } else if ((*s)->IsGenbank() && (*s)->GetGenbank().IsSetAccession()) {
  227.             accession = (*s)->GetGenbank().GetAccession();
  228.         } else if ((*s)->IsSwissprot() && (*s)->GetSwissprot().IsSetAccession()) {
  229.             accession = (*s)->GetSwissprot().GetAccession();
  230.         } else if ((*s)->IsOther() && (*s)->GetOther().IsSetAccession()) {
  231.             accession = (*s)->GetOther().GetAccession();
  232.         } else if ((*s)->IsEmbl() && (*s)->GetEmbl().IsSetAccession()) {
  233.             accession = (*s)->GetEmbl().GetAccession();
  234.         } else if ((*s)->IsDdbj() && (*s)->GetDdbj().IsSetAccession()) {
  235.             accession = (*s)->GetDdbj().GetAccession();
  236.         }
  237.     }
  238.     if (gi == MoleculeIdentifier::VALUE_NOT_SET && pdbID.size() == 0 && accession.size() == 0) {
  239.         ERRORMSG("Sequence::Sequence() - can't parse SeqId");
  240.         return;
  241.     }
  243.     if (bioseq.IsSetDescr()) {
  244.         string defline, taxid;
  245.         CSeq_descr::Tdata::const_iterator d, de = bioseq.GetDescr().Get().end();
  246.         for (d=bioseq.GetDescr().Get().begin(); d!=de; ++d) {
  248.             // get "defline" from title or compound
  249.             if ((*d)->IsTitle()) {              // prefer title over compound
  250.                 defline = (*d)->GetTitle();
  251.             } else if (defline.size() == 0 && (*d)->IsPdb() && (*d)->GetPdb().GetCompound().size() > 0) {
  252.                 defline = (*d)->GetPdb().GetCompound().front();
  253.             }
  255.             // get taxonomy
  256.             if ((*d)->IsSource()) {
  257.                 if ((*d)->GetSource().GetOrg().IsSetTaxname())
  258.                     taxid = (*d)->GetSource().GetOrg().GetTaxname();
  259.                 else if ((*d)->GetSource().GetOrg().IsSetCommon())
  260.                     taxid = (*d)->GetSource().GetOrg().GetCommon();
  261.             }
  262.         }
  263.         if (taxid.size() > 0)
  264.             description = string("[") + taxid + ']';
  265.         if (defline.size() > 0) {
  266.             if (taxid.size() > 0) description += ' ';
  267.             description += defline;
  268.         }
  269.     }
  271.     // get link to MMDB id - mainly for CDD's where Biostrucs have to be loaded separately
  272.     if (bioseq.IsSetAnnot()) {
  273.         CBioseq::TAnnot::const_iterator a, ae = bioseq.GetAnnot().end();
  274.         for (a=bioseq.GetAnnot().begin(); a!=ae; ++a) {
  275.             if (a->GetObject().GetData().IsIds()) {
  276.                 CSeq_annot::C_Data::TIds::const_iterator i, ie = a->GetObject().GetData().GetIds().end();
  277.                 for (i=a->GetObject().GetData().GetIds().begin(); i!=ie; ++i) {
  278.                     if (i->GetObject().IsGeneral() &&
  279.                         i->GetObject().GetGeneral().GetDb() == "mmdb" &&
  280.                         i->GetObject().GetGeneral().GetTag().IsId()) {
  281.                         mmdbID = i->GetObject().GetGeneral().GetTag().GetId();
  282.                         break;
  283.                     }
  284.                 }
  285.                 if (i != ie) break;
  286.             }
  287.         }
  288.     }
  289.     if (mmdbID != MoleculeIdentifier::VALUE_NOT_SET)
  290.         TRACEMSG("sequence gi " << gi << ", PDB '" << pdbID << "' chain '" << (char) pdbChain <<
  291.             "', is from MMDB id " << mmdbID);
  293.     // get sequence string
  294.     if (bioseq.GetInst().GetRepr() == CSeq_inst::eRepr_raw && bioseq.GetInst().IsSetSeq_data()) {
  296.         // protein formats
  297.         if (bioseq.GetInst().GetSeq_data().IsNcbieaa()) {
  298.             sequenceString = bioseq.GetInst().GetSeq_data().GetNcbieaa().Get();
  299.             isProtein = true;
  300.         } else if (bioseq.GetInst().GetSeq_data().IsIupacaa()) {
  301.             sequenceString = bioseq.GetInst().GetSeq_data().GetIupacaa().Get();
  302.             isProtein = true;
  303.         } else if (bioseq.GetInst().GetSeq_data().IsNcbistdaa()) {
  304.             StringFromStdaa(bioseq.GetInst().GetSeq_data().GetNcbistdaa().Get(), &sequenceString);
  305.             isProtein = true;
  306.         }
  308.         // nucleotide formats
  309.         else if (bioseq.GetInst().GetSeq_data().IsIupacna()) {
  310.             sequenceString = bioseq.GetInst().GetSeq_data().GetIupacna().Get();
  311.             // convert 'T' to 'U' for RNA
  312.             if (bioseq.GetInst().GetMol() == CSeq_inst::eMol_rna) {
  313.                 for (int i=0; i<sequenceString.size(); ++i) {
  314.                     if (sequenceString[i] == 'T')
  315.                         sequenceString[i] = 'U';
  316.                 }
  317.             }
  318.         } else if (bioseq.GetInst().GetSeq_data().IsNcbi4na()) {
  319.             StringFrom4na(bioseq.GetInst().GetSeq_data().GetNcbi4na().Get(), &sequenceString,
  320.                 (bioseq.GetInst().GetMol() == CSeq_inst::eMol_dna));
  321.         } else if (bioseq.GetInst().GetSeq_data().IsNcbi8na()) {  // same repr. for non-X as 4na
  322.             StringFrom4na(bioseq.GetInst().GetSeq_data().GetNcbi8na().Get(), &sequenceString,
  323.                 (bioseq.GetInst().GetMol() == CSeq_inst::eMol_dna));
  324.         } else if (bioseq.GetInst().GetSeq_data().IsNcbi2na()) {
  325.             StringFrom2na(bioseq.GetInst().GetSeq_data().GetNcbi2na().Get(), &sequenceString,
  326.                 (bioseq.GetInst().GetMol() == CSeq_inst::eMol_dna));
  327.             if (bioseq.GetInst().IsSetLength() && bioseq.GetInst().GetLength() < sequenceString.length())
  328.                 sequenceString.resize(bioseq.GetInst().GetLength());
  329.         }
  331.         else {
  332.             ERRORMSG("Sequence::Sequence() - sequence " << gi
  333.                 << ": confused by sequence string format");
  334.             return;
  335.         }
  337.         // check length
  338.         if (bioseq.GetInst().IsSetLength() && bioseq.GetInst().GetLength() != sequenceString.length()) {
  339.             ERRORMSG("Sequence::Sequence() - sequence string length mismatch");
  340.             return;
  341.         }
  343.         // force uppercase
  344.         for (int i=0; i<sequenceString.length(); ++i)
  345.             sequenceString[i] = toupper(sequenceString[i]);
  347.     } else {
  348.         ERRORMSG("Sequence::Sequence() - sequence " << gi << ": confused by sequence representation");
  349.         return;
  350.     }
  352.     // get identifier (may be NULL if there's a problem!)
  353.     identifier = MoleculeIdentifier::GetIdentifier(this, pdbID, pdbChain, mmdbID, gi, accession);
  354. }
  356. void Sequence::AddMMDBAnnotTag(int mmdbID) const
  357. {
  358.     CBioseq::TAnnot::const_iterator a, ae = bioseqASN->GetAnnot().end();
  359.     CSeq_annot::C_Data::TIds::const_iterator i, ie;
  360.     bool found = false;
  361.     for (a=bioseqASN->GetAnnot().begin(); a!=ae; ++a) {
  362.         if ((*a)->GetData().IsIds()) {
  363.             for (i=(*a)->GetData().GetIds().begin(), ie=(*a)->GetData().GetIds().end(); i!=ie; ++i) {
  364.                 if ((*i)->IsGeneral() && (*i)->GetGeneral().GetDb() == "mmdb" &&
  365.                     (*i)->GetGeneral().GetTag().IsId())
  366.                 {
  367.                     found = true;
  368.                     TRACEMSG("mmdb link already present in sequence " << identifier->ToString());
  369.                     if ((*i)->GetGeneral().GetTag().GetId() != mmdbID ||
  370.                             (identifier->mmdbID != MoleculeIdentifier::VALUE_NOT_SET &&
  371.                                 identifier->mmdbID != mmdbID))
  372.                         ERRORMSG("Sequence::AddMMDBAnnotTag() - mmdbID mismatch");
  373.                     break;
  374.                 }
  375.             }
  376.         }
  377.         if (found) break;
  378.     }
  379.     if (!found) {
  380.         CRef < CSeq_id > seqid(new CSeq_id());
  381.         seqid->SetGeneral().SetDb("mmdb");
  382.         seqid->SetGeneral().SetTag().SetId(mmdbID);
  383.         CRef < CSeq_annot > annot(new CSeq_annot());
  384.         annot->SetData().SetIds().push_back(seqid);
  385.         (const_cast<Sequence*>(this))->bioseqASN->SetAnnot().push_back(annot);
  386.     }
  387. }
  389. CSeq_id * Sequence::CreateSeqId(void) const
  390. {
  391.     CSeq_id *sid = new CSeq_id();
  392.     FillOutSeqId(sid);
  393.     return sid;
  394. }
  396. void Sequence::FillOutSeqId(ncbi::objects::CSeq_id *sid) const
  397. {
  398.     sid->Reset();
  399.     CBioseq::TId::const_iterator i, ie = bioseqASN->GetId().end();
  401.     // use pdb id if present
  402.     for (i=bioseqASN->GetId().begin(); i!=ie; ++i) {
  403.         if ((*i)->IsPdb()) {
  404.             sid->Assign(**i);
  405.             return;
  406.         }
  407.     }
  409.     // use gi if present
  410.     for (i=bioseqASN->GetId().begin(); i!=ie; ++i) {
  411.         if ((*i)->IsGi()) {
  412.             sid->Assign(**i);
  413.             return;
  414.         }
  415.     }
  417.     // otherwise, just use the first one
  418.     if (bioseqASN->GetId().size() > 0)
  419.         sid->Assign(bioseqASN->GetId().front().GetObject());
  420.     else
  421.         ERRORMSG("Sequence::FillOutSeqId() - can't do Seq-id on sequence " << identifier->ToString());
  423.     // dangerous to create new Seq-id's...
  424. //    if (identifier->pdbID.size() > 0 && identifier->pdbChain != MoleculeIdentifier::VALUE_NOT_SET) {
  425. //        sid->SetPdb().SetMol().Set(identifier->pdbID);
  426. //        if (identifier->pdbChain != ' ') sid->SetPdb().SetChain(identifier->pdbChain);
  427. //    } else if (identifier->gi != MoleculeIdentifier::VALUE_NOT_SET) { // use gi
  428. //        sid->SetGi(identifier->gi);
  429. //    } else if (identifier->accession.size() > 0) {
  430. //        CObject_id *oid = new CObject_id();
  431. //        oid->SetStr(identifier->accession);
  432. //        sid->SetLocal(*oid);
  433. }
  435. void Sequence::AddCSeqId(SeqIdPtr *id, bool addAllTypes) const
  436. {
  437.     if (identifier->pdbID.size() > 0) {
  438.         PDBSeqIdPtr pdbid = PDBSeqIdNew();
  439.         pdbid->mol = StrSave(identifier->pdbID.c_str());
  440.         pdbid->chain = (Uint1) identifier->pdbChain;
  441.         ValNodeAddPointer(id, SEQID_PDB, pdbid);
  442.         if (!addAllTypes) return;
  443.     }
  444.     if (identifier->gi != MoleculeIdentifier::VALUE_NOT_SET) {
  445.         ValNodeAddInt(id, SEQID_GI, identifier->gi);
  446.         if (!addAllTypes) return;
  447.     }
  448.     if (identifier->accession.size() > 0) {
  449.         ObjectIdPtr local = ObjectIdNew();
  450.         local->str = StrSave(identifier->accession.c_str());
  451.         ValNodeAddPointer(id, SEQID_LOCAL, local);
  452.         if (!addAllTypes) return;
  453.     }
  454. }
  456. int Sequence::GetOrSetMMDBLink(void) const
  457. {
  458.     if (molecule) {
  459.         const StructureObject *object;
  460.         if (!molecule->GetParentOfType(&object)) return identifier->mmdbID;
  461.         if (identifier->mmdbID != MoleculeIdentifier::VALUE_NOT_SET &&
  462.             identifier->mmdbID != object->mmdbID)
  463.             ERRORMSG("Sequence::GetOrSetMMDBLink() - mismatched MMDB ID: identifier says "
  464.                 << identifier->mmdbID << ", StructureObject says " << object->mmdbID);
  465.         else
  466.             const_cast<MoleculeIdentifier*>(identifier)->mmdbID = object->mmdbID;
  467.     }
  468.     return identifier->mmdbID;
  469. }
  471. void Sequence::LaunchWebBrowserWithInfo(void) const
  472. {
  473.     string db = isProtein ? "Protein" : "Nucleotide";
  474.     string opt = isProtein ? "GenPept" : "GenBank";
  475.     CNcbiOstrstream oss;
  476.     oss << "http://www.ncbi.nlm.nih.gov/entrez/query.fcgi?cmd=Search&doptcmdl=" << opt
  477.         << "&db=" << db << "&term=";
  478.     // prefer gi's, since accessions can be outdated
  479.     if (identifier->gi != MoleculeIdentifier::VALUE_NOT_SET) {
  480.         oss << identifier->gi;
  481.     } else if (identifier->pdbID.size() > 0) {
  482.         oss << identifier->pdbID.c_str();
  483.         if (identifier->pdbChain != ' ')
  484.             oss << (char) identifier->pdbChain;
  485.     } else if (identifier->accession.size() > 0) {
  486.         if (identifier->accession == "query" || identifier->accession == "consensus") return;
  487.         oss << identifier->accession.c_str();
  488.     }
  489.     oss << '';
  490.     LaunchWebPage(oss.str());
  491.     delete oss.str();
  492. }
  494. static bool Prosite2Regex(const string& prosite, string *regex, int *nGroups)
  495. {
  496.     try {
  497.         // check allowed characters
  498.         static const string allowed = "-ABCDEFGHIJKLMNOPQRSTUVWXYZ0123456789[],(){}<>.";
  499.         int i;
  500.         for (i=0; i<prosite.size(); ++i)
  501.             if (allowed.find(toupper(prosite[i])) == string::npos) break;
  502.         if (i != prosite.size()) throw "invalid ProSite character";
  503.         if (prosite[prosite.size() - 1] != '.') throw "ProSite pattern must end with '.'";
  505.         // translate into real regex syntax;
  506.         regex->erase();
  507.         *nGroups = 0;
  509.         bool inGroup = false;
  510.         for (int i=0; i<prosite.size(); ++i) {
  512.             // handle grouping and termini
  513.             bool characterHandled = true;
  514.             switch (prosite[i]) {
  515.                 case '-': case '.': case '>':
  516.                     if (inGroup) {
  517.                         *regex += ')';
  518.                         inGroup = false;
  519.                     }
  520.                     if (prosite[i] == '>') *regex += '$';
  521.                     break;
  522.                 case '<':
  523.                     *regex += '^';
  524.                     break;
  525.                 default:
  526.                     characterHandled = false;
  527.                     break;
  528.             }
  529.             if (characterHandled) continue;
  530.             if (!inGroup && (
  531.                     (isalpha(prosite[i]) && toupper(prosite[i]) != 'X') ||
  532.                     prosite[i] == '[' || prosite[i] == '{')) {
  533.                 *regex += '(';
  534.                 (*nGroups)++;
  535.                 inGroup = true;
  536.             }
  538.             // translate syntax
  539.             switch (prosite[i]) {
  540.                 case '(':
  541.                     *regex += '{';
  542.                     break;
  543.                 case ')':
  544.                     *regex += '}';
  545.                     break;
  546.                 case '{':
  547.                     *regex += "[^";
  548.                     break;
  549.                 case '}':
  550.                     *regex += ']';
  551.                     break;
  552.                 case 'X': case 'x':
  553.                     *regex += '.';
  554.                     break;
  555.                 default:
  556.                     *regex += toupper(prosite[i]);
  557.                     break;
  558.             }
  559.         }
  560.     }
  562.     catch (const char *err) {
  563.         ERRORMSG("Prosite2Regex() - " << err);
  564.         return false;
  565.     }
  567.     return true;
  568. }
  570. bool Sequence::HighlightPattern(const string& prositePattern) const
  571. {
  572.     // setup regex syntax
  573.     reg_syntax_t newSyntax = RE_CONTEXT_INDEP_ANCHORS | RE_CONTEXT_INVALID_OPS | RE_INTERVALS |
  574.         RE_LIMITED_OPS | RE_NO_BK_BRACES | RE_NO_BK_PARENS | RE_NO_EMPTY_RANGES;
  575.     reg_syntax_t oldSyntax = re_set_syntax(newSyntax);
  577.     bool retval = true;
  578.     try {
  579.         // allocate structures
  580.         static re_pattern_buffer *patternBuffer = NULL;
  581.         static re_registers *registers = NULL;
  582.         if (!patternBuffer) {
  583.             // new pattern initialized to zero
  584.             patternBuffer = (re_pattern_buffer *) calloc(1, sizeof(re_pattern_buffer));
  585.             if (!patternBuffer) throw "can't allocate pattern buffer";
  586.             patternBuffer->fastmap = (char *) calloc(256, sizeof(char));
  587.             if (!patternBuffer->fastmap) throw "can't allocate fastmap";
  588.             registers = (re_registers *) calloc(1, sizeof(re_registers));
  589.             if (!registers) throw "can't allocate registers";
  590.             patternBuffer->regs_allocated = REGS_UNALLOCATED;
  591.         }
  593.         // update pattern buffer if not the same pattern as before
  594.         static string previousPrositePattern;
  595.         static int nGroups;
  596.         int i;
  597.         if (prositePattern != previousPrositePattern) {
  599.             // convert from ProSite syntax
  600.             string regexPattern;
  601.             if (!Prosite2Regex(prositePattern, &regexPattern, &nGroups))
  602.                 throw "error converting ProSite to regex syntax";
  604.             // create pattern buffer
  605.             TRACEMSG("compiling pattern '" << regexPattern << "'");
  606.             const char *error = re_compile_pattern(regexPattern.c_str(), regexPattern.size(), patternBuffer);
  607.             if (error) throw error;
  609.             // optimize pattern buffer
  610.             int err = re_compile_fastmap(patternBuffer);
  611.             if (err) throw "re_compile_fastmap internal error";
  613.             previousPrositePattern = prositePattern;
  614.         }
  616.         // do the search, finding all non-overlapping matches
  617.         int start = 0;
  618.         while (start < Length()) {
  620.             int result = re_search(patternBuffer, sequenceString.c_str(),
  621.                 Length(), start, Length(), registers);
  622.             if (result == -1)
  623.                 break;
  624.             else if (result == -2)
  625.                 throw "re_search internal error";
  627.             // re_search gives the longest hit, but we want the shortest; so try to find the
  628.             // shortest hit within the hit already found by limiting the length of the string
  629.             // allowed to be included in the search. (This isn't very efficient! but
  630.             // the regex API doesn't have an option for finding the shortest hit...)
  631.             int stringSize = start + 1;
  632.             while (stringSize <= Length()) {
  633.                 result = re_search(patternBuffer, sequenceString.c_str(),
  634.                     stringSize, start, stringSize, registers);
  635.                 if (result >= 0) break;
  636.                 ++stringSize;
  637.             }
  639.             // parse the match registers, highlight ranges
  640. //            TESTMSG("found match starting at " << identifier->ToString() << " loc " << result+1);
  641.             int lastMatched = result;
  642.             for (i=1; i<registers->num_regs; ++i) {
  643.                 int from = registers->start[i], to = registers->end[i] - 1;
  644.                 if (from >= 0 && to >= 0) {
  645.                     if (to > lastMatched) lastMatched = to;
  647.                     // highlight this ranage
  648. //                    TESTMSG("register " << i << ": from " << from+1 << " to " << to+1);
  649.                     GlobalMessenger()->AddHighlights(this, from, to);
  650.                 }
  651.             }
  653.             start = lastMatched + 1;
  654.         }
  656.     } catch (const char *err) {
  657.         ERRORMSG("Sequence::HighlightPattern() - " << err);
  658.         retval = false;
  659.     }
  661.     // cleanup
  662.     re_set_syntax(oldSyntax);
  664.     return retval;
  665. }
  667. END_SCOPE(Cn3D)
  669. /*
  670. * ---------------------------------------------------------------------------
  671. * $Log: sequence_set.cpp,v $
  672. * Revision 1000.4  2004/06/01 18:29:07  gouriano
  673. * PRODUCTION: UPGRADED [GCC34_MSVC7] Dev-tree R1.68
  674. *
  675. * Revision 1.68  2004/05/21 21:41:39  gorelenk
  676. * Added PCH ncbi_pch.hpp
  677. *
  678. * Revision 1.67  2004/05/21 17:29:51  thiessen
  679. * allow conversion of mime to cdd data
  680. *
  681. * Revision 1.66  2004/03/15 18:27:12  thiessen
  682. * prefer prefix vs. postfix ++/-- operators
  683. *
  684. * Revision 1.65  2004/03/01 22:56:09  thiessen
  685. * convert 'T' to 'U' for RNA with IUPACna
  686. *
  687. * Revision 1.64  2004/02/19 17:05:06  thiessen
  688. * remove cn3d/ from include paths; add pragma to disable annoying msvc warning
  689. *
  690. * Revision 1.63  2004/01/05 17:09:16  thiessen
  691. * abort import and warn if same accession different gi
  692. *
  693. * Revision 1.62  2003/11/26 20:37:54  thiessen
  694. * prefer gi for URLs
  695. *
  696. * Revision 1.61  2003/11/20 22:08:49  thiessen
  697. * update Entrez url
  698. *
  699. * Revision 1.60  2003/09/03 18:14:01  thiessen
  700. * hopefully fix Seq-id issues
  701. *
  702. * Revision 1.59  2003/08/30 14:01:15  thiessen
  703. * use existing CSeq_id instead of creating new one
  704. *
  705. * Revision 1.58  2003/07/14 18:37:08  thiessen
  706. * change GetUngappedAlignedBlocks() param types; other syntax changes
  707. *
  708. * Revision 1.57  2003/03/13 14:26:18  thiessen
  709. * add file_messaging module; split cn3d_main_wxwin into cn3d_app, cn3d_glcanvas, structure_window, cn3d_tools
  710. *
  711. * Revision 1.56  2003/02/03 19:20:05  thiessen
  712. * format changes: move CVS Log to bottom of file, remove std:: from .cpp files, and use new diagnostic macros
  713. *
  714. * Revision 1.55  2002/12/12 15:04:11  thiessen
  715. * fix for report launch on nucleotides
  716. *
  717. * Revision 1.54  2002/12/09 16:25:04  thiessen
  718. * allow negative score threshholds
  719. *
  720. * Revision 1.53  2002/11/22 19:54:29  thiessen
  721. * fixes for wxMac/OSX
  722. *
  723. * Revision 1.52  2002/11/19 21:19:44  thiessen
  724. * more const changes for objects; fix user vs default style bug
  725. *
  726. * Revision 1.51  2002/11/06 00:18:10  thiessen
  727. * fixes for new CRef/const rules in objects
  728. *
  729. * Revision 1.50  2002/10/11 17:21:39  thiessen
  730. * initial Mac OSX build
  731. *
  732. * Revision 1.49  2002/09/13 14:21:45  thiessen
  733. * finish hooking up browser launch on unix
  734. *
  735. * Revision 1.48  2002/09/06 18:56:48  thiessen
  736. * add taxonomy to description
  737. *
  738. * Revision 1.47  2002/09/05 17:14:14  thiessen
  739. * create new netscape window on unix
  740. *
  741. * Revision 1.46  2002/07/12 13:24:10  thiessen
  742. * fixes for PSSM creation to agree with cddumper/RPSBLAST
  743. *
  744. * Revision 1.45  2002/01/24 20:08:17  thiessen
  745. * fix local id problem
  746. *
  747. * Revision 1.44  2002/01/19 02:34:46  thiessen
  748. * fixes for changes in asn serialization API
  749. *
  750. * Revision 1.43  2002/01/11 15:49:01  thiessen
  751. * update for Mac CW7
  752. *
  753. * Revision 1.42  2001/11/30 14:59:55  thiessen
  754. * add netscape launch for Mac
  755. *
  756. * Revision 1.41  2001/11/30 14:02:05  thiessen
  757. * progress on sequence imports to single structures
  758. *
  759. * Revision 1.40  2001/11/27 16:26:09  thiessen
  760. * major update to data management system
  761. *
  762. * Revision 1.39  2001/09/27 15:37:59  thiessen
  763. * decouple sequence import and BLAST
  764. *
  765. * Revision 1.38  2001/08/21 01:10:31  thiessen
  766. * fix lit. launch for nucleotides
  767. *
  768. * Revision 1.37  2001/08/09 19:07:13  thiessen
  769. * add temperature and hydrophobicity coloring
  770. *
  771. * Revision 1.36  2001/07/24 15:02:59  thiessen
  772. * use ProSite syntax for pattern searches
  773. *
  774. * Revision 1.35  2001/07/23 20:09:23  thiessen
  775. * add regex pattern search
  776. *
  777. * Revision 1.34  2001/07/19 19:14:38  thiessen
  778. * working CDD alignment annotator ; misc tweaks
  779. *
  780. * Revision 1.33  2001/07/10 16:39:55  thiessen
  781. * change selection control keys; add CDD name/notes dialogs
  782. *
  783. * Revision 1.32  2001/06/21 02:02:34  thiessen
  784. * major update to molecule identification and highlighting ; add toggle highlight (via alt)
  785. *
  786. * Revision 1.31  2001/06/18 21:48:27  thiessen
  787. * add [PACC] to Entrez query for PDB id's
  788. *
  789. * Revision 1.30  2001/06/02 17:22:46  thiessen
  790. * fixes for GCC
  791. *
  792. * Revision 1.29  2001/05/31 18:47:09  thiessen
  793. * add preliminary style dialog; remove LIST_TYPE; add thread single and delete all; misc tweaks
  794. *
  795. * Revision 1.28  2001/05/31 14:32:03  thiessen
  796. * better netscape launch for unix
  797. *
  798. * Revision 1.27  2001/05/25 01:38:16  thiessen
  799. * minor fixes for compiling on SGI
  800. *
  801. * Revision 1.26  2001/05/24 13:32:32  thiessen
  802. * further tweaks for GTK
  803. *
  804. * Revision 1.25  2001/05/15 23:48:37  thiessen
  805. * minor adjustments to compile under Solaris/wxGTK
  806. *
  807. * Revision 1.24  2001/05/02 16:35:15  thiessen
  808. * launch entrez web page on sequence identifier
  809. *
  810. * Revision 1.23  2001/04/20 18:03:22  thiessen
  811. * add ncbistdaa parsing
  812. *
  813. * Revision 1.22  2001/04/18 15:46:53  thiessen
  814. * show description, length, and PDB numbering in status line
  815. *
  816. * Revision 1.21  2001/03/28 23:02:17  thiessen
  817. * first working full threading
  818. *
  819. * Revision 1.20  2001/02/16 00:40:01  thiessen
  820. * remove unused sequences from asn data
  821. *
  822. * Revision 1.19  2001/02/13 01:03:57  thiessen
  823. * backward-compatible domain ID's in output; add ability to delete rows
  824. *
  825. * Revision 1.18  2001/02/08 23:01:50  thiessen
  826. * hook up C-toolkit stuff for threading; working PSSM calculation
  827. *
  828. * Revision 1.17  2001/01/25 20:21:18  thiessen
  829. * fix ostrstream memory leaks
  830. *
  831. * Revision 1.16  2001/01/09 21:45:00  thiessen
  832. * always use pdbID as title if known
  833. *
  834. * Revision 1.15  2001/01/04 18:20:52  thiessen
  835. * deal with accession seq-id
  836. *
  837. * Revision 1.14  2000/12/29 19:23:39  thiessen
  838. * save row order
  839. *
  840. * Revision 1.13  2000/12/28 20:43:09  vakatov
  841. * Do not use "set" as identifier.
  842. * Do not include <strstream>. Use "CNcbiOstrstream" instead of "ostrstream".
  843. *
  844. * Revision 1.12  2000/12/21 23:42:15  thiessen
  845. * load structures from cdd's
  846. *
  847. * Revision 1.11  2000/12/20 23:47:48  thiessen
  848. * load CDD's
  849. *
  850. * Revision 1.10  2000/12/15 15:51:47  thiessen
  851. * show/hide system installed
  852. *
  853. * Revision 1.9  2000/11/17 19:48:13  thiessen
  854. * working show/hide alignment row
  855. *
  856. * Revision 1.8  2000/11/11 21:15:54  thiessen
  857. * create Seq-annot from BlockMultipleAlignment
  858. *
  859. * Revision 1.7  2000/09/15 19:24:22  thiessen
  860. * allow repeated structures w/o different local id
  861. *
  862. * Revision 1.6  2000/09/03 18:46:49  thiessen
  863. * working generalized sequence viewer
  864. *
  865. * Revision 1.5  2000/08/30 23:46:27  thiessen
  866. * working alignment display
  867. *
  868. * Revision 1.4  2000/08/30 19:48:41  thiessen
  869. * working sequence window
  870. *
  871. * Revision 1.3  2000/08/28 23:47:19  thiessen
  872. * functional denseg and dendiag alignment parsing
  873. *
  874. * Revision 1.2  2000/08/28 18:52:42  thiessen
  875. * start unpacking alignments
  876. *
  877. * Revision 1.1  2000/08/27 18:52:22  thiessen
  878. * extract sequence information
  879. *
  880. */