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

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

生物技术

开发平台：

C/C++

parse.cpp：源码内容
  1. /*
  2.  * ===========================================================================
  3.  * PRODUCTION $Log: parse.cpp,v $
  4.  * PRODUCTION Revision 1000.2  2004/06/01 18:08:37  gouriano
  5.  * PRODUCTION PRODUCTION: UPGRADED [GCC34_MSVC7] Dev-tree R1.3
  6.  * PRODUCTION
  7.  * ===========================================================================
  8.  */
  10. /*  $Id: parse.cpp,v 1000.2 2004/06/01 18:08:37 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:  Alexandre Souvorov
  36.  *
  37.  * File Description:
  38.  *
  39.  */
  41. #include <ncbi_pch.hpp>
  42. #include "gene_finder.hpp"
  44. BEGIN_NCBI_SCOPE
  46. template<class T> inline void PushState(ParseVec<T>* vecp, int strand, int point)
  47. {
  48.     for (int phase = 0;  phase < 3;  ++phase) {
  49.         vecp[phase].push_back(T(strand,phase,point));
  50.         ++vecp[phase].num;
  51.     }
  52. }
  55. template<class T> inline void PushState(ParseVec<T>& vec, int strand, int point)
  56. {
  57.     vec.push_back(T(strand,point));
  58.     ++vec.num;
  59. }
  62. inline void PushState(ParseVec<LastExon>& vec, int strand, int point)
  63. {
  64.     vec.push_back(LastExon(strand,2,point));
  65.     ++vec.num;
  66. }
  69. inline void PushState(ParseVec<FirstExon>& vec, int strand, int point)
  70. {
  71.     vec.push_back(FirstExon(strand,0,point));
  72.     ++vec.num;
  73. }
  76. template<class L, class R> inline 
  77. bool TooFar(const L& left, const R& right, double score) 
  78. { 
  79.     if (left.MScore() == BadScore) {
  80.         return true;
  81.     }
  82.     int len = right.Stop() - left.Stop();
  83.     return (len > TooFarLen  &&  score + left.MScore() < right.Score());
  84. }
  87. struct RState
  88. {
  89.     RState(const HMM_State& l, const HMM_State& r)
  90.     {
  91.         lsp = r.LeftState();
  92.         rsp = const_cast<HMM_State*>(&r);
  93.         rsp->UpdateLeftState(l);
  94.     }
  95.     ~RState() { rsp->UpdateLeftState(*lsp); } 
  97.     const HMM_State* lsp;
  98.     HMM_State* rsp;
  99. };
  101. template<class Left, class Right> inline
  102. bool EvaluateNewScore(const Left& left, const Right& right, double& rscore, 
  103.                       bool& openrgn, bool& inalign)
  104. {
  106.     rscore = BadScore;
  108.     RState rst(left,right);
  110.     int len = right.Stop() - left.Stop();
  111.     if (len > right.MaxLen()) {
  112.         return false;
  113.     }
  114.     if ( !right.NoRightEnd()  &&  len < right.MinLen() ) {
  115.         return true;
  116.     }
  118.     double scr, score = 0;
  119.     if (left.isPlus()) {
  120.         scr = left.BranchScore(right);
  121.         if (scr == BadScore) {
  122.             return true;
  123.         }
  124.     } else {
  125.         scr = right.BranchScore(left);
  126.         if (scr == BadScore) {
  127.             return true;
  128.         }
  129.         scr += right.DenScore() - left.DenScore();
  130.     }
  131.     score += scr;
  133.     // this check is frame - dependent and MUST be after BranchScore
  134.     if (right.StopInside()) {
  135.         return false;
  136.     }     
  138.     if (right.NoRightEnd()) {
  139.         scr = right.ClosingLengthScore();
  140.     } else {
  141.         scr = right.LengthScore();
  142.     }
  143.     if (scr == BadScore) {
  144.         return true;
  145.     }
  146.     score += scr;
  148.     scr = right.RgnScore();
  149.     if (scr == BadScore) {
  150.         return true;
  151.     }
  152.     score += scr;
  154.     if ( !right.NoRightEnd() ) {
  155.         scr = right.TermScore();
  156.         if (scr == BadScore) {
  157.             return true;
  158.         }
  159.         score += scr;
  160.     }
  162.     openrgn = right.OpenRgn();
  163.     inalign = right.InAlignment();
  164.     rscore = score;
  165.     return true;
  166. }
  169. class RightIsExon {};
  170. class RightIsNotExon {};
  172. template<class L, class R> // template for all right exons
  173. inline bool ForwardStep(const L& left, R& right, RightIsExon rie)
  174. {
  175.     double score;
  176.     bool openrgn, inalign;
  177.     if ( !EvaluateNewScore(left,right,score,openrgn,inalign) ) {
  178.         return false;
  179.     } else if (score == BadScore) {
  180.         return true;
  181.     }
  183.     if (left.Score() != BadScore  &&  openrgn  &&  !inalign) {
  184.         double scr = score + left.Score();
  185.         if (scr > right.Score()) {
  186.             right.UpdateLeftState(left);
  187.             right.UpdateScore(scr);
  188.         }
  189.     }
  191.     if ( !openrgn ) {
  192.         return false;
  193.     }
  195.     return true;
  196. }
  199. template<class L, class R> // template for all right non - exons
  200. inline bool ForwardStep(const L& left, R& right, RightIsNotExon rine)
  201. {
  202.     double score;
  203.     bool openrgn, inalign;
  204.     if ( !EvaluateNewScore(left,right,score,openrgn,inalign) ) {
  205.         return false;
  206.     } else if (score == BadScore) {
  207.         return true;
  208.     }
  210.     if (left.Score() != BadScore  &&  openrgn  &&  !inalign) {
  211.         double scr = score + left.Score();
  212.         if (scr > right.Score()) {
  213.             right.UpdateLeftState(left);
  214.             right.UpdateScore(scr);
  215.         }
  216.     }
  218.     if ( !openrgn  ||  TooFar(left, right, score) ) {
  219.         return false;
  220.     }
  222.     return true;
  223. }
  226. template<class L, class R> // template for all right exons
  227. inline void MakeStep(ParseVec<L>& lvec, ParseVec<R>& rvec, RightIsExon rie)
  228. {
  229.     if (lvec.num < 0) {
  230.         return;
  231.     }
  232.     int i = lvec.num;
  233.     if (lvec[i].Stop() == rvec[rvec.num].Stop()) {
  234.         --i;
  235.     }
  236.     while (i >= 0  &&  ForwardStep(lvec[i--],rvec[rvec.num],rie)) {
  237.     }
  239.     if (rvec.num > 0) {
  240.         rvec[rvec.num].UpdatePrevExon(rvec[rvec.num - 1]);
  241.     }
  242. }
  245. template<class L, class R> // template for all right non - exons
  246. inline void MakeStep(ParseVec<L>& lvec, ParseVec<R>& rvec, RightIsNotExon rine)
  247. {
  248.     if (lvec.num < 0) {
  249.         return;
  250.     }
  251.     R& right = rvec[rvec.num];
  252.     int i = lvec.num;
  253.     int rlimit = right.Stop();
  254.     if (lvec[i].Stop() == rlimit) {
  255.         --i;
  256.     }
  257.     int nearlimit = max(0,rlimit - TooFarLen);
  258.     while (i >= 0  &&  lvec[i].Stop() >= nearlimit) {
  259.         if ( !ForwardStep(lvec[i--],right,rine) ) {
  260.             return;
  261.         }
  262.     }
  263.     if (i < 0) {
  264.         return;
  265.     }
  266.     for (const L* p = &lvec[i];  p !=0;  p = p->PrevExon()) {
  267.         if ( !ForwardStep(*p,right,rine) ) {
  268.             return;
  269.         }
  270.     }
  271. }
  273. // R / L indicates that left / right parameter is a 3 - element array
  275. template<class L, class R, class D> 
  276. inline void MakeStepR(L& lvec, R rvec, D d)
  277. {
  278.     for (int phr = 0;  phr < 3;  ++phr) {
  279.         MakeStep(lvec, rvec[phr], d);
  280.     }
  281. }
  284. template<class L, class R, class D> 
  285. inline void MakeStepL(L lvec, R& rvec, D d)
  286. {
  287.     for (int phl = 0;  phl < 3;  ++phl) {
  288.         MakeStep(lvec[phl], rvec, d);
  289.     }
  290. }
  293. template<class L, class R, class D> 
  294. inline void MakeStepLR(L lvec, R rvec, D d)
  295. {
  296.     for (int phl = 0;  phl < 3;  ++phl) {
  297.         for (int phr = 0;  phr < 3;  ++phr)
  298.         {
  299.             MakeStep(lvec[phl], rvec[phr], d);
  300.         }
  301.     }
  302. }
  305. template<class L, class R, class D> 
  306. inline void MakeStepLR(L lvec, R rvec, D d, int shift)
  307. {
  308.     for (int phl = 0;  phl < 3;  ++phl) {
  309.         int phr = (shift + phl)%3;
  310.         MakeStep(lvec[phl], rvec[phr], d);
  311.     }
  312. }
  315. Parse::Parse(const SeqScores& ss) : seqscr(ss)
  316. {
  317.     try {
  318.         igplus.reserve (seqscr.StartNumber(Plus) + 1);
  319.         igminus.reserve(seqscr.StopNumber (Minus) + 1);
  320.         feminus.reserve(seqscr.StartNumber(Minus) + 1);
  321.         leplus.reserve (seqscr.StopNumber (Plus) + 1);
  322.         seplus.reserve (seqscr.StopNumber (Plus) + 1);
  323.         seminus.reserve(seqscr.StartNumber(Minus) + 1);
  324.         for (int phase = 0;  phase < 3;  ++phase) {
  325.             inplus[phase ].reserve(seqscr.AcceptorNumber(Plus) + 1);
  326.             inminus[phase].reserve(seqscr.DonorNumber   (Minus) + 1);
  327.             feplus[phase ].reserve(seqscr.DonorNumber   (Plus) + 1);
  328.             ieplus[phase ].reserve(seqscr.DonorNumber   (Plus) + 1);
  329.             ieminus[phase].reserve(seqscr.AcceptorNumber(Minus) + 1);
  330.             leminus[phase].reserve(seqscr.AcceptorNumber(Minus) + 1);
  331.         }
  332.     }
  333.     catch(bad_alloc&) {
  334.         NCBI_THROW(CGnomonException, eGenericError,
  335.                    "Not enough memory in Parse");
  336.     }
  338.     int len = seqscr.SeqLen();
  340.     RightIsExon rie;
  341.     RightIsNotExon rine;
  343.     for (int i = 0;  i < len;  ++i) {
  344.         if (seqscr.AcceptorScore(i,Plus) != BadScore) {
  345.             PushState (inplus,Plus,i);
  346.             MakeStepLR(ieplus,inplus,rine,1);
  347.             MakeStepLR(feplus,inplus,rine,1);
  348.         }
  350.         if (seqscr.AcceptorScore(i,Minus) != BadScore) {
  351.             PushState (ieminus,Minus,i);
  352.             PushState (leminus,Minus,i);
  353.             MakeStepLR(inminus,ieminus,rie);
  354.             MakeStepR(igminus,leminus,rie);
  355.         }
  357.         if (seqscr.DonorScore(i,Plus) != BadScore) {
  358.             PushState (feplus,Plus,i);
  359.             PushState (ieplus,Plus,i);
  360.             MakeStepLR(inplus,ieplus,rie);
  361.             MakeStepR(igplus,feplus,rie);
  362.         }
  364.         if (seqscr.DonorScore(i,Minus) != BadScore) {
  365.             PushState (inminus,Minus,i);
  366.             MakeStepLR(leminus,inminus,rine,0);
  367.             MakeStepLR(ieminus,inminus,rine,0);
  368.         }
  370.         if (seqscr.StartScore(i,Plus) != BadScore) {
  371.             PushState(igplus,Plus,i);
  372.             MakeStep (seplus,igplus,rine);
  373.             MakeStep (seminus,igplus,rine);
  374.             MakeStep (leplus,igplus,rine);
  375.             MakeStep (feminus,igplus,rine);
  376.         }
  378.         if (seqscr.StartScore(i,Minus) != BadScore) {
  379.             PushState(feminus,Minus,i);
  380.             PushState(seminus,Minus,i);
  381.             MakeStepL(inminus,feminus,rie);
  382.             MakeStep (igminus,seminus,rie);
  383.         }
  385.         if (seqscr.StopScore(i,Plus) != BadScore) {
  386.             PushState(leplus,Plus,i);
  387.             PushState(seplus,Plus,i);
  388.             MakeStepL(inplus,leplus,rie);
  389.             MakeStep (igplus,seplus,rie);
  390.         }
  392.         if (seqscr.StopScore(i,Minus) != BadScore) {
  393.             PushState(igminus,Minus,i);
  394.             MakeStep (seplus,igminus,rine);
  395.             MakeStep (seminus,igminus,rine);
  396.             MakeStep (leplus,igminus,rine);
  397.             MakeStep (feminus,igminus,rine);
  398.         }
  399.     }
  401.     PushState (inplus,Plus,-1);
  402.     MakeStepLR(ieplus,inplus,rine,1);
  403.     MakeStepLR(feplus,inplus,rine,1);
  405.     PushState (ieminus,Minus,-1);
  406.     MakeStepLR(inminus,ieminus,rie);
  408.     PushState(leminus,Minus,-1);
  409.     MakeStepR(igminus,leminus,rie);
  411.     PushState (ieplus,Plus,-1);
  412.     MakeStepLR(inplus,ieplus,rie);
  414.     PushState(feplus,Plus,-1);
  415.     MakeStepR(igplus,feplus,rie);
  417.     PushState (inminus,Minus,-1);
  418.     MakeStepLR(leminus,inminus,rine,0);
  419.     MakeStepLR(ieminus,inminus,rine,0);
  421.     PushState(igplus,Plus,-1);
  422.     MakeStep (seplus,igplus,rine);
  423.     MakeStep (seminus,igplus,rine);
  424.     MakeStep (leplus,igplus,rine);
  425.     MakeStep (feminus,igplus,rine);
  427.     PushState(feminus,Minus,-1);
  428.     MakeStepL(inminus,feminus,rie);
  430.     PushState(seminus,Minus,-1);
  431.     MakeStep (igminus,seminus,rie);
  433.     PushState(leplus,Plus,-1);
  434.     MakeStepL(inplus,leplus,rie);
  436.     PushState(seplus,Plus,-1);
  437.     MakeStep(igplus,seplus,rie);
  439.     PushState(igminus,Minus,-1);
  440.     MakeStep (seplus,igminus,rine);
  441.     MakeStep (seminus,igminus,rine);
  442.     MakeStep (leplus,igminus,rine);
  443.     MakeStep (feminus,igminus,rine);
  445.     igplus.back().UpdateScore(AddProbabilities(igplus.back().Score(),igminus.back().Score()));
  447.     const HMM_State* p = &igplus.back();
  448.     if (feminus.back().Score() > p->Score()) {
  449.         p = &feminus.back();
  450.     }
  451.     if (leplus.back().Score() > p->Score()) {
  452.         p = &leplus.back();
  453.     }
  454.     if (seplus.back().Score() > p->Score()) {
  455.         p = &seplus.back();
  456.     }
  457.     if (seminus.back().Score() > p->Score()) {
  458.         p = &seminus.back();
  459.     }
  460.     for (int ph = 0;  ph < 3;  ++ph) {
  461.         if (inplus[ph].back().Score() > p->Score()) {
  462.             p = &inplus[ph].back();
  463.         }
  464.         if (inminus[ph].back().Score() > p->Score()) {
  465.             p = &inminus[ph].back();
  466.         }
  467.         if (ieplus[ph].back().Score() > p->Score()) {
  468.             p = &ieplus[ph].back();
  469.         }
  470.         if (ieminus[ph].back().Score() > p->Score()) {
  471.             p = &ieminus[ph].back();
  472.         }
  473.         if (leminus[ph].back().Score() > p->Score()) {
  474.             p = &leminus[ph].back();
  475.         }
  476.         if (feplus[ph].back().Score() > p->Score()) {
  477.             p = &feplus[ph].back();
  478.         }
  479.     }
  481.     path = p;
  482. }
  485. template<class T> void Out(T t, int w, CNcbiOstream& to = cout)
  486. {
  487.     to.width(w);
  488.     to.setf(IOS_BASE::right,IOS_BASE::adjustfield);
  489.     to << t;
  490. }
  493. void Out(double t, int w, CNcbiOstream& to = cout, int prec = 1)
  494. {
  495.     to.width(w);
  496.     to.setf(IOS_BASE::right,IOS_BASE::adjustfield);
  497.     to.setf(IOS_BASE::fixed,IOS_BASE::floatfield);
  498.     to.precision(prec);
  500.     if (t > 1000000000) {
  501.         to << "+Inf";
  502.     } else if (t < -1000000000) {
  503.         to << "-Inf";
  504.     } else {
  505.         to << t;
  506.     }
  507. }
  510. inline char toACGT(int c)
  511. {
  512.     switch (c) {
  513.     case nA: 
  514.         return 'A';
  515.     case nC: 
  516.         return 'C';
  517.     case nG: 
  518.         return 'G';
  519.     case nT: 
  520.         return 'T';
  521.     case nN: 
  522.         return 'N';
  523.     }
  524. }
  527. int Parse::PrintGenes(CNcbiOstream& to, CNcbiOstream& toprot, bool complete) const
  528. {
  529.     enum {DNA_Align = 1, Prot_Align = 2 };
  531.     list<Gene> genes = GetGenes();
  533.     int right = seqscr.SeqMap(seqscr.SeqLen() - 1,true);
  534.     if (complete  &&  !genes.empty()  &&  !genes.back().RightComplete()) {
  535.         int partial_start = genes.back().front().Start();
  536.         genes.pop_back();
  538.         if ( !genes.empty() ) // end of the last complete gene
  539.         {
  540.             right = genes.back().back().Stop();
  541.         } else {
  542.             if (partial_start > seqscr.SeqMap(0,true) + 1000) {
  543.                 right = partial_start - 100;
  544.             } else {
  545.                 return -1;   // calling program MUST be aware of this!!!!!!!
  546.             }
  547.         }
  548.     }
  550.     to << "n" << seqscr.Contig() << '_' << seqscr.SeqMap(0,true) << '_' << right << "nn";
  551.     to << genes.size() << " genes foundn";
  553.     set<int> chain_id, prot_id;
  554.     for (list<Gene>::iterator it = genes.begin();  it != genes.end();  ++it) {
  555.         const Gene& gene = *it;
  556.         for (int i = 0;  i < gene.size();  ++i) {
  557.             ITERATE(set<int>, iter, gene[i].ChainID()) {
  558.                 chain_id.insert(*iter);
  559.             }
  561.             ITERATE(set<int>, iter, gene[i].ProtID()) {
  562.                 prot_id.insert(*iter);
  563.             }
  565.             //chain_id.insert(gene[i].ChainID().begin(), gene[i].ChainID().end());
  566.             //prot_id.insert(gene[i].ProtID().begin(),  gene[i].ProtID().end());
  567.         }
  568.     }
  569.     to << (chain_id.size() + prot_id.size()) << " alignments usedn";
  571.     int num = 0;
  572.     for (list<Gene>::iterator it = genes.begin();  it != genes.end();  ++it) {
  573.         ++num;
  574.         to << "nn";
  575.         Out(" ",11,to);
  576.         Out("Start",10,to);
  577.         Out("Stop",10,to);
  578.         Out("Length",7,to);
  579.         Out("Align",10,to);
  580.         Out("Prot",12,to);
  581.         Out("FShift",7,to);
  582.         to << endl;
  584.         const Gene& gene = *it;
  585.         int gene_start = -1;
  586.         int gene_stop;
  587.         int align_type = 0;
  589.         for (int i = 0;  i < gene.size();  ++i) {
  590.             const ExonData& exon = gene[i];
  591.             const TFrameShifts& fshifts = exon.ExonFrameShifts();
  592.             int estart = exon.Start();
  593.             if (gene_start < 0  &&  exon.Type() == ExonData::Cds) {
  594.                 gene_start = estart;
  595.             }
  596.             if ( !exon.ChainID().empty() ) {
  597.                 align_type = align_type|DNA_Align;
  598.             }
  599.             if ( !exon.ProtID().empty() ) {
  600.                 align_type = align_type|Prot_Align;
  601.             }
  603.             for (int k = 0;  k < fshifts.size();  ++k) {
  604.                 int estop = fshifts[k].Loc() - 1;
  606.                 Out(num,4,to);
  607.                 Out((exon.Type() == ExonData::Cds) ? "CDS" : "UTR",5,to); 
  608.                 Out((gene.Strand() == Plus) ? '+' : '-',2,to);
  609.                 Out(estart,10,to);
  610.                 Out(estop,10,to);
  611.                 Out(estop - estart + 1,7,to);
  612.                 if (exon.ChainID().empty()) {
  613.                     Out("-",10,to);
  614.                 } else {
  615.                     Out(*exon.ChainID().begin(),10,to);
  616.                 }
  617.                 if (exon.ProtID().empty()) {
  618.                     Out("-",12,to);
  619.                 } else {
  620.                     Out(*exon.ProtID().begin(),12,to);
  621.                 }
  622.                 if (fshifts[k].IsInsertion()) {
  623.                     Out("+",7,to);
  624.                     estart = estop;
  625.                 } else {
  626.                     Out("-",7,to);
  627.                     estart = estop + 2;
  628.                 }
  629.                 to << endl;
  630.             }
  631.             int estop = exon.Stop();
  632.             Out(num,4,to);
  633.             Out((exon.Type() == ExonData::Cds) ? "CDS" : "UTR",5,to); 
  634.             Out((gene.Strand() == Plus) ? '+' : '-',2,to);
  635.             Out(estart,10,to);
  636.             Out(estop,10,to);
  637.             Out(estop - estart + 1,7,to);
  638.             if (exon.ChainID().empty()) {
  639.                 Out("-",10,to);
  640.             } else {
  641.                 Out(*exon.ChainID().begin(),10,to);
  642.             }
  643.             if (exon.ProtID().empty()) {
  644.                 Out("-",12,to);
  645.             } else {
  646.                 Out(*exon.ProtID().begin(),12,to);
  647.             }
  648.             Out("*",7,to);
  649.             to << endl;
  650.             if (exon.Type() == ExonData::Cds) {
  651.                 gene_stop = estop;
  652.             }
  653.         }   
  655.         toprot << '>' << seqscr.Contig() << "_" 
  656.             << gene_start << "_" 
  657.             << gene_stop << "_"
  658.             << ((gene.Strand() == Plus) ? "plus" : "minus") << "_"
  659.             << align_type << endl;
  661.         const IVec& cds = gene.CDS();
  662.         int i;
  663.         for (i = 0;  i < cds.size() - 2;  i +=3) {
  664.             toprot << aa_table[cds[i]*25 + cds[i + 1]*5 + cds[i + 2]];
  665.             if ((i + 3)%150 == 0) {
  666.                 toprot << endl;
  667.             }
  668.         }
  669.         if (i%150 != 0) {
  670.             toprot << endl;
  671.         }
  673.         //      for (i = 0;  i < cds.size();  ++i)
  674.         //      {
  675.         //          toprot << toACGT(cds[i]);
  676.         //          if ((i + 3)%50 == 0) toprot << endl;
  677.         //      }
  678.         //      if (i%50 != 0) toprot << endl;
  679.     }
  681.     return right;
  682. }
  685. /*
  686. int Parse::PrintGenes(CNcbiOstream& to, CNcbiOstream& toprot, bool complete) const
  687. {
  688.     const char* aa_table = "KNKNXTTTTTRSRSXIIMIXXXXXXQHQHXPPPPPRRRRRLLLLLXXXXXEDEDXAAAAAGGGGGVVVVVXXXXX * Y*YXSSSSS * CWCXLFLFXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX";
  690.     const HMM_State* plast = Path();
  691.     if (complete  &&  dynamic_cast<const Intron*>(plast)) {
  692.         while (!dynamic_cast<const Intergenic*>(plast)) {
  693.             plast = plast->LeftState();
  694.         }
  696.         if ( !plast->LeftState() ) {
  697.             return plast->Start() + seqscr.SeqMap(0,true) - 1;
  698.         } else {
  699.             plast = plast->LeftState();
  700.         }
  701.     }
  703.     const CClusterSet& algns = seqscr.Alignments();
  705.     int right = seqscr.SeqMap(plast->Stop(),false);  // either end of chunk or and of last complete gene and its alignments
  706.     for (CClusterSet::ConstIt it = algns.begin();  it != algns.end();  ++it) {
  707.         if (it->Limits().Include(right)) {
  708.             right = it->Limits().second;
  709.         }
  710.     }
  712.     to << "n" << seqscr.Contig() << '_' << seqscr.SeqMap(0,true) << '_' << right << "nn";
  714.     vector<const Exon*> states;
  715.     const Exon* pe;
  716.     int gennum = 0;
  717.     for (const HMM_State* p = plast;  p != 0;  p = p->LeftState()) {
  718.         if (pe = dynamic_cast<const Exon*>(p)) {
  719.             states.push_back(pe);
  720.         }
  722.         if ( !states.empty() ) {
  723.             if ( !p->LeftState()  ||  dynamic_cast<const Intergenic*>(p) ) {
  724.                 ++gennum;
  725.             }
  726.         }
  727.     }
  728.     reverse(states.begin(),states.end());
  730.     int algnum = 0, algused = 0;
  731.     IPair chunk(seqscr.SeqMap(0,true),seqscr.SeqMap(seqscr.SeqLen() - 1,false));
  732.     for (CClusterSet::ConstIt it = algns.begin();  it != algns.end();  ++it) {
  733.         IPair lim = it->Limits();
  734.         if (lim.Intersect(chunk)) {
  735.             algnum += it->size();
  736.         }
  737.         for (int i = 0;  i < states.size();  ++i) {
  738.             int a = seqscr.SeqMap(states[i]->Start(),true);
  739.             int b = seqscr.SeqMap(states[i]->Stop(),false);
  740.             IPair exon(a,b);
  741.             if (exon.Intersect(lim)) {
  742.                 algused += it->size();
  743.                 break;
  744.             }
  745.         }
  746.     }
  747.     IVec ids(states.size(),0), prot_ids(states.size(),0);
  748.     for (int i = 0;  i < states.size();  ++i) {
  749.         int a = seqscr.SeqMap(states[i]->Start(),true);
  750.         int b = seqscr.SeqMap(states[i]->Stop(),false);
  751.         IPair exon(a,b);
  752.         for (CClusterSet::ConstIt cls_it = algns.begin();  cls_it != algns.end();  ++cls_it) {
  753.             for (CCluster::ConstIt it = cls_it->begin();  it != cls_it->end();  ++it)
  754.             {
  755.                 const AlignVec& al = *it;
  756.                 for (int j = 0;  j < al.size();  ++j) {
  757.                     if (exon.Intersect(al[j])) {
  758.                         if (al.Type() == AlignVec::Prot) {
  759.                             prot_ids[i] = al.ID();
  760.                         } else {
  761.                             ids[i] = al.ID();
  762.                             if (exon.first != al[j].first  ||  exon.second != al[j].second) {
  763.                                 ids[i] = -ids[i];
  764.                             }
  765.                         }
  766.                     }
  767.                 }
  768.             }
  769.         }
  770.     }
  773.     to << gennum << " genes foundn";
  774.     to << algused << " alignments used out of " << algnum << endl;
  775.     if (gennum == 0) {
  776.         return seqscr.SeqMap(plast->Stop(),false);
  777.     }
  779.     IVec cds;
  780.     int num = 0;
  781.     double gene_fscore, left_score;
  782.     int gene_start, gene_stop, align_type;
  783.     enum {DNA_Align = 1, Prot_Align = 2 };
  784.     for (int i = 0;  i < states.size();  ++i) {
  785.         pe = states[i];
  786.         int shift = 0;
  788.         if (i == 0  ||  
  789.             dynamic_cast<const SingleExon*>(pe)  ||  
  790.             (pe->isPlus()  &&  dynamic_cast<const FirstExon*>(pe))  ||  
  791.             (pe->isMinus()  &&  dynamic_cast<const LastExon*>(pe)))
  792.         {
  793.             gene_start = seqscr.SeqMap(pe->Start(),true);
  794.             ++num;
  795.             to << "nn";
  796.             Out(" ",19,to);
  797.             Out("Start",10,to);
  798.             Out("Stop",10,to);
  799.             Out("Length",7,to);
  800.             Out("Align",10,to);
  801.             Out("Prot",12,to);
  802.             Out("FShift",7,to);
  803.             to << endl;
  805.             left_score = pe->LeftState()->Score();
  807.             cds.clear();
  808.             const Intron* pi = dynamic_cast<const Intron*>(pe->LeftState());
  809.             if (pi != 0) {
  810.                 int ph = pi->Phase();
  811.                 if (pe->isPlus()) {
  812.                     shift = (3 - ph)%3;
  813.                 } else {
  814.                     shift = ph;
  815.                 }
  816.             }
  818.             align_type = 0;
  819.         }
  821.         //If first / last base is insertion, this information is lost
  822.         vector<IPair> subexons;
  823.         int a = seqscr.SeqMap(pe->Start(),true);
  824.         int b = a;
  825.         for (int k = pe->Start() + 1;  k <= pe->Stop();  ++k) {
  826.             int c = seqscr.SeqMap(k,false);
  827.             if (c == b) {
  828.                 subexons.push_back(IPair(a,b));
  829.                 a = b;
  830.             }
  831.             if (c == b + 2) {
  832.                 subexons.push_back(IPair(a,b));
  833.                 a = c;
  834.             }
  835.             b = c;
  836.         }
  837.         subexons.push_back(IPair(a,b));
  839.         for (int k = 0;  k < subexons.size();  ++k) {
  840.             int a = subexons[k].first;
  841.             int b = subexons[k].second;
  842.             Out(num,4,to);
  843.             Out(pe->GetStateName(),13,to);
  844.             Out(pe->isPlus() ? '+' : '-',2,to);
  846.             Out(a,10,to);
  847.             Out(b,10,to);
  848.             Out(b - a+1,7,to);
  850.             if (ids[i]) {
  851.                 align_type = align_type|DNA_Align;
  852.                 Out(abs(ids[i]),9,to);
  853.                 to << ((ids[i] > 0) ? ' ' : '*');
  854.             } else {
  855.                 Out("-",9,to);
  856.                 to << ' ';
  857.             }
  858.             if (prot_ids[i]) {
  859.                 align_type = align_type|Prot_Align;
  860.                 Out(prot_ids[i],12,to);
  861.             } else {
  862.                 Out("-",12,to);
  863.             }
  865.             if (k != 0) {
  866.                 if (a == subexons[k - 1].second) {
  867.                     Out("+",7,to);
  868.                 } else if (a == subexons[k - 1].second + 2) {
  869.                     Out("-",7,to);
  870.                 } else {
  871.                     Out("*",7,to);
  872.                 }
  873.             } else {
  874.                 Out("*",7,to);
  875.             }
  877.             to << endl;
  878.         }
  880.         bool last = false;
  882.         if (i == states.size() - 1  ||  
  883.             dynamic_cast<const SingleExon*>(pe)  ||  
  884.             (pe->isPlus()  &&  dynamic_cast<const LastExon*>(pe))  ||  
  885.             (pe->isMinus()  &&  dynamic_cast<const FirstExon*>(pe)))
  886.         {
  887.             last = true;
  888.             gene_stop = seqscr.SeqMap(pe->Stop(),false);
  889.         }
  891.         copy(seqscr.SeqPtr(pe->Start() + shift,Plus),
  892.              seqscr.SeqPtr(pe->Stop() + 1,Plus),back_inserter(cds));
  894.         if (last) {
  895.             to << endl;
  897.             cds.resize(cds.size() - cds.size()%3);
  899.             if (pe->isMinus()) {
  900.                 for (int i = 0;  i < cds.size();  ++i) cds[i] = toMinus[cds[i]];
  901.                 reverse(cds.begin(),cds.end());
  902.             }
  904.             toprot << '>' << seqscr.Contig() << "_" 
  905.                 << gene_start << "_" 
  906.                 << gene_stop << "_"
  907.                 << (pe->isPlus() ? "plus" : "minus") << "_"
  908.                 << align_type << endl;
  910.             int i;
  911.             for (i = 0;  i < cds.size() - 2;  i +=3) {
  912.                 toprot << aa_table[cds[i]*25 + cds[i + 1]*5 + cds[i + 2]];
  913.                 if ((i + 3)%150 == 0) {
  914.                     toprot << endl;
  915.                 }
  916.             }
  917.             if (i%150 != 0) {
  918.                 toprot << endl;
  919.             }
  920.         }
  921.     }
  923.     return right;
  924. }
  925. */
  928. list<Gene> Parse::GetGenes() const
  929. {
  930.     list<Gene> genes;
  932.     const CClusterSet& cls = seqscr.Alignments();
  934.     if (dynamic_cast<const Intron*>(Path())) {
  935.         genes.push_front(Gene(Path()->Strand(),false,false));
  936.     }
  938.     for (const HMM_State* p = Path();  p != 0;  p = p->LeftState()) {
  939.         if (const Exon* pe = dynamic_cast<const Exon*>(p)) {
  940.             bool stopcodon = false;
  941.             bool startcodon = false;
  943.             if (dynamic_cast<const SingleExon*>(pe)) {
  944.                 genes.push_front(Gene(pe->Strand(),true,true));
  945.                 startcodon = true;
  946.                 stopcodon = true;
  947.             } else if (dynamic_cast<const LastExon*>(pe)) {
  948.                 if (pe->isPlus()) {
  949.                     genes.push_front(Gene(pe->Strand(),false,true));
  950.                 } else {
  951.                     genes.front().leftend = true;
  952.                 }
  953.                 stopcodon = true;
  954.             } else if (dynamic_cast<const FirstExon*>(pe)) {
  955.                 if (pe->isMinus()) {
  956.                     genes.push_front(Gene(pe->Strand(),false,true));
  957.                 } else {
  958.                     genes.front().leftend = true;
  959.                 }
  960.                 startcodon = true;
  961.             }
  963.             Gene& curgen = genes.front();
  965.             int estart = pe->Start();
  966.             int estop = pe->Stop();
  968.             if (pe->isPlus()) {
  969.                 int ph = (pe->Phase() - estop + estart)%3;
  970.                 if (ph < 0) {
  971.                     ph += 3;
  972.                 }
  973.                 curgen.cds_shift = ph;
  974.             } else if (curgen.empty()) {
  975.                 curgen.cds_shift = pe->Phase();
  976.             }
  978.             if (pe->isMinus()  &&  startcodon) {
  979.                 curgen.cds.push_back(nT);
  980.                 curgen.cds.push_back(nA);
  981.                 curgen.cds.push_back(nC);
  982.             }
  983.             for (int i = estop;  i >= estart;  --i) {
  984.                 curgen.cds.push_back(*seqscr.SeqPtr(i,Plus));
  985.             }
  986.             if (pe->isPlus()  &&  startcodon) {
  987.                 curgen.cds.push_back(nG);
  988.                 curgen.cds.push_back(nT);
  989.                 curgen.cds.push_back(nA);
  990.             }
  992.             int a = seqscr.SeqMap(estart,true);
  993.             bool a_insertion = (estart != seqscr.RevSeqMap(a,true));
  994.             int b = seqscr.SeqMap(estop,false);
  995.             bool b_insertion = (estop != seqscr.RevSeqMap(b,true));
  997.             pair<CClusterSet::ConstIt,CClusterSet::ConstIt> lim = cls.equal_range(CCluster(a,b));
  998.             CClusterSet::ConstIt first = lim.first;
  999.             CClusterSet::ConstIt last = lim.second;
  1000.             if (lim.first != lim.second) {
  1001.                 --last;
  1002.             }
  1004.             int lastid = 0;
  1005.             int margin = numeric_limits<int>::min();   // difference between end of exon and end of alignment exon
  1006.             bool rightexon = false;
  1007.             if (((startcodon  &&  pe->isMinus())  ||  (stopcodon  &&  pe->isPlus()))  &&  lim.first != lim.second)   // rightside UTR
  1008.             {
  1009.                 for (CCluster::ConstIt it = last->begin();  it != last->end();  ++it) {
  1010.                     const AlignVec& algn = *it;
  1011.                     if (algn.Type() == AlignVec::Prot) {
  1012.                         continue;
  1013.                     }
  1015.                     for (int i = algn.size() - 1;  i >= 0;  --i) {
  1016.                         if (b < algn[i].first) {
  1017.                             curgen.push_back(ExonData(algn[i].first,algn[i].second,ExonData::Utr));
  1018.                             curgen.back().chain_id.insert(algn.ID()); 
  1019.                             rightexon = true;
  1020.                         } else if (algn[i].second >= a  &&  algn[i].second - b > margin) {
  1021.                             margin = algn[i].second - b;
  1022.                             lastid = algn.ID();
  1023.                         }
  1024.                     }
  1025.                 }
  1026.             }
  1029.             int remain = 0;
  1030.             if ((pe->isMinus()  &&  startcodon)  ||  (pe->isPlus()  &&  stopcodon))   // start / stop position correction
  1031.             {
  1032.                 if (margin >= 3) {
  1033.                     b += 3;
  1034.                     margin -= 3;
  1035.                 } else if (margin >= 0  &&  rightexon) {
  1036.                     b += margin;
  1037.                     remain = 3 - margin;
  1038.                     margin = 0;
  1039.                 } else {
  1040.                     b += 3;
  1041.                     margin = 0;
  1042.                 }
  1043.             }
  1044.             if (margin > 0) {
  1045.                 curgen.push_back(ExonData(b + 1,b + margin,ExonData::Utr)); 
  1046.                 curgen.back().chain_id.insert(lastid); 
  1047.             } else if (remain > 0) {
  1048.                 ExonData& e = curgen.back();
  1049.                 int aa = e.start;
  1050.                 e.start += remain;
  1051.                 int bb = aa + remain - 1;
  1052.                 ExonData er(aa,bb,ExonData::Cds);
  1053.                 er.chain_id = e.chain_id;
  1054.                 curgen.push_back(er);
  1055.             }
  1057.             curgen.push_back(ExonData(a,b,ExonData::Cds));
  1058.             ExonData& exon = curgen.back();
  1059.             margin = numeric_limits<int>::min();
  1060.             bool leftexon = false;
  1061.             lastid = 0;
  1062.             for (CClusterSet::ConstIt cls_it = lim.first;  cls_it != lim.second;  ++cls_it) {
  1063.                 for (CCluster::ConstIt it = cls_it->begin();  it != cls_it->end();  ++it)
  1064.                 {
  1065.                     const AlignVec& algn = *it;
  1067.                     for (int i = 0;  i < algn.size();  ++i) {
  1068.                         if (exon.stop < algn[i].first) {
  1069.                             break;
  1070.                         } else if (algn[i].second < exon.start) {
  1071.                             leftexon = true;
  1072.                             continue;
  1073.                         } else if (algn.Type() == AlignVec::Prot) {
  1074.                             exon.prot_id.insert(algn.ID()); 
  1075.                         } else {
  1076.                             exon.chain_id.insert(algn.ID());
  1077.                             if (a - algn[i].first > margin) {
  1078.                                 margin = a - algn[i].first;
  1079.                                 lastid = algn.ID();
  1080.                             }
  1081.                         }
  1082.                     }
  1083.                 }
  1084.             }
  1086.             const TFrameShifts& fs = seqscr.SeqFrameShifts();
  1087.             for (int i = 0;  i < fs.size();  ++i) {
  1088.                 if ((fs[i].Loc() == a  &&  a_insertion)  ||  // exon starts from insertion
  1089.                     (fs[i].Loc() == b + 1  &&  b_insertion)  ||  // exon ends by insertion
  1090.                     (fs[i].Loc() <= b  &&  fs[i].Loc() > a)) {
  1091.                     exon.fshifts.push_back(fs[i]);
  1092.                 }
  1093.             }
  1095.             remain = 0;
  1096.             if ((pe->isPlus()  &&  startcodon)  ||  (pe->isMinus()  &&  stopcodon))   // start / stop position correction
  1097.             {
  1098.                 if (margin >= 3) {
  1099.                     exon.start -= 3;
  1100.                     margin -= 3;
  1101.                 } else if (margin >= 0  &&  leftexon) {
  1102.                     exon.start -= margin;
  1103.                     remain = 3 - margin;
  1104.                     margin = 0;
  1105.                 } else {
  1106.                     exon.start -= 3;
  1107.                     margin = 0;
  1108.                 }
  1109.             }
  1111.             if (((startcodon  &&  pe->isPlus())  ||  (stopcodon  &&  pe->isMinus()))  &&  lim.first != lim.second)   // lefttside UTR
  1112.             {
  1113.                 int estart = exon.start;
  1115.                 if (margin > 0) {
  1116.                     curgen.push_back(ExonData(estart - margin,estart - 1,ExonData::Utr)); 
  1117.                     curgen.back().chain_id.insert(lastid); 
  1118.                 }
  1120.                 for (CCluster::ConstIt it = first->begin();  it != first->end();  ++it) {
  1121.                     const AlignVec& algn = *it;
  1122.                     if (algn.Type() == AlignVec::Prot) {
  1123.                         continue;
  1124.                     }
  1126.                     for (int i = algn.size() - 1;  i >= 0;  --i) {
  1127.                         if (estart > algn[i].second) {
  1128.                             int aa = algn[i].first;
  1129.                             int bb = algn[i].second;
  1130.                             if (remain > 0) {
  1131.                                 curgen.push_back(ExonData(bb - remain + 1,bb,ExonData::Cds));
  1132.                                 curgen.back().chain_id.insert(algn.ID()); 
  1133.                                 bb -= remain;
  1134.                                 remain = 0;
  1135.                             }
  1136.                             curgen.push_back(ExonData(aa,bb,ExonData::Utr)); 
  1137.                             curgen.back().chain_id.insert(algn.ID()); 
  1138.                         }
  1139.                     }
  1140.                 }
  1141.             }
  1142.         }
  1143.     }
  1145.     for (list<Gene>::iterator it = genes.begin();  it != genes.end();  ++it) {
  1146.         Gene& g = *it;
  1147.         reverse(g.begin(),g.end());
  1149.         if (g.Strand() == Plus) {
  1150.             reverse(g.cds.begin(),g.cds.end());
  1151.         } else {
  1152.             for (int i = 0;  i < g.cds.size();  ++i) {
  1153.                 g.cds[i] = toMinus[g.cds[i]];
  1154.             }
  1155.         }
  1157.     }
  1159.     return genes;
  1160. }
  1163. void Parse::PrintInfo() const
  1164. {
  1165.     vector<const HMM_State*> states;
  1166.     for (const HMM_State* p = Path();  p != 0;  p = p->LeftState()) states.push_back(p);
  1167.     reverse(states.begin(),states.end());
  1169.     Out(" ",15);
  1170.     Out("Start",8);
  1171.     Out("Stop",8);
  1172.     Out("Score",8);
  1173.     Out("BrScr",8);
  1174.     Out("LnScr",8);
  1175.     Out("RgnScr",8);
  1176.     Out("TrmScr",8);
  1177.     Out("AccScr",8);
  1178.     cout << endl;
  1180.     for (int i = 0;  i < states.size();  ++i) {
  1181.         const HMM_State* p = states[i];
  1182.         if (dynamic_cast<const Intergenic*>(p)) {
  1183.             cout << endl;
  1184.         }
  1186.         Out(p->GetStateName(),13);
  1187.         Out(p->isPlus() ? '+' : '-',2);
  1188.         int a = seqscr.SeqMap(p->Start(),true);
  1189.         int b = seqscr.SeqMap(p->Stop(),false);
  1190.         Out(a,8);
  1191.         Out(b,8);
  1192.         StateScores sc = p->GetStateScores();
  1193.         Out(sc.score,8);
  1194.         Out(sc.branch,8);
  1195.         Out(sc.length,8);
  1196.         Out(sc.region,8);
  1197.         Out(sc.term,8);
  1198.         Out(p->Score(),8);
  1199.         cout << endl;
  1200.     }
  1201. }
  1204. void LogicalCheck(const HMM_State& st, const SeqScores& ss)
  1205. {
  1206.     typedef const Intergenic* Ig;
  1207.     typedef const Intron* I;
  1208.     typedef const FirstExon* FE;
  1209.     typedef const InternalExon* IE;
  1210.     typedef const LastExon* LE;
  1211.     typedef const SingleExon* SE;
  1212.     typedef const Exon* E;
  1214.     bool out;
  1215.     int phase,strand;
  1216.     vector< pair<int,int> > exons;
  1217.     for (const HMM_State* state = &st;  state != 0;  state = state->LeftState()) {
  1218.         const HMM_State* left = state->LeftState();
  1219.         if (left == 0) {
  1220.             if ( !dynamic_cast<Ig>(state)  &&  !dynamic_cast<I>(state) ) {
  1221.                 NCBI_THROW(CGnomonException, eGenericError,
  1222.                            "GNOMON: Logical error 1");
  1223.             }
  1225.             out = true;
  1226.         } else if (dynamic_cast<Ig>(state)) {
  1227.             if ( !dynamic_cast<SE>(left)  &&  
  1228.                  !(dynamic_cast<LE>(left)  &&  left->isPlus())  &&  
  1229.                  !(dynamic_cast<FE>(left)  &&  left->isMinus()) ) {
  1230.                 NCBI_THROW(CGnomonException, eGenericError,
  1231.                            "GNOMON: Logical error 2");
  1232.             }
  1234.             out = true;
  1235.         } else if (I ps = dynamic_cast<I>(state)) {
  1236.             E p;
  1237.             int ph = -1;
  1238.             if (((p = dynamic_cast<FE>(left))  &&  left->isPlus())  ||  
  1239.                 ((p = dynamic_cast<LE>(left))  &&  left->isMinus())  ||  
  1240.                 (p = dynamic_cast<IE>(left))) {
  1241.                 ph = p->Phase();
  1242.             }
  1244.             if ((ph < 0)  ||  (ps->Strand() != left->Strand())  ||  
  1245.                 (ps->isPlus()  &&  (ph + 1)%3 != ps->Phase())  ||  
  1246.                 (ps->isMinus()  &&  ph != ps->Phase())) {
  1247.                 NCBI_THROW(CGnomonException, eGenericError,
  1248.                            "GNOMON: Logical error 3");
  1249.             }
  1251.             out = false;
  1252.         } else if (SE ps = dynamic_cast<SE>(state)) {
  1253.             if ( !(dynamic_cast<Ig>(left)  &&  ps->Strand() == left->Strand()) ) {
  1254.                 NCBI_THROW(CGnomonException, eGenericError,
  1255.                            "GNOMON: Logical error 4");
  1256.             }
  1258.             int a = ps->Start(), b = ps->Stop();
  1259.             if ((b - a+1)%3  ||  
  1260.                 (ps->isPlus()  &&  (!ss.isStart(a,Plus)  ||  !ss.isStop(b + 1,Plus)))  ||  
  1261.                 (ps->isMinus()  &&  (!ss.isStart(b,Minus)  ||  !ss.isStop(a - 1,Minus))))  {
  1262.                 NCBI_THROW(CGnomonException, eGenericError,
  1263.                            "GNOMON: Logical error 5");
  1264.             }
  1265.         } else if (FE ps = dynamic_cast<FE>(state)) {
  1266.             I p;
  1267.             if ((ps->Strand() != left->Strand())  ||  
  1268.                 (!(dynamic_cast<Ig>(left)  &&  ps->isPlus())  &&  
  1269.                  !((p = dynamic_cast<I>(left))  &&  ps->isMinus()))) {
  1270.                 NCBI_THROW(CGnomonException, eGenericError,
  1271.                            "GNOMON: Logical error 6");
  1272.             }
  1274.             int a = ps->Start(), b = ps->Stop();
  1275.             if (ps->isPlus()) {
  1276.                 if (ps->Phase() != (b - a)%3) {
  1277.                     NCBI_THROW(CGnomonException, eGenericError,
  1278.                                "GNOMON: Logical error 7");
  1279.                 }
  1280.             } else {
  1281.                 if (p->Phase() != (b + 1-a)%3) {
  1282.                     NCBI_THROW(CGnomonException, eGenericError,
  1283.                                "GNOMON: Logical error 8");
  1284.                 }
  1285.             }
  1287.             //          if ((ps->isPlus()  &&  (!ss.isStart(a,Plus)  ||  !ss.isGT(b + 1,Plus)))  ||  
  1288.             //             (ps->isMinus()  &&  (!ss.isStart(b,Minus)  ||  !ss.isGT(a - 1,Minus)))) 
  1289.             //          { cout << "Logical error9n"; exit(1); }
  1290.         } else if (IE ps = dynamic_cast<IE>(state)) {
  1291.             I p;
  1292.             if ((ps->Strand() != left->Strand())  ||  
  1293.                 !(p = dynamic_cast<I>(left))) {
  1294.                 NCBI_THROW(CGnomonException, eGenericError,
  1295.                            "GNOMON: Logical error 10");
  1296.             }
  1298.             int a = ps->Start(), b = ps->Stop();
  1299.             if (ps->isPlus()) {
  1300.                 if (ps->Phase() != (p->Phase() + b-a)%3) {
  1301.                     NCBI_THROW(CGnomonException, eGenericError,
  1302.                                "GNOMON: Logical error 11");
  1303.                 }
  1304.             } else {
  1305.                 if (p->Phase() != (ps->Phase() + b+1 - a)%3) {
  1306.                     NCBI_THROW(CGnomonException, eGenericError,
  1307.                                "GNOMON: Logical error 12");
  1308.                 }
  1309.             }
  1311.             //          if ((ps->isPlus()  &&  (!ss.isAG(a - 1,Plus)  ||  !ss.isGT(b + 1,Plus)))  ||  
  1312.             //             (ps->isMinus()  &&  (!ss.isAG(b + 1,Minus)  ||  !ss.isGT(a - 1,Minus)))) 
  1313.             //          { cout << "Logical error13n"; exit(1); }
  1314.         } else if (LE ps = dynamic_cast<LE>(state)) {
  1315.             I p;
  1316.             if ((ps->Strand() != left->Strand())  ||  
  1317.                 (!(dynamic_cast<Ig>(left)  &&  ps->isMinus())  &&  
  1318.                  !((p = dynamic_cast<I>(left))  &&  ps->isPlus()))) {
  1319.                 NCBI_THROW(CGnomonException, eGenericError,
  1320.                            "GNOMON: Logical error 14");
  1321.             }
  1323.             int a = ps->Start(), b = ps->Stop();
  1324.             if (ps->isPlus()) {
  1325.                 if ((p->Phase() + b - a) % 3 != 2) {
  1326.                     NCBI_THROW(CGnomonException, eGenericError,
  1327.                                "GNOMON: Logical error 15");
  1328.                 }
  1329.             } else {
  1330.                 if ((ps->Phase() + b-a)%3 != 2) {
  1331.                     NCBI_THROW(CGnomonException, eGenericError,
  1332.                                "GNOMON: Logical error 16");
  1333.                 }
  1334.             }
  1336.             //          if ((ps->isPlus()  &&  (!ss.isAG(a - 1,Plus)  ||  !ss.isStop(b + 1,Plus)))  ||  
  1337.             //             (ps->isMinus()  &&  (!ss.isStop(a - 1,Minus)  ||  !ss.isAG(b + 1,Minus)))) 
  1338.             //          { cout << "Logical error17n"; exit(1); }
  1339.         }
  1341.         if (E ps = dynamic_cast<E>(state)) {
  1342.             phase = ps->Phase();
  1343.             strand = ps->Strand();
  1344.             exons.push_back(make_pair(ps->Start(),ps->Stop()));
  1345.             out = false;
  1346.         }
  1348.         if (out  &&  !exons.empty()) {
  1349.             IVec seq;
  1351.             if (strand == Plus) {
  1352.                 for (int i = exons.size() - 1;  i >= 0;  --i) {
  1353.                     int a = exons[i].first;
  1354.                     int b = exons[i].second;
  1355.                     seq.insert(seq.end(),ss.SeqPtr(a,Plus),ss.SeqPtr(b,Plus) + 1);
  1356.                 }
  1357.                 phase = (phase - exons.back().second + exons.back().first)%3;
  1358.                 if (phase < 0) {
  1359.                     phase += 3;
  1360.                 }
  1361.             } else {
  1362.                 for (int i = 0;  i < exons.size();  ++i) {
  1363.                     int b = exons[i].first;
  1364.                     int a = exons[i].second;
  1365.                     seq.insert(seq.end(),ss.SeqPtr(a,Minus),ss.SeqPtr(b,Minus) + 1);
  1366.                 }
  1367.                 int l = seq.size() - 1;
  1368.                 phase = (phase + exons.back().second - exons.back().first - l)%3;
  1369.                 if (phase < 0) {
  1370.                     phase += 3;
  1371.                 }
  1372.             }
  1374.             for (int i = (3 - phase) % 3;  i + 2 < seq.size();  i += 3) {
  1375.                 if ((seq[i] == nT  &&  seq[i + 1] == nA  &&  seq[i + 2] == nA)  ||  
  1376.                     (seq[i] == nT  &&  seq[i + 1] == nA  &&  seq[i + 2] == nG)  ||  
  1377.                     (seq[i] == nT  &&  seq[i + 1] == nG  &&  seq[i + 2] == nA))  {
  1378.                     NCBI_THROW(CGnomonException, eGenericError,
  1379.                                "GNOMON: Logical error 18");
  1380.                 }
  1381.             }
  1383.             exons.clear();
  1384.             //          cout << "Gene is clearedn";
  1385.         }
  1387.     }
  1388. }
  1391. END_NCBI_SCOPE
  1394. /*
  1395.  * ===========================================================================
  1396.  * $Log: parse.cpp,v $
  1397.  * Revision 1000.2  2004/06/01 18:08:37  gouriano
  1398.  * PRODUCTION: UPGRADED [GCC34_MSVC7] Dev-tree R1.3
  1399.  *
  1400.  * Revision 1.3  2004/05/21 21:41:03  gorelenk
  1401.  * Added PCH ncbi_pch.hpp
  1402.  *
  1403.  * Revision 1.2  2003/11/06 15:02:21  ucko
  1404.  * Use iostream interface from ncbistre.hpp for GCC 2.95 compatibility.
  1405.  *
  1406.  * Revision 1.1  2003/10/24 15:07:25  dicuccio
  1407.  * Initial revision
  1408.  *
  1409.  * ===========================================================================
  1410.  */