gnomon.cpp
上传用户:yhdzpy8989
上传日期:2007-06-13
资源大小:13604k
文件大小:9k
源码类别:

生物技术

开发平台:

C/C++

  1. /*
  2.  * ===========================================================================
  3.  * PRODUCTION $Log: gnomon.cpp,v $
  4.  * PRODUCTION Revision 1000.2  2004/06/01 18:08:33  gouriano
  5.  * PRODUCTION PRODUCTION: UPGRADED [GCC34_MSVC7] Dev-tree R1.4
  6.  * PRODUCTION
  7.  * ===========================================================================
  8.  */
  9. /*  $Id: gnomon.cpp,v 1000.2 2004/06/01 18:08:33 gouriano Exp $
  10.  * ===========================================================================
  11.  *
  12.  *                            PUBLIC DOMAIN NOTICE
  13.  *               National Center for Biotechnology Information
  14.  *
  15.  *  This software/database is a "United States Government Work" under the
  16.  *  terms of the United States Copyright Act.  It was written as part of
  17.  *  the author's official duties as a United States Government employee and
  18.  *  thus cannot be copyrighted.  This software/database is freely available
  19.  *  to the public for use. The National Library of Medicine and the U.S.
  20.  *  Government have not placed any restriction on its use or reproduction.
  21.  *
  22.  *  Although all reasonable efforts have been taken to ensure the accuracy
  23.  *  and reliability of the software and data, the NLM and the U.S.
  24.  *  Government do not and cannot warrant the performance or results that
  25.  *  may be obtained by using this software or data. The NLM and the U.S.
  26.  *  Government disclaim all warranties, express or implied, including
  27.  *  warranties of performance, merchantability or fitness for any particular
  28.  *  purpose.
  29.  *
  30.  *  Please cite the author in any work or product based on this material.
  31.  *
  32.  * ===========================================================================
  33.  *
  34.  * Authors:  Mike DiCuccio
  35.  *
  36.  * File Description:
  37.  *
  38.  */
  39. #include <ncbi_pch.hpp>
  40. #include <algo/gnomon/gnomon.hpp>
  41. #include "gene_finder.hpp"
  42. #include <objects/seqloc/Seq_loc.hpp>
  43. #include <objects/seqloc/Seq_id.hpp>
  44. #include <objects/seqloc/Packed_seqint.hpp>
  45. #include <objects/seqloc/Seq_interval.hpp>
  46. #include <objects/seqfeat/Seq_feat.hpp>
  47. #include <objects/seqfeat/RNA_ref.hpp>
  48. #include <objects/seqfeat/Cdregion.hpp>
  49. BEGIN_NCBI_SCOPE
  50. USING_SCOPE(ncbi::objects);
  51. CGnomon::CGnomon()
  52. : m_Repeats(false)
  53. , m_ScanRange(CRange<TSeqPos>::GetWhole())
  54. {
  55. }
  56. CGnomon::~CGnomon()
  57. {
  58. }
  59. //
  60. // set our sequence information
  61. //
  62. void CGnomon::SetSequence(const vector<char>& seq)
  63. {
  64.     m_Seq = seq;
  65. }
  66. void CGnomon::SetSequence(const string& str)
  67. {
  68.     m_Seq.clear();
  69.     m_Seq.reserve(str.length());
  70.     ITERATE (string, iter, str) {
  71.         m_Seq.push_back(*iter);
  72.     }
  73. }
  74. void CGnomon::SetSequence(const CSeqVector& vec)
  75. {
  76.     CSeqVector my_vec(vec);
  77.     my_vec.SetCoding(CSeq_data::e_Iupacna);
  78.     m_Seq.clear();
  79.     m_Seq.reserve(my_vec.size());
  80.     CSeqVector_CI iter(my_vec);
  81.     for ( ;  iter;  ++iter) {
  82.         m_Seq.push_back(*iter);
  83.     }
  84. }
  85. //
  86. // set the model data file
  87. // we save this as a string, as it is read multiple times
  88. //
  89. void CGnomon::SetModelData(const string& data_file)
  90. {
  91.     m_ModelFile = data_file;
  92. }
  93. //
  94. // set the a priori information
  95. //
  96. void CGnomon::SetAprioriInfo(const string& file)
  97. {
  98.     m_Clusters.reset(new CClusterSet());
  99.     if (file.empty()) {
  100.         return;
  101.     }
  102.     CNcbiIfstream from(file.c_str());
  103.     from >> *m_Clusters;
  104. }
  105. //
  106. // set the frame shift information
  107. //
  108. void CGnomon::SetFrameShiftInfo(const string& file)
  109. {
  110.     m_Fshifts.reset(new TFrameShifts());
  111.     if (file.empty()) {
  112.         return;
  113.     }
  114.     CNcbiIfstream from(file.c_str());
  115.     int loc;
  116.     while(from >> loc)
  117.     {
  118.         char is_i;
  119.         char c = 'N';
  120.         from >> is_i;
  121.         if(is_i == '+') {
  122.             from >> c;
  123.         }
  124.         m_Fshifts->push_back(CFrameShiftInfo(loc, (is_i == '+'), c));
  125.     }
  126. }
  127. //
  128. // set the repeats flag
  129. //
  130. void CGnomon::SetRepeats(bool b)
  131. {
  132.     m_Repeats = b;
  133. }
  134. //
  135. // set our scan range
  136. //
  137. void CGnomon::SetScanRange(const CRange<TSeqPos>& range)
  138. {
  139.     m_ScanRange = range;
  140. }
  141. //
  142. // set the scan start position
  143. //
  144. void CGnomon::SetScanFrom(TSeqPos from)
  145. {
  146.     m_ScanRange.SetFrom(from);
  147. }
  148. //
  149. // set the scan stop position
  150. //
  151. void CGnomon::SetScanTo(TSeqPos to)
  152. {
  153.     m_ScanRange.SetTo(to);
  154. }
  155. //
  156. // helper function - convert a vector<IPair> into a compact CSeq_loc
  157. //
  158. static inline
  159. CRef<CSeq_loc> s_ExonDataToLoc(const vector<IPair>& vec,
  160.                                ENa_strand strand, CSeq_id& id)
  161. {
  162.     CRef<CSeq_loc> loc(new CSeq_loc());
  163.     CPacked_seqint::Tdata data;
  164.     ITERATE (vector<IPair>, iter, vec) {
  165.         CRef<CSeq_interval> ival(new CSeq_interval());
  166.         ival->SetId(id);
  167.         ival->SetStrand(strand);
  168.         ival->SetFrom(iter->first);
  169.         ival->SetTo  (iter->second);
  170.         data.push_back(ival);
  171.     }
  172.     if (data.size() == 1) {
  173.         loc->SetInt(*data.front());
  174.     } else {
  175.         loc->SetPacked_int().Set().swap(data);
  176.     }
  177.     return loc;
  178. }
  179. //
  180. // Run() performs all the housekeeping necessary for GNOMON to work
  181. //
  182. void CGnomon::Run(void)
  183. {
  184.     // compute the GC content of the sequence
  185.     TSeqPos gc_content = 0;
  186.     ITERATE (vector<char>, iter, m_Seq) {
  187.         int c = toupper(*iter);
  188.         if (c == 'C'  ||  c == 'G') {
  189.             ++gc_content;
  190.         }
  191.     }
  192.     gc_content = static_cast<int>
  193.         (gc_content*100.0 / (m_Seq.size()) + 0.4999999);
  194.     //
  195.     // determine our scan range
  196.     //
  197.     TSeqPos from = m_ScanRange.GetFrom();
  198.     TSeqPos to   = m_ScanRange.GetTo();
  199.     if (m_ScanRange == CRange<TSeqPos>::GetWhole()) {
  200.         from = 0;
  201.         to   = m_Seq.size();
  202.     }
  203.     // TSeqPos is unsigned, so from is always > 0
  204.     to = min(to, (TSeqPos)m_Seq.size());
  205.     if (from > to) {
  206.         swap (from, to);
  207.     }
  208.     //
  209.     // GNOMON setup
  210.     //
  211.     auto_ptr<Terminal> donorp;
  212.     if(m_ModelFile.find("Human.inp") == m_ModelFile.length() - 9) {
  213.         donorp.reset(new MDD_Donor(m_ModelFile, gc_content));
  214.     } else {
  215.         donorp.reset(new WAM_Donor<2>(m_ModelFile,gc_content));
  216.     }
  217.     WAM_Acceptor<2>       acceptor      (m_ModelFile, gc_content);
  218.     WMM_Start             start         (m_ModelFile, gc_content);
  219.     WAM_Stop              stop          (m_ModelFile, gc_content);
  220.     MC3_CodingRegion<5>   cdr           (m_ModelFile, gc_content);
  221.     MC_NonCodingRegion<5> intron_reg    (m_ModelFile, gc_content);
  222.     MC_NonCodingRegion<5> intergenic_reg(m_ModelFile, gc_content);
  223.     Intron::Init    (m_ModelFile, gc_content, to - from + 1);
  224.     Intergenic::Init(m_ModelFile, gc_content, to - from + 1);
  225.     Exon::Init      (m_ModelFile, gc_content);
  226.     bool leftwall = true;
  227.     bool rightwall = true;
  228.     if ( !m_Clusters.get() ) {
  229.         m_Clusters.reset(new CClusterSet());
  230.     }
  231.     if ( !m_Fshifts.get() ) {
  232.         m_Fshifts.reset(new TFrameShifts());
  233.     }
  234.     SeqScores ss(acceptor, *donorp, start, stop, cdr, intron_reg, 
  235.                  intergenic_reg, m_Seq, from, to - from + 1,
  236.                  *m_Clusters, *m_Fshifts, m_Repeats,
  237.                  leftwall, rightwall, m_ModelFile);
  238.     HMM_State::SetSeqScores(ss);
  239.     m_Annot.Reset(new CSeq_annot());
  240.     m_Annot->AddName("GNOMON gene scan output");
  241.     Parse parse(ss);
  242.     list<Gene> genes = parse.GetGenes();
  243.     CRef<CSeq_id> id(new CSeq_id("lcl|gnomon"));
  244.     CSeq_annot::C_Data::TFtable& ftable = m_Annot->SetData().SetFtable();
  245.     ITERATE (list<Gene>, it, genes) {
  246.         const Gene& igene = *it;
  247.         int strand = igene.Strand();
  248.         vector<IPair> mrna_vec;
  249.         vector<IPair> cds_vec;
  250.         for (int j = 0;  j < igene.size();  ++j) {
  251.             const ExonData& exon = igene[j];
  252.             int a = exon.Start();
  253.             int b = exon.Stop();
  254.             if (j == 0  ||  igene[j-1].Stop()+1 != a) {
  255.                 // new exon
  256.                 mrna_vec.push_back(IPair(a,b));
  257.             } else {
  258.                 // attaching cds-part to utr-part
  259.                 mrna_vec.back().second = b;
  260.             } 
  261.             if (exon.Type() == ExonData::Cds) {
  262.                 cds_vec.push_back(IPair(a,b));
  263.             }
  264.         }
  265.         // stop-codon removal from cds
  266.         if (strand == Plus  &&  igene.RightComplete()) {
  267.             cds_vec.back().second -= 3;
  268.         }
  269.         if (strand == Minus  &&  igene.LeftComplete()) {
  270.             cds_vec.front().first += 3;
  271.         }
  272.         //
  273.         // create our mRNA
  274.         CRef<CSeq_feat> feat_mrna(new CSeq_feat());
  275.         feat_mrna->SetData().SetRna().SetType(CRNA_ref::eType_mRNA);
  276.         feat_mrna->SetLocation
  277.             (*s_ExonDataToLoc(mrna_vec,
  278.                               (strand == Plus ? eNa_strand_plus : eNa_strand_minus), *id));
  279.         //
  280.         // create the accompanying CDS
  281.         CRef<CSeq_feat> feat_cds(new CSeq_feat);
  282.         feat_cds->SetData().SetCdregion();
  283.         feat_cds->SetLocation
  284.             (*s_ExonDataToLoc(cds_vec,
  285.                               (strand == Plus ? eNa_strand_plus : eNa_strand_minus), *id));
  286.         ftable.push_back(feat_mrna);
  287.         ftable.push_back(feat_cds);
  288.     }
  289. }
  290. CRef<CSeq_annot> CGnomon::GetAnnot(void) const
  291. {
  292.     return m_Annot;
  293. }
  294. END_NCBI_SCOPE
  295. /*
  296.  * ===========================================================================
  297.  * $Log: gnomon.cpp,v $
  298.  * Revision 1000.2  2004/06/01 18:08:33  gouriano
  299.  * PRODUCTION: UPGRADED [GCC34_MSVC7] Dev-tree R1.4
  300.  *
  301.  * Revision 1.4  2004/05/21 21:41:03  gorelenk
  302.  * Added PCH ncbi_pch.hpp
  303.  *
  304.  * Revision 1.3  2003/11/06 00:51:05  ucko
  305.  * Adjust usage of min for platforms with 64-bit size_t.
  306.  *
  307.  * Revision 1.2  2003/11/05 21:15:24  ucko
  308.  * Fixed usage of ITERATE.  (m_Seq is a vector, not a string.)
  309.  *
  310.  * Revision 1.1  2003/10/24 15:07:25  dicuccio
  311.  * Initial revision
  312.  *
  313.  * ===========================================================================
  314.  */