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

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

生物技术

开发平台：

C/C++

cav_seqset.cpp：源码内容
  1. /*
  2.  * ===========================================================================
  3.  * PRODUCTION $Log: cav_seqset.cpp,v $
  4.  * PRODUCTION Revision 1000.2  2004/06/01 19:41:21  gouriano
  5.  * PRODUCTION PRODUCTION: UPGRADED [GCC34_MSVC7] Dev-tree R1.5
  6.  * PRODUCTION
  7.  * ===========================================================================
  8.  */
  10. /*  $Id: cav_seqset.cpp,v 1000.2 2004/06/01 19:41:21 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. #include <ncbi_pch.hpp>
  44. #include <corelib/ncbistd.hpp>
  45. #include <corelib/ncbistre.hpp>
  47. #include <objects/seq/Bioseq.hpp>
  48. #include <objects/seqloc/Seq_id.hpp>
  49. #include <objects/seqloc/PDB_seq_id.hpp>
  50. #include <objects/seqloc/PDB_mol_id.hpp>
  51. #include <objects/general/Object_id.hpp>
  52. #include <objects/seqset/Bioseq_set.hpp>
  53. #include <objects/seq/Seq_inst.hpp>
  54. #include <objects/seq/Seq_data.hpp>
  55. #include <objects/seq/NCBIeaa.hpp>
  56. #include <objects/seq/IUPACaa.hpp>
  57. #include <objects/seq/NCBI4na.hpp>
  58. #include <objects/seq/NCBI8na.hpp>
  59. #include <objects/seq/NCBI2na.hpp>
  60. #include <objects/seq/IUPACna.hpp>
  61. #include <objects/seq/NCBIstdaa.hpp>
  62. #include <objects/seq/Seq_annot.hpp>
  63. #include <objects/general/Dbtag.hpp>
  64. #include <objects/seqloc/Textseq_id.hpp>
  65. #include <objects/seq/Seq_descr.hpp>
  66. #include <objects/seq/Seqdesc.hpp>
  67. #include <objects/seqblock/PDB_block.hpp>
  69. #include <memory>
  71. #include <objtools/cddalignview/cav_seqset.hpp>
  72. #include <objtools/cddalignview/cddalignview.h>
  75. BEGIN_NCBI_SCOPE
  76. USING_SCOPE(objects);
  78. void SequenceSet::UnpackSeqSet(const CBioseq_set& bss)
  79. {
  80.     CBioseq_set::TSeq_set::const_iterator q, qe = bss.GetSeq_set().end();
  81.     for (q=bss.GetSeq_set().begin(); q!=qe; ++q) {
  82.         if (q->GetObject().IsSeq()) {
  84.             // only store amino acid or nucleotide sequences
  85.             if (q->GetObject().GetSeq().GetInst().GetMol() != CSeq_inst::eMol_aa &&
  86.                 q->GetObject().GetSeq().GetInst().GetMol() != CSeq_inst::eMol_dna &&
  87.                 q->GetObject().GetSeq().GetInst().GetMol() != CSeq_inst::eMol_rna &&
  88.                 q->GetObject().GetSeq().GetInst().GetMol() != CSeq_inst::eMol_na)
  89.                 continue;
  91.             const Sequence *sequence = new Sequence(q->GetObject().GetSeq());
  92.             if (!sequence || sequence->Status() != CAV_SUCCESS) {
  93.                 status = sequence->Status();
  94.                 return;
  95.             }
  96.             sequences.push_back(sequence);
  98.         } else { // Bioseq-set
  99.             UnpackSeqSet(q->GetObject().GetSet());
  100.         }
  101.     }
  102. }
  104. void SequenceSet::UnpackSeqEntry(const CSeq_entry& seqEntry)
  105. {
  106.     if (seqEntry.IsSeq()) {
  107.         const Sequence *sequence = new Sequence(seqEntry.GetSeq());
  108.         if (!sequence || sequence->Status() != CAV_SUCCESS) {
  109.             status = sequence->Status();
  110.             return;
  111.         }
  112.         sequences.push_back(sequence);
  113.     } else { // Bioseq-set
  114.         UnpackSeqSet(seqEntry.GetSet());
  115.     }
  116. }
  118. SequenceSet::SequenceSet(const CSeq_entry& seqEntry) :
  119.     master(NULL), status(CAV_SUCCESS)
  120. {
  121.     UnpackSeqEntry(seqEntry);
  122.     ERR_POST(Info << "number of sequences: " << sequences.size());
  123. }
  125. SequenceSet::SequenceSet(const SeqEntryList& seqEntries) :
  126.     master(NULL), status(CAV_SUCCESS)
  127. {
  128.     SeqEntryList::const_iterator s, se = seqEntries.end();
  129.     for (s=seqEntries.begin(); s!=se; ++s)
  130.         UnpackSeqEntry(s->GetObject());
  131.     ERR_POST(Info << "number of sequences: " << sequences.size());
  132. }
  134. SequenceSet::~SequenceSet(void)
  135. {
  136.     SequenceList::iterator s, se = sequences.end();
  137.     for (s=sequences.begin(); s!=se; ++s) delete *s;
  138. }
  141. const int Sequence::NOT_SET = -1;
  143. #define FIRSTOF2(byte) (((byte) & 0xF0) >> 4)
  144. #define SECONDOF2(byte) ((byte) & 0x0F)
  146. static void StringFrom4na(const vector< char >& vec, string *str, bool isDNA)
  147. {
  148.     if (SECONDOF2(vec.back()) > 0)
  149.         str->resize(vec.size() * 2);
  150.     else
  151.         str->resize(vec.size() * 2 - 1);
  153.     // first, extract 4-bit values
  154.     int i;
  155.     for (i=0; i<vec.size(); ++i) {
  156.         str->at(2*i) = FIRSTOF2(vec[i]);
  157.         if (SECONDOF2(vec[i]) > 0) str->at(2*i + 1) = SECONDOF2(vec[i]);
  158.     }
  160.     // then convert 4-bit values to ascii characters
  161.     for (i=0; i<str->size(); ++i) {
  162.         switch (str->at(i)) {
  163.             case 1: str->at(i) = 'A'; break;
  164.             case 2: str->at(i) = 'C'; break;
  165.             case 4: str->at(i) = 'G'; break;
  166.             case 8: isDNA ? str->at(i) = 'T' : str->at(i) = 'U'; break;
  167.             default:
  168.                 str->at(i) = 'X';
  169.         }
  170.     }
  171. }
  173. #define FIRSTOF4(byte) (((byte) & 0xC0) >> 6)
  174. #define SECONDOF4(byte) (((byte) & 0x30) >> 4)
  175. #define THIRDOF4(byte) (((byte) & 0x0C) >> 2)
  176. #define FOURTHOF4(byte) ((byte) & 0x03)
  178. static void StringFrom2na(const vector< char >& vec, string *str, bool isDNA)
  179. {
  180.     str->resize(vec.size() * 4);
  182.     // first, extract 4-bit values
  183.     int i;
  184.     for (i=0; i<vec.size(); ++i) {
  185.         str->at(4*i) = FIRSTOF4(vec[i]);
  186.         str->at(4*i + 1) = SECONDOF4(vec[i]);
  187.         str->at(4*i + 2) = THIRDOF4(vec[i]);
  188.         str->at(4*i + 3) = FOURTHOF4(vec[i]);
  189.     }
  191.     // then convert 4-bit values to ascii characters
  192.     for (i=0; i<str->size(); ++i) {
  193.         switch (str->at(i)) {
  194.             case 0: str->at(i) = 'A'; break;
  195.             case 1: str->at(i) = 'C'; break;
  196.             case 2: str->at(i) = 'G'; break;
  197.             case 3: isDNA ? str->at(i) = 'T' : str->at(i) = 'U'; break;
  198.         }
  199.     }
  200. }
  202. static void StringFromStdaa(const vector < char >& vec, std::string *str)
  203. {
  204.     static const char *stdaaMap = "-ABCDEFGHIKLMNPQRSTVWXYZU*";
  206.     str->resize(vec.size());
  207.     for (int i=0; i<vec.size(); ++i)
  208.         str->at(i) = stdaaMap[vec[i]];
  209. }
  211. Sequence::Sequence(const CBioseq& bioseq) :
  212.     gi(NOT_SET), pdbChain(' '), mmdbLink(NOT_SET), status(CAV_ERROR_SEQUENCES),
  213.     bioseqASN(&bioseq)
  214. {
  215.     // get Seq-id info
  216.     CBioseq::TId::const_iterator s, se = bioseq.GetId().end();
  217.     for (s=bioseq.GetId().begin(); s!=se; ++s) {
  218.         if (s->GetObject().IsGi()) {
  219.             gi = s->GetObject().GetGi();
  220.         } else if (s->GetObject().IsPdb()) {
  221.             pdbID = s->GetObject().GetPdb().GetMol().Get();
  222.             if (s->GetObject().GetPdb().IsSetChain())
  223.                 pdbChain = s->GetObject().GetPdb().GetChain();
  224.         } else if (s->GetObject().IsLocal() && s->GetObject().GetLocal().IsStr()) {
  225.             pdbID = s->GetObject().GetLocal().GetStr();
  226.         } else if (s->GetObject().IsGenbank() && s->GetObject().GetGenbank().IsSetAccession()) {
  227.             accession = s->GetObject().GetGenbank().GetAccession();
  228.         } else if (s->GetObject().IsSwissprot() && s->GetObject().GetSwissprot().IsSetAccession()) {
  229.             accession = s->GetObject().GetSwissprot().GetAccession();
  230.         }
  231.     }
  232.     if (gi == NOT_SET && pdbID.size() == 0 && accession.size() == 0) {
  233.         ERR_POST(Error << "Sequence::Sequence() - can't parse SeqId");
  234.         return;
  235.     }
  237.     // try to get description from title or compound
  238.     if (bioseq.IsSetDescr()) {
  239.         CSeq_descr::Tdata::const_iterator d, de = bioseq.GetDescr().Get().end();
  240.         for (d=bioseq.GetDescr().Get().begin(); d!=de; ++d) {
  241.             if (d->GetObject().IsTitle()) {
  242.                 description = d->GetObject().GetTitle();
  243.                 break;
  244.             } else if (d->GetObject().IsPdb() && d->GetObject().GetPdb().GetCompound().size() > 0) {
  245.                 description = d->GetObject().GetPdb().GetCompound().front();
  246.                 break;
  247.             }
  248.         }
  249.     }
  251.     // get link to MMDB id - mainly for CDD's where Biostrucs have to be loaded separately
  252.     if (bioseq.IsSetAnnot()) {
  253.         CBioseq::TAnnot::const_iterator a, ae = bioseq.GetAnnot().end();
  254.         for (a=bioseq.GetAnnot().begin(); a!=ae; ++a) {
  255.             if (a->GetObject().GetData().IsIds()) {
  256.                 CSeq_annot::C_Data::TIds::const_iterator i, ie = a->GetObject().GetData().GetIds().end();
  257.                 for (i=a->GetObject().GetData().GetIds().begin(); i!=ie; ++i) {
  258.                     if (i->GetObject().IsGeneral() &&
  259.                         i->GetObject().GetGeneral().GetDb() == "mmdb" &&
  260.                         i->GetObject().GetGeneral().GetTag().IsId()) {
  261.                         mmdbLink = i->GetObject().GetGeneral().GetTag().GetId();
  262.                         break;
  263.                     }
  264.                 }
  265.                 if (i != ie) break;
  266.             }
  267.         }
  268.     }
  269.     if (mmdbLink != NOT_SET)
  270.         ERR_POST(Info << "sequence gi " << gi << ", PDB '" << pdbID << "' chain '" << (char) pdbChain <<
  271.             "', is from MMDB id " << mmdbLink);
  273.     // get sequence string
  274.     if (bioseq.GetInst().GetRepr() == CSeq_inst::eRepr_raw && bioseq.GetInst().IsSetSeq_data()) {
  276.         // protein formats
  277.         if (bioseq.GetInst().GetSeq_data().IsNcbieaa()) {
  278.             sequenceString = bioseq.GetInst().GetSeq_data().GetNcbieaa().Get();
  279.         } else if (bioseq.GetInst().GetSeq_data().IsIupacaa()) {
  280.             sequenceString = bioseq.GetInst().GetSeq_data().GetIupacaa().Get();
  281.         } else if (bioseq.GetInst().GetSeq_data().IsNcbistdaa()) {
  282.             StringFromStdaa(bioseq.GetInst().GetSeq_data().GetNcbistdaa().Get(), &sequenceString);
  283.         }
  285.         // nucleotide formats
  286.         else if (bioseq.GetInst().GetSeq_data().IsIupacna()) {
  287.             sequenceString = bioseq.GetInst().GetSeq_data().GetIupacna().Get();
  288.         } else if (bioseq.GetInst().GetSeq_data().IsNcbi4na()) {
  289.             StringFrom4na(bioseq.GetInst().GetSeq_data().GetNcbi4na().Get(), &sequenceString,
  290.                 (bioseq.GetInst().GetMol() == CSeq_inst::eMol_dna));
  291.         } else if (bioseq.GetInst().GetSeq_data().IsNcbi8na()) {  // same repr. for non-X as 4na
  292.             StringFrom4na(bioseq.GetInst().GetSeq_data().GetNcbi8na().Get(), &sequenceString,
  293.                 (bioseq.GetInst().GetMol() == CSeq_inst::eMol_dna));
  294.         } else if (bioseq.GetInst().GetSeq_data().IsNcbi2na()) {
  295.             StringFrom2na(bioseq.GetInst().GetSeq_data().GetNcbi2na().Get(), &sequenceString,
  296.                 (bioseq.GetInst().GetMol() == CSeq_inst::eMol_dna));
  297.             if (bioseq.GetInst().IsSetLength() && bioseq.GetInst().GetLength() < sequenceString.length())
  298.                 sequenceString.resize(bioseq.GetInst().GetLength());
  299.         }
  301.         else {
  302.             ERR_POST(Critical << "Sequence::Sequence() - sequence " << gi
  303.                 << ": confused by sequence string format");
  304.             return;
  305.         }
  306.         if (bioseq.GetInst().IsSetLength() && bioseq.GetInst().GetLength() != sequenceString.length()) {
  307.             ERR_POST(Critical << "Sequence::Sequence() - sequence string length mismatch");
  308.             return;
  309.         }
  310.     } else {
  311.         ERR_POST(Critical << "Sequence::Sequence() - sequence " << gi
  312.                 << ": confused by sequence representation");
  313.         return;
  314.     }
  316.     status = CAV_SUCCESS;
  317. }
  319. CSeq_id * Sequence::CreateSeqId(void) const
  320. {
  321.     CSeq_id *sid = new CSeq_id();
  322.     if (pdbID.size() > 0) {
  323.         sid->SetPdb().SetMol().Set(pdbID);
  324.         if (pdbChain != ' ') sid->SetPdb().SetChain(pdbChain);
  325.     } else { // use gi
  326.         sid->SetGi(gi);
  327.     }
  328.     return sid;
  329. }
  331. string Sequence::GetTitle(void) const
  332. {
  333.     CNcbiOstrstream oss;
  334.     if (pdbID.size() > 0) {
  335.         oss << pdbID;
  336.         if (pdbChain != ' ') {
  337.             oss <<  '_' << (char) pdbChain;
  338.         }
  339.     } else if (gi != NOT_SET)
  340.         oss << "gi " << gi;
  341.     else if (accession.size() > 0)
  342.         oss << "acc " << accession;
  343.     else
  344.         oss << '?';
  345.     oss << '';
  346. auto_ptr<char> chars(oss.str()); // deletes memory upon function return
  347.     return string(oss.str());
  348. }
  350. END_NCBI_SCOPE
  353. /*
  354. * ---------------------------------------------------------------------------
  355. * $Log: cav_seqset.cpp,v $
  356. * Revision 1000.2  2004/06/01 19:41:21  gouriano
  357. * PRODUCTION: UPGRADED [GCC34_MSVC7] Dev-tree R1.5
  358. *
  359. * Revision 1.5  2004/05/21 21:42:51  gorelenk
  360. * Added PCH ncbi_pch.hpp
  361. *
  362. * Revision 1.4  2004/05/08 10:56:20  thiessen
  363. * better handling of repeated sequences
  364. *
  365. * Revision 1.3  2004/03/15 18:51:27  thiessen
  366. * prefer prefix vs. postfix ++/-- operators
  367. *
  368. * Revision 1.2  2003/06/02 16:06:41  dicuccio
  369. * Rearranged src/objects/ subtree.  This includes the following shifts:
  370. *     - src/objects/asn2asn --> arc/app/asn2asn
  371. *     - src/objects/testmedline --> src/objects/ncbimime/test
  372. *     - src/objects/objmgr --> src/objmgr
  373. *     - src/objects/util --> src/objmgr/util
  374. *     - src/objects/alnmgr --> src/objtools/alnmgr
  375. *     - src/objects/flat --> src/objtools/flat
  376. *     - src/objects/validator --> src/objtools/validator
  377. *     - src/objects/cddalignview --> src/objtools/cddalignview
  378. * In addition, libseq now includes six of the objects/seq... libs, and libmmdb
  379. * replaces the three libmmdb? libs.
  380. *
  381. * Revision 1.1  2003/03/19 19:04:12  thiessen
  382. * move again
  383. *
  384. * Revision 1.2  2003/03/19 16:06:28  thiessen
  385. * add <memory>
  386. *
  387. * Revision 1.1  2003/03/19 05:33:43  thiessen
  388. * move to src/app/cddalignview
  389. *
  390. * Revision 1.9  2003/02/03 17:52:03  thiessen
  391. * move CVS Log to end of file
  392. *
  393. * Revision 1.8  2002/05/28 17:55:48  thiessen
  394. * cavpr
  395. *
  396. * Revision 1.7  2001/07/12 21:32:13  thiessen
  397. * update for swissprot accessions
  398. *
  399. * Revision 1.6  2001/01/29 18:13:34  thiessen
  400. * split into C-callable library + main
  401. *
  402. * Revision 1.5  2001/01/25 00:51:20  thiessen
  403. * add command-line args; can read asn data from stdin
  404. *
  405. * Revision 1.4  2001/01/23 20:42:01  thiessen
  406. * add description
  407. *
  408. * Revision 1.3  2001/01/23 17:34:12  thiessen
  409. * add HTML output
  410. *
  411. * Revision 1.2  2001/01/22 15:55:11  thiessen
  412. * correctly set up ncbi namespacing
  413. *
  414. * Revision 1.1  2001/01/22 13:15:24  thiessen
  415. * initial checkin
  416. *
  417. */