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

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

生物技术

开发平台：

C/C++

gene_finder.hpp：源码内容
  1. /*
  2.  * ===========================================================================
  3.  * PRODUCTION $Log: gene_finder.hpp,v $
  4.  * PRODUCTION Revision 1000.1  2003/11/21 21:31:53  gouriano
  5.  * PRODUCTION PRODUCTION: UPGRADED [ORIGINAL] Dev-tree R1.2
  6.  * PRODUCTION
  7.  * ===========================================================================
  8.  */
  10. #ifndef __GENE_FINDER__HPP
  11. #define __GENE_FINDER__HPP
  13. /*  $Id: gene_finder.hpp,v 1000.1 2003/11/21 21:31:53 gouriano Exp $
  14.  * ===========================================================================
  15.  *
  16.  *                            PUBLIC DOMAIN NOTICE
  17.  *               National Center for Biotechnology Information
  18.  *
  19.  *  This software/database is a "United States Government Work" under the
  20.  *  terms of the United States Copyright Act.  It was written as part of
  21.  *  the author's official duties as a United States Government employee and
  22.  *  thus cannot be copyrighted.  This software/database is freely available
  23.  *  to the public for use. The National Library of Medicine and the U.S.
  24.  *  Government have not placed any restriction on its use or reproduction.
  25.  *
  26.  *  Although all reasonable efforts have been taken to ensure the accuracy
  27.  *  and reliability of the software and data, the NLM and the U.S.
  28.  *  Government do not and cannot warrant the performance or results that
  29.  *  may be obtained by using this software or data. The NLM and the U.S.
  30.  *  Government disclaim all warranties, express or implied, including
  31.  *  warranties of performance, merchantability or fitness for any particular
  32.  *  purpose.
  33.  *
  34.  *  Please cite the author in any work or product based on this material.
  35.  *
  36.  * ===========================================================================
  37.  *
  38.  * Authors:  Alexandre Souvorov
  39.  *
  40.  * File Description:
  41.  *
  42.  */
  44. #include <corelib/ncbistd.hpp>
  45. #include <corelib/ncbi_limits.hpp>
  46. #include <algo/gnomon/gnomon_exception.hpp>
  47. #include <string>
  48. #include <vector>
  49. #include <list>
  50. #include <set>
  51. #include <algorithm>
  52. #include <math.h>
  55. BEGIN_NCBI_SCOPE
  58. typedef vector<char> CVec;
  59. typedef vector<int> IVec;
  60. typedef vector<double> DVec;
  62. struct IPair : pair<int,int> 
  63. {
  64.     IPair(int f, int s) : pair<int,int>(f,s) {}
  65.     bool operator<(const IPair& p) const { return second < p.first; }
  66.     bool operator>(const IPair& p) const { return first > p.second; }
  67.     bool Intersect(const IPair& p) const { return !(*this < p || *this > p); }
  68.     bool Include(int i) const { return (i >= first && i <= second); }
  69. };
  71. enum Nucleotides { nA, nC, nG, nT, nN };
  72. enum { Plus, Minus };
  74. extern const double BadScore;
  75. extern const double LnHalf;
  76. extern const double LnThree;
  77. extern const int toMinus[5];
  78. extern const int TooFarLen;
  79. extern const char* aa_table;
  81. template<int order> class MarkovChain
  82. {
  83. public:
  84.     void InitScore(CNcbiIfstream& from);
  85.     double Score(const int* seq) const { return next[*(seq-order)].Score(seq); }
  86. private:
  87.     typedef MarkovChain<order> Type;
  88.     friend class MarkovChain<order+1>;
  89.     void Init(CNcbiIfstream& from);
  90.     void Average(Type& mc0, Type& mc1, Type& mc2, Type& mc3);
  91.     void toScore();
  93.     MarkovChain<order-1> next[5];
  94. };
  96. template<> class MarkovChain<0>
  97. {
  98. public:
  99.     void InitScore(CNcbiIfstream& from);
  100.     double Score(const int* seq) const { return score[*seq]; }
  101. private:
  102.     typedef MarkovChain<0> Type;
  103.     friend class MarkovChain<1>;
  104.     void Init(CNcbiIfstream& from);
  105.     void Average(Type& mc0, Type& mc1, Type& mc2, Type& mc3);
  106.     void toScore();
  108.     double score[5];
  109. };
  111. template<int order> class MarkovChainArray
  112. {
  113. public:
  114.     void InitScore(int l, CNcbiIfstream& from);
  115.     double Score(const int* seq) const;
  116. private:
  117.     int length;
  118.     vector< MarkovChain<order> > mc;
  119. };
  121. class InputModel
  122. {
  123. public:
  124.     virtual ~InputModel() = 0;
  126. protected:
  127.     static void Error(const string& label) {
  128.         NCBI_THROW(CGnomonException, eGenericError, label);
  129.     }
  130.     static pair<int,int> FindContent(CNcbiIfstream& from, const string& label, int cgcontent);
  131. };
  133. //Terminal's score is located on the last position of the left state
  134. class Terminal : public InputModel
  135. {
  136. public:
  137.     int InExon() const { return inexon; }
  138.     int InIntron() const { return inintron; }
  139.     int Left() const { return left; }
  140.     int Right() const { return right; }
  141.     virtual double Score(const IVec& seq, int i) const = 0;
  142.     ~Terminal() {}
  144. protected:
  145.     int inexon, inintron, left, right;
  146. };
  148. class MDD_Donor : public Terminal
  149. {
  150. public:
  151.     MDD_Donor(const string& file, int cgcontent);
  152.     double Score(const IVec& seq, int i) const;
  154. private:
  155.     IVec position, consensus;
  156.     vector< MarkovChainArray<0> > matrix;
  157. };
  159. template<int order> class WAM_Donor : public Terminal
  160. {
  161. public:
  162.     WAM_Donor(const string& file, int cgcontent);
  163.     double Score(const IVec& seq, int i) const;
  165. private:
  166.     MarkovChainArray<order> matrix;
  167. };
  169. template<int order> class WAM_Acceptor : public Terminal
  170. {
  171. public:
  172.     WAM_Acceptor(const string& file, int cgcontent);
  173.     double Score(const IVec& seq, int i) const;
  175. private:
  176.     MarkovChainArray<order> matrix;
  178. };
  180. class WMM_Start : public Terminal
  181. {
  182. public:
  183.     WMM_Start(const string& file, int cgcontent);
  184.     double Score(const IVec& seq, int i) const;
  186. private:
  187.     MarkovChainArray<0> matrix;
  188. };
  190. class WAM_Stop : public Terminal
  191. {
  192. public:
  193.     WAM_Stop(const string& file, int cgcontent);
  194.     double Score(const IVec& seq, int i) const;
  196. private:
  197.     MarkovChainArray<1> matrix;
  198. };
  200. class CodingRegion : public InputModel
  201. {
  202. public:
  203.     virtual double Score(const IVec& seq, int i, int codonshift) const = 0;
  204.     ~CodingRegion() {}
  205. };
  207. template<int order> class MC3_CodingRegion : public CodingRegion
  208. {
  209. public:
  210.     MC3_CodingRegion(const string& file, int cgcontent);
  211.     double Score(const IVec& seq, int i, int codonshift) const;
  213. private:
  214.     MarkovChain<order> matrix[3];
  215. };
  217. class NonCodingRegion : public InputModel
  218. {
  219. public:
  220.     virtual double Score(const IVec& seq, int i) const = 0;
  221.     ~NonCodingRegion() {}
  222. };
  224. template<int order> class MC_NonCodingRegion : public NonCodingRegion
  225. {
  226. public:
  227.     MC_NonCodingRegion(const string& file, int cgcontent);
  228.     double Score(const IVec& seq, int i) const;
  230. private:
  231.     MarkovChain<order> matrix;
  232. };
  234. class NullRegion : public NonCodingRegion
  235. {
  236. public:
  237.     double Score(const IVec& seq, int i) const { return 0; };
  238. };
  240. class AlignVec : public vector<IPair>
  241. {
  242. public:
  243.     typedef vector<IPair>::iterator It;
  244.     typedef vector<IPair>::const_iterator ConstIt;
  246.     enum { Prot, EST, mRNA, RefSeq, RefSeqBest, Wall};
  247.     AlignVec(int s = Plus, int i = 0, int t = EST, IPair cdl = IPair(-1,-1)) : type(t), strand(s), id(i), 
  248.     limits(numeric_limits<int>::max(),0), cds_limits(cdl), score(BadScore) {}
  249.     void Insert(IPair p);
  250.     void Erase(int i);
  251.     IPair Limits() const { return limits; }
  252.     IPair CdsLimits() const { return cds_limits; }
  253.     void SetCdsLimits(IPair p) { cds_limits = p; }
  254.     bool Intersect(const AlignVec& a) const 
  255.     {
  256.         return limits.Intersect(a.limits); 
  257.     }
  258.     void SetStrand(int s) { strand = s; }
  259.     int Strand() const { return strand; }
  260.     void SetType(int t) { type = t; }
  261.     int Type() const { return type; }
  262.     void SetID(int i) { id = i; }
  263.     int ID() const { return id; }
  264.     bool operator<(const AlignVec& a) const { return limits < a.limits; }
  265.     void Init();
  266.     void SetScore(double s) { score = s; }
  267.     double Score() const { return score; }
  268.     void Extend(const AlignVec& a); 
  270. private:
  271.     int type, strand, id;
  272.     IPair limits, cds_limits;
  273.     double score;
  274. };
  276. CNcbiIstream& operator>>(CNcbiIstream& s, AlignVec& a);
  277. CNcbiOstream& operator<<(CNcbiOstream& s, const AlignVec& a);
  280. class CCluster : public list<AlignVec>
  281. {
  282. public:
  283.     typedef list<AlignVec>::iterator It;
  284.     typedef list<AlignVec>::const_iterator ConstIt;
  285.     CCluster(int f = numeric_limits<int>::max(), int s = 0, int t = AlignVec::Prot) : limits(f,s), type(t) {}
  286.     void Insert(const AlignVec& a);
  287.     void Insert(const CCluster& c);
  288.     IPair Limits() const { return limits; }
  289.     int Type() const { return type; }
  290.     bool operator<(const CCluster& c) const { return limits < c.limits; }
  291.     void Init(int first, int second, int t);
  293. private:
  294.     IPair limits;
  295.     int type;
  296. };
  298. CNcbiIstream& operator>>(CNcbiIstream& s, CCluster& c);
  299. CNcbiOstream& operator<<(CNcbiOstream& s, const CCluster& c);
  302. class CClusterSet : public set<CCluster>
  303. {
  304. public:
  305.     typedef set<CCluster>::iterator It;
  306.     typedef set<CCluster>::const_iterator ConstIt;
  307.     CClusterSet() {}
  308.     CClusterSet(string c) : contig(c) {}
  309.     void InsertAlignment(const AlignVec& a);
  310.     void InsertCluster(CCluster c);
  311.     string Contig() const { return contig; }
  312.     void Init(string cnt);
  314. private:
  315.     string contig;
  316. };
  318. CNcbiIstream& operator>>(CNcbiIstream& s, CClusterSet& cls);
  319. CNcbiOstream& operator<<(CNcbiOstream& s, const CClusterSet& cls);
  322. class CFrameShiftInfo
  323. {
  324. public:
  325.     CFrameShiftInfo(int l, bool is_i, char i_v = 0) : loc(l), is_insert(is_i), insert_value(i_v) {}
  326.     int Loc() const { return loc; }
  327.     bool IsInsertion() const { return is_insert; }
  328.     bool IsDeletion() const { return !is_insert; }
  329.     char InsertValue() const { return insert_value; }
  330.     bool operator<(const CFrameShiftInfo& fsi) const { return loc < fsi.loc; }
  331. private:
  332.     int loc;  // location for deletion, insertion before location
  333.     // incertion - when I INCERT a new base in sequence
  334.     bool is_insert;
  335.     char insert_value;
  336. };
  338. typedef vector<CFrameShiftInfo> TFrameShifts;
  341. class SeqScores
  342. {
  343. public:
  344.     SeqScores (Terminal& a, Terminal& d, Terminal& stt, Terminal& stp, 
  345.                CodingRegion& cr, NonCodingRegion& ncr, NonCodingRegion& ing, 
  346.                CVec& sequence, int from, int to, const CClusterSet& cls, 
  347.                const TFrameShifts& initial_fshifts, bool repeats, bool leftwall, 
  348.                bool rightwall, string cntg);
  350.     int Shift() const { return shift; }
  351.     int AcceptorNumber(int strand) const { return anum[strand]; }
  352.     int DonorNumber(int strand) const { return dnum[strand]; }
  353.     int StartNumber(int strand) const { return sttnum[strand]; }
  354.     int StopNumber(int strand) const { return stpnum[strand]; }
  355.     double AcceptorScore(int i, int strand) const { return ascr[strand][i]; }
  356.     double DonorScore(int i, int strand) const { return dscr[strand][i]; }
  357.     double StartScore(int i, int strand) const { return sttscr[strand][i]; }
  358.     double StopScore(int i, int strand) const { return stpscr[strand][i]; }
  359.     Terminal& Acceptor() const { return acceptor; }
  360.     Terminal& Donor() const { return donor; }
  361.     Terminal& Start() const { return start; }
  362.     Terminal& Stop() const { return stop; }
  363.     const CClusterSet& Alignments() const { return cluster_set; }
  364.     const TFrameShifts& SeqFrameShifts() const { return fshifts; }
  365.     string Contig() const { return contig; }
  366.     bool StopInside(int a, int b, int strand, int frame) const;
  367.     bool OpenCodingRegion(int a, int b, int strand, int frame) const;
  368.     double CodingScore(int a, int b, int strand, int frame) const;
  369.     bool OpenNonCodingRegion(int a, int b, int strand) const;
  370.     double NonCodingScore(int a, int b, int strand) const;
  371.     bool OpenIntergenicRegion(int a, int b) const;
  372.     bool InAlignment(int a, int b) const;
  373.     double IntergenicScore(int a, int b, int strand) const;
  374.     int SeqLen() const { return seq[0].size(); }
  375.     bool SplittedStop(int id, int ia, int strand, int ph) const 
  376.     { return (dsplit[strand][ph][id]&asplit[strand][ph][ia]) ? true : false; }
  377.     bool isStart(int i, int strand) const;
  378.     bool isStop(int i, int strand) const;
  379.     bool isAG(int i, int strand) const;
  380.     bool isGT(int i, int strand) const;
  381.     bool isConsensusIntron(int i, int j, int strand) const;
  382.     const int* SeqPtr(int i, int strand) const;
  383.     int SeqMap(int i, bool forwrd) const;  // maps new coordinates to old coordinates,
  384.     // if insertion gives next or previous point
  385.     // depending on forwrd
  386.     int RevSeqMap(int i, bool forwrd) const; // maps old coordinates to new coordinates, 
  387.     // if deletion gives next or previous point
  388.     // depending on forwrd
  390. private:
  391.     Terminal &acceptor, &donor, &start, &stop;
  392.     CodingRegion &cdr;
  393.     NonCodingRegion &ncdr, &intrg;
  394.     const CClusterSet& cluster_set;
  395.     TFrameShifts fshifts;
  396.     IVec seq[2], laststop[2][3], notinexon[2][3], notinintron[2], notining;
  397.     IVec seq_map, rev_seq_map;
  398.     DVec ascr[2], dscr[2], sttscr[2], stpscr[2], ncdrscr[2], ingscr[2], cdrscr[2][3];
  399.     IVec asplit[2][2], dsplit[2][2];
  400.     IVec inalign;
  401.     int anum[2], dnum[2], sttnum[2], stpnum[2];
  402.     int shift;
  403.     string contig;
  404. };
  406. struct StateScores
  407. {
  408.     double score,branch,length,region,term;
  409. };
  411. template<class State> StateScores CalcStateScores(const State& st)
  412. {
  413.     StateScores sc;
  415.     if(st.NoLeftEnd())
  416.     {
  417.         if(st.NoRightEnd()) sc.length = st.ThroughLengthScore();
  418.         else sc.length = st.InitialLengthScore();
  419.     }
  420.     else
  421.     {
  422.         if(st.NoRightEnd()) sc.length = st.ClosingLengthScore();
  423.         else sc.length = st.LengthScore();
  424.     }
  426.     sc.region = st.RgnScore();
  427.     sc.term = st.TermScore();
  428.     if(sc.term == BadScore) sc.term = 0;
  429.     sc.score = st.Score();
  430.     if(st.LeftState()) sc.score -= st.LeftState()->Score();
  431.     sc.branch = sc.score-sc.length-sc.region-sc.term;
  433.     return sc;
  434. }
  436. class Lorentz
  437. {
  438. public:
  439.     bool Init(CNcbiIstream& from, const string& label);
  440.     double Score(int l) const { return score[(l-1)/step]; }
  441.     double ClosingScore(int l) const;
  442.     int MinLen() const { return minl; }
  443.     int MaxLen() const { return maxl; }
  444.     double AvLen() const { return avlen; }
  445.     double Through(int seqlen) const;
  447. private:
  448.     int minl, maxl, step;
  449.     double A, L, avlen, lnthrough;
  450.     DVec score, clscore;
  451. };
  453. class HMM_State : public InputModel
  454. {
  455. public:
  456.     HMM_State(int strn, int point);
  457.     const HMM_State* LeftState() const { return leftstate; }
  458.     const Terminal* TerminalPtr() const { return terminal; }
  459.     void UpdateLeftState(const HMM_State& left) { leftstate = &left; }
  460.     void UpdateScore(double scr) { score = scr; }
  461.     int MaxLen() const { return numeric_limits<int>::max(); };
  462.     int MinLen() const;
  463.     bool StopInside() const { return false; }
  464.     bool InAlignment() const { return false; }
  465.     int Strand() const { return strand; }
  466.     bool isPlus() const { return (strand == Plus); }
  467.     bool isMinus() const { return (strand == Minus); }
  468.     double Score() const { return score; }
  469.     int Start() const { return leftstate ? leftstate->stop+1 : 0; }
  470.     bool NoRightEnd() const { return stop < 0; }
  471.     bool NoLeftEnd() const { return leftstate == 0; }
  472.     int Stop() const  { return NoRightEnd() ? seqscr->SeqLen()-1 : stop; }
  473.     int RegionStart() const;
  474.     int RegionStop() const;
  475.     virtual StateScores GetStateScores() const = 0;
  476.     virtual string GetStateName() const = 0;
  477.     static void SetSeqScores(const SeqScores& s) { seqscr = &s; }
  479. protected:
  480.     int stop, strand;
  481.     double score;
  482.     const HMM_State* leftstate;
  483.     const Terminal* terminal;
  485.     static const SeqScores* seqscr;
  486. };
  488. class Intron; 
  489. class Intergenic;
  491. class Exon : public HMM_State
  492. {
  493. public:
  494.     static void Init(const string& file, int cgcontent);
  495.     Exon(int strn, int point, int ph) : HMM_State(strn,point), phase(ph), 
  496.     prevexon(0), mscore(BadScore) {}
  497.     int Phase() const { return phase; }
  498.     bool StopInside() const;
  499.     bool OpenRgn() const;
  500.     double RgnScore() const;
  501.     double DenScore() const { return 0; }
  502.     double ThroughLengthScore() const { return BadScore; } 
  503.     double InitialLengthScore() const { return BadScore; }
  504.     double ClosingLengthScore() const { return BadScore; }
  505.     void UpdatePrevExon(const Exon& e);
  506.     double MScore() const { return mscore; }
  509. protected:
  510.     int phase;
  511.     const Exon* prevexon;
  512.     double mscore;
  514.     static double firstphase[3], internalphase[3][3];
  515.     static Lorentz firstlen, internallen, lastlen, singlelen; 
  516.     static bool initialised;
  517. };
  519. class SingleExon : public Exon
  520. {
  521. public:
  522.     ~SingleExon() {}
  523.     SingleExon(int strn, int point);
  524.     int MaxLen() const { return singlelen.MaxLen(); }
  525.     int MinLen() const { return singlelen.MinLen(); }
  526.     const SingleExon* PrevExon() const { return static_cast<const SingleExon*>(prevexon); }
  527.     double LengthScore() const; 
  528.     double TermScore() const;
  529.     double BranchScore(const HMM_State& next) const { return BadScore; }
  530.     double BranchScore(const Intergenic& next) const;
  531.     StateScores GetStateScores() const { return CalcStateScores(*this); }
  532.     string GetStateName() const { return "SingleExon"; }
  533. };
  535. class FirstExon : public Exon
  536. {
  537. public:
  538.     ~FirstExon() {}
  539.     FirstExon(int strn, int ph, int point);
  540.     int MaxLen() const { return firstlen.MaxLen(); }
  541.     int MinLen() const { return firstlen.MinLen(); }
  542.     const FirstExon* PrevExon() const { return static_cast<const FirstExon*>(prevexon); }
  543.     double LengthScore() const;
  544.     double TermScore() const;
  545.     double BranchScore(const HMM_State& next) const { return BadScore; }
  546.     double BranchScore(const Intron& next) const;
  547.     StateScores GetStateScores() const { return CalcStateScores(*this); }
  548.     string GetStateName() const { return "FirstExon"; }
  549. };
  551. class InternalExon : public Exon
  552. {
  553. public:
  554.     ~InternalExon() {}
  555.     InternalExon(int strn, int ph, int point);
  556.     int MaxLen() const { return internallen.MaxLen(); }
  557.     int MinLen() const { return internallen.MinLen(); }
  558.     const InternalExon* PrevExon() const { return static_cast<const InternalExon*>(prevexon); }
  559.     double LengthScore() const; 
  560.     double TermScore() const;
  561.     double BranchScore(const HMM_State& next) const { return BadScore; }
  562.     double BranchScore(const Intron& next) const;
  563.     StateScores GetStateScores() const { return CalcStateScores(*this); }
  564.     string GetStateName() const { return "InternalExon"; }
  565. };
  567. class LastExon : public Exon
  568. {
  569. public:
  570.     ~LastExon() {}
  571.     LastExon(int strn, int ph, int point);
  572.     int MaxLen() const { return lastlen.MaxLen(); }
  573.     int MinLen() const { return lastlen.MinLen(); }
  574.     const LastExon* PrevExon() const { return static_cast<const LastExon*>(prevexon); }
  575.     double LengthScore() const; 
  576.     double TermScore() const;
  577.     double BranchScore(const HMM_State& next) const { return BadScore; }
  578.     double BranchScore(const Intergenic& next) const;
  579.     StateScores GetStateScores() const { return CalcStateScores(*this); }
  580.     string GetStateName() const { return "LastExon"; }
  581. };
  583. class Intron : public HMM_State
  584. {
  585. public:
  586.     static void Init(const string& file, int cgcontent, int seqlen);
  588.     ~Intron() {}
  589.     Intron(int strn, int ph, int point);
  590.     static int MinIntron() { return intronlen.MinLen(); }   // used for introducing frameshifts
  591.     int MinLen() const { return intronlen.MinLen(); }
  592.     int MaxLen() const { return intronlen.MaxLen(); }
  593.     int Phase() const { return phase; }
  594.     bool OpenRgn() const;
  595.     double RgnScore() const { return 0; }   // Intron scores are substructed from all others
  596.     double TermScore() const;
  597.     double DenScore() const { return lnDen[Phase()]; }
  598.     double LengthScore() const;
  599.     double ClosingLengthScore() const;
  600.     double ThroughLengthScore() const  { return lnThrough[Phase()]; }
  601.     double InitialLengthScore() const;
  602.     double BranchScore(const HMM_State& next) const { return BadScore; }
  603.     double BranchScore(const LastExon& next) const;
  604.     double BranchScore(const InternalExon& next) const;
  605.     bool SplittedStop() const;
  607.     StateScores GetStateScores() const { return CalcStateScores(*this); }
  608.     string GetStateName() const { return "Intron"; }
  610. protected:
  611.     int phase;
  612.     static double lnThrough[3], lnDen[3];
  613.     static double lnTerminal, lnInternal;
  614.     static Lorentz intronlen;
  615.     static bool initialised;
  616. };
  618. class Intergenic : public HMM_State
  619. {
  620. public:
  621.     static void Init(const string& file, int cgcontent, int seqlen);
  623.     ~Intergenic() {}
  624.     Intergenic(int strn, int point);
  625.     bool OpenRgn() const;
  626.     double RgnScore() const;
  627.     double TermScore() const;
  628.     double DenScore() const { return lnDen; }
  629.     double LengthScore() const { return intergeniclen.Score(Stop()-Start()+1); }
  630.     double ClosingLengthScore() const { return intergeniclen.ClosingScore(Stop()-Start()+1); }
  631.     double ThroughLengthScore() const { return lnThrough; }
  632.     double InitialLengthScore() const { return lnDen+ClosingLengthScore(); }
  633.     double BranchScore(const HMM_State& next) const { return BadScore; }
  634.     double BranchScore(const FirstExon& next) const; 
  635.     double BranchScore(const SingleExon& next) const;
  636.     bool InAlignment() const; 
  637.     StateScores GetStateScores() const { return CalcStateScores(*this); }
  638.     string GetStateName() const { return "Intergenic"; }
  640. protected:
  641.     static double lnThrough, lnDen;
  642.     static double lnSingle, lnMulti;
  643.     static Lorentz intergeniclen;
  644.     static bool initialised;
  645. };
  647. template<class T> class ParseVec : public vector<T>
  648. {
  649. public:
  650.     ParseVec() : num(-1) {}
  651.     int num;
  652. };
  654. class Gene;
  656. class Parse
  657. {
  658. public:
  659.     Parse(const SeqScores& ss);
  660.     const HMM_State* Path() const { return path; }
  661.     int PrintGenes(CNcbiOstream& to = cout, CNcbiOstream& toprot = cout,
  662.                    bool complete = false) const;
  663.     void PrintInfo() const;
  664.     list<Gene> GetGenes() const;
  665.     //typedef list<Gene>::iterator GenIt;
  667. private:
  668.     const SeqScores& seqscr;
  669.     const HMM_State* path;
  671.     ParseVec<Intergenic> igplus, igminus;
  672.     ParseVec<Intron> inplus[3], inminus[3];
  673.     ParseVec<FirstExon> feplus[3], feminus;
  674.     ParseVec<InternalExon> ieplus[3], ieminus[3];
  675.     ParseVec<LastExon> leplus, leminus[3];
  676.     ParseVec<SingleExon> seplus, seminus;
  677. };
  679. class ExonData
  680. {
  681.     friend list<Gene> Parse::GetGenes() const;
  682. public:
  683.     ExonData(int stt, int stp, int tp) : start(stt), stop(stp), type(tp) {} 
  684.     int Start() const { return start; }
  685.     int Stop() const { return stop; }
  686.     int Type() const { return type; }
  687.     const set<int>& ChainID() const { return chain_id; }
  688.     const set<int>& ProtID() const { return prot_id; }
  689.     const TFrameShifts& ExonFrameShifts() const { return fshifts; }
  690.     bool Identical(const ExonData& ed) const { return (ed.start == start && ed.stop == stop); }
  691.     bool operator<(const ExonData& ed) const { return (stop < ed.start); }
  692.     enum {Cds, Utr};
  693. private:
  694.     int start, stop, type;
  695.     set<int> chain_id, prot_id;
  696.     TFrameShifts fshifts;
  697. };
  699. class Gene : public vector<ExonData>
  700. {
  701.     friend list<Gene> Parse::GetGenes() const;
  702. public:
  703.     Gene(int s, bool l = true, bool r = true, int csf = 0) : 
  704.         strand(s), leftend(l), rightend(r), cds_shift(csf) {}
  705.     int Strand() const { return strand; }
  706.     int CDS_Shift() const { return cds_shift; }
  707.     const IVec& CDS() const { return cds; }
  708.     bool LeftComplete() const { return leftend; }
  709.     bool RightComplete() const { return rightend; }
  710.     bool Complete() const { return (leftend && rightend); }
  711. private:
  712.     int strand, cds_shift;
  713.     bool leftend, rightend;
  714.     IVec cds;
  715. };
  717. void LogicalCheck(const HMM_State& st, const SeqScores& ss);
  720. END_NCBI_SCOPE
  722. #include "models.hpp"
  723. #include "states.hpp"
  725. /*
  726.  * ===========================================================================
  727.  * $Log: gene_finder.hpp,v $
  728.  * Revision 1000.1  2003/11/21 21:31:53  gouriano
  729.  * PRODUCTION: UPGRADED [ORIGINAL] Dev-tree R1.2
  730.  *
  731.  * Revision 1.2  2003/11/06 15:02:21  ucko
  732.  * Use iostream interface from ncbistre.hpp for GCC 2.95 compatibility.
  733.  *
  734.  * Revision 1.1  2003/10/24 15:07:25  dicuccio
  735.  * Initial revision
  736.  *
  737.  * ===========================================================================
  738.  */
  739. #endif  // __GENE_FINDER__HPP