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

生物技术

开发平台:

C/C++

  1. /*
  2.  * ===========================================================================
  3.  * PRODUCTION $Log: models.hpp,v $
  4.  * PRODUCTION Revision 1000.0  2003/10/29 19:38:58  gouriano
  5.  * PRODUCTION PRODUCTION: IMPORTED [ORIGINAL] Dev-tree R1.1
  6.  * PRODUCTION
  7.  * ===========================================================================
  8.  */
  9. #ifndef __MODELS__HPP
  10. #define __MODELS__HPP
  11. /*  $Id: models.hpp,v 1000.0 2003/10/29 19:38:58 gouriano Exp $
  12.  * ===========================================================================
  13.  *
  14.  *                            PUBLIC DOMAIN NOTICE
  15.  *               National Center for Biotechnology Information
  16.  *
  17.  *  This software/database is a "United States Government Work" under the
  18.  *  terms of the United States Copyright Act.  It was written as part of
  19.  *  the author's official duties as a United States Government employee and
  20.  *  thus cannot be copyrighted.  This software/database is freely available
  21.  *  to the public for use. The National Library of Medicine and the U.S.
  22.  *  Government have not placed any restriction on its use or reproduction.
  23.  *
  24.  *  Although all reasonable efforts have been taken to ensure the accuracy
  25.  *  and reliability of the software and data, the NLM and the U.S.
  26.  *  Government do not and cannot warrant the performance or results that
  27.  *  may be obtained by using this software or data. The NLM and the U.S.
  28.  *  Government disclaim all warranties, express or implied, including
  29.  *  warranties of performance, merchantability or fitness for any particular
  30.  *  purpose.
  31.  *
  32.  *  Please cite the author in any work or product based on this material.
  33.  *
  34.  * ===========================================================================
  35.  *
  36.  * Authors:  Alexandre Souvorov
  37.  *
  38.  * File Description:
  39.  *
  40.  */
  41. #include <corelib/ncbistd.hpp>
  42. BEGIN_NCBI_SCOPE
  43. template<int order> void MarkovChain<order>::InitScore(CNcbiIfstream& from)
  44. {
  45.     Init(from);
  46.     if(from) toScore();
  47. }
  48. template<int order> void MarkovChain<order>::Init(CNcbiIfstream& from)
  49. {
  50.     next[nA].Init(from);
  51.     next[nC].Init(from);
  52.     next[nG].Init(from);
  53.     next[nT].Init(from);
  54.     if(from) next[nN].Average(next[nA],next[nC],next[nG],next[nT]);
  55. }
  56. template<int order> void MarkovChain<order>::Average(Type& mc0, Type& mc1, Type& mc2, Type& mc3)
  57. {
  58.     next[nA].Average(mc0.next[nA],mc1.next[nA],mc2.next[nA],mc3.next[nA]);
  59.     next[nC].Average(mc0.next[nC],mc1.next[nC],mc2.next[nC],mc3.next[nC]);
  60.     next[nG].Average(mc0.next[nG],mc1.next[nG],mc2.next[nG],mc3.next[nG]);
  61.     next[nT].Average(mc0.next[nT],mc1.next[nT],mc2.next[nT],mc3.next[nT]);
  62.     next[nN].Average(next[nA],next[nC],next[nG],next[nT]);
  63. }
  64. template<int order> void MarkovChain<order>::toScore()
  65. {
  66.     next[nA].toScore();
  67.     next[nC].toScore();
  68.     next[nG].toScore();
  69.     next[nT].toScore();
  70.     next[nN].toScore();
  71. }
  72. template<int order> inline double MarkovChainArray<order>::Score(const int* seq) const
  73. {
  74.     double score = 0;
  75.     for(int i = 0; i < length; ++i)
  76.     {
  77.         double s = mc[i].Score(seq+i);
  78.         if(s == BadScore) return BadScore;
  79.         score += s;
  80.     }
  81.     return score;
  82. }
  83. template<int order> void MarkovChainArray<order>::InitScore(int l, CNcbiIfstream& from)
  84. {
  85.     length = l;
  86.     mc.resize(length);
  87.     for(int i = 0; i < length; ++i) mc[i].InitScore(from);
  88. }
  89. template<int order>
  90. double WAM_Donor<order>::Score(const IVec& seq, int i) const
  91. {
  92.     int first = i-left+1;
  93.     int last = i+right;
  94.     if(first-order < 0 || last >= seq.size()) return BadScore;    // out of range
  95.     if(seq[i+1] != nG || seq[i+2] != nT) return BadScore;   // no GT
  96.     return matrix.Score(&seq[first]);
  97. }
  98. template<int order>
  99. WAM_Donor<order>::WAM_Donor(const string& file, int cgcontent)
  100. {
  101.     CNcbiOstrstream ost;
  102.     ost << "[WAM_Donor_" << order << "]";
  103.     string label = CNcbiOstrstreamToString(ost);
  104.     CNcbiIfstream from(file.c_str());
  105.     pair<int,int> cgrange = FindContent(from,label,cgcontent);
  106.     if(cgrange.first < 0) Error(label);
  107.     string str;
  108.     from >> str;
  109.     if(str != "InExon:") Error(label);
  110.     from >> inexon;
  111.     if(!from) Error(label);
  112.     from >> str;
  113.     if(str != "InIntron:") Error(label);
  114.     from >> inintron;
  115.     if(!from) Error(label);
  116.     left = inexon;
  117.     right = inintron;
  118.     matrix.InitScore(inexon+inintron,from);
  119.     if(!from) Error(label);
  120. }
  121. template<int order>
  122. double WAM_Acceptor<order>::Score(const IVec& seq, int i) const
  123. {
  124.     int first = i-left+1;
  125.     int last = i+right;
  126.     if(first-order < 0 || last >= seq.size()) return BadScore;  // out of range
  127.     if(seq[i-1] != nA || seq[i] != nG) return BadScore;   // no AG
  128.     return matrix.Score(&seq[first]);
  129. }
  130. template<int order>
  131. WAM_Acceptor<order>::WAM_Acceptor(const string& file, int cgcontent)
  132. {
  133.     CNcbiOstrstream ost;
  134.     ost << "[WAM_Acceptor_" << order << "]";
  135.     string label = CNcbiOstrstreamToString(ost);
  136.     CNcbiIfstream from(file.c_str());
  137.     pair<int,int> cgrange = FindContent(from,label,cgcontent);
  138.     if(cgrange.first < 0) Error(label);
  139.     string str;
  140.     from >> str;
  141.     if(str != "InExon:") Error(label);
  142.     from >> inexon;
  143.     if(!from) Error(label);
  144.     from >> str;
  145.     if(str != "InIntron:") Error(label);
  146.     from >> inintron;
  147.     if(!from) Error(label);
  148.     left = inintron;
  149.     right = inexon;
  150.     matrix.InitScore(inexon+inintron,from);
  151.     if(!from) Error(label);
  152. }
  153. template<int order>
  154. MC3_CodingRegion<order>::MC3_CodingRegion(const string& file, int cgcontent)
  155. {
  156.     CNcbiOstrstream ost;
  157.     ost << "[MC3_CodingRegion_" << order << "]";
  158.     string label = CNcbiOstrstreamToString(ost);
  159.     CNcbiIfstream from(file.c_str());
  160.     pair<int,int> cgrange = FindContent(from,label,cgcontent);
  161.     if(cgrange.first < 0) Error(label);
  162.     matrix[0].InitScore(from);
  163.     matrix[1].InitScore(from);
  164.     matrix[2].InitScore(from);
  165.     if(!from) Error(label);
  166. }
  167. template<int order>
  168. double MC3_CodingRegion<order>::Score(const IVec& seq, int i, int codonshift) const
  169. {
  170.     if(i-order < 0) return BadScore;  // out of range
  171.     else if(seq[i] == nN) return -2.;  // discourage making exons of Ns
  172.     else return matrix[codonshift].Score(&seq[i]);
  173. }
  174. template<int order>
  175. MC_NonCodingRegion<order>::MC_NonCodingRegion(const string& file, int cgcontent)
  176. {
  177.     CNcbiOstrstream ost;
  178.     ost << "[MC_NonCodingRegion_" << order << "]";
  179.     string label = CNcbiOstrstreamToString(ost);
  180.     CNcbiIfstream from(file.c_str());
  181.     pair<int,int> cgrange = FindContent(from,label,cgcontent);
  182.     if(cgrange.first < 0) Error(label);
  183.     matrix.InitScore(from);
  184.     if(!from) Error(label);
  185. }
  186. template<int order>
  187. double MC_NonCodingRegion<order>::Score(const IVec& seq, int i) const
  188. {
  189.     if(i-order < 0) return BadScore;  // out of range
  190.     return matrix.Score(&seq[i]);
  191. }
  192. inline bool SeqScores::isStart(int i, int strand) const
  193. {
  194.     const IVec& ss = seq[strand];
  195.     int ii = (strand == Plus) ? i : SeqLen()-1-i;
  196.     if(ii < 0 || ii+2 >= SeqLen()) return false;  //out of range
  197.     else if(ss[ii] != nA || ss[ii+1] != nT || ss[ii+2] != nG) return false;
  198.     else return true;
  199. }
  200. inline bool SeqScores::isStop(int i, int strand) const
  201. {
  202.     const IVec& ss = seq[strand];
  203.     int ii = (strand == Plus) ? i : SeqLen()-1-i;
  204.     if(ii < 0 || ii+2 >= SeqLen()) return false;  //out of range
  205.     if((ss[ii] != nT || ss[ii+1] != nA || ss[ii+2] != nA) &&
  206.        (ss[ii] != nT || ss[ii+1] != nA || ss[ii+2] != nG) &&
  207.        (ss[ii] != nT || ss[ii+1] != nG || ss[ii+2] != nA)) return false;
  208.     else return true;
  209. }
  210. inline bool SeqScores::isAG(int i, int strand) const
  211. {
  212.     const IVec& ss = seq[strand];
  213.     int ii = (strand == Plus) ? i : SeqLen()-1-i;
  214.     if(ii-1 < 0 || ii >= SeqLen()) return false;  //out of range
  215.     if(ss[ii-1] != nA || ss[ii] != nG) return false;
  216.     else return true;
  217. }
  218. inline bool SeqScores::isGT(int i, int strand) const
  219. {
  220.     const IVec& ss = seq[strand];
  221.     int ii = (strand == Plus) ? i : SeqLen()-1-i;
  222.     if(ii < 0 || ii+1 >= SeqLen()) return false;  //out of range
  223.     if(ss[ii] != nG || ss[ii+1] != nT) return false;
  224.     else return true;
  225. }
  226. inline bool SeqScores::isConsensusIntron(int i, int j, int strand) const
  227. {
  228.     if(strand == Plus) return (dscr[Plus][i-1] != BadScore) && (ascr[Plus][j] != BadScore);
  229.     else return (ascr[Minus][i-1] != BadScore) && (dscr[Minus][j] != BadScore);
  230.     // if(strand == Plus) return isGT(i,strand) && isAG(j,strand);
  231.     // else return isAG(i,strand) && isGT(j,strand);
  232. }
  233. inline const int* SeqScores::SeqPtr(int i, int strand) const
  234. {
  235.     const IVec& ss = seq[strand];
  236.     int ii = (strand == Plus) ? i : SeqLen()-1-i;
  237.     return &ss.front()+ii;
  238. }
  239. inline bool SeqScores::StopInside(int a, int b, int strand, int frame) const
  240. {
  241.     return (a <= laststop[strand][frame][b]);
  242. }
  243. inline bool SeqScores::OpenCodingRegion(int a, int b, int strand, int frame) const
  244. {
  245.     return (a > notinexon[strand][frame][b]);
  246. }
  247. inline double SeqScores::CodingScore(int a, int b, int strand, int frame) const
  248. {
  249.     if(a > b) return 0; // for splitted start/stop
  250.     double score = cdrscr[strand][frame][b];
  251.     if(a > 0) score -= cdrscr[strand][frame][a-1];
  252.     return score;
  253. }
  254. inline bool SeqScores::OpenNonCodingRegion(int a, int b, int strand) const
  255. {
  256.     return (a > notinintron[strand][b]);
  257. }
  258. inline double SeqScores::NonCodingScore(int a, int b, int strand) const
  259. {
  260.     double score = ncdrscr[strand][b];
  261.     if(a > 0) score -= ncdrscr[strand][a-1];
  262.     return score;
  263. }
  264. inline bool SeqScores::OpenIntergenicRegion(int a, int b) const
  265. {
  266.     return (a > notining[b]);
  267. inline bool SeqScores::InAlignment(int a, int b) const
  268. {
  269.     return (inalign[b] >=0 && a >= inalign[b]);
  270. inline double SeqScores::IntergenicScore(int a, int b, int strand) const
  271. {
  272.     double score = ingscr[strand][b];
  273.     if(a > 0) score -= ingscr[strand][a-1];
  274.     return score;
  275. }
  276. inline int SeqScores::SeqMap(int i, bool forwrd) const 
  277.     int l = seq_map[i];
  278.     if(l < 0)
  279.     {
  280.         if(forwrd) return seq_map[i+1];
  281.         else return seq_map[i-1];
  282.     }
  283.     else return l;
  284. }
  285. inline int SeqScores::RevSeqMap(int i, bool forwrd) const 
  286.     int l = rev_seq_map[i-shift];
  287.     if(l < 0)
  288.     {
  289.         if(forwrd) return rev_seq_map[i-shift+1];
  290.         else return rev_seq_map[i-shift-1];
  291.     }
  292.     else return l;
  293. }
  294. inline int ACGT(int c)
  295. {
  296.     switch(c)
  297.     {
  298.     case 'A': 
  299.     case 'a': 
  300.         return nA;
  301.     case 'C': 
  302.     case 'c': 
  303.         return nC;
  304.     case 'G': 
  305.     case 'g': 
  306.         return nG;
  307.     case 'T': 
  308.     case 't': 
  309.         return nT;
  310.     default:
  311.         return nN;
  312.     }
  313. }
  314. END_NCBI_SCOPE
  315. /*
  316.  * ===========================================================================
  317.  * $Log: models.hpp,v $
  318.  * Revision 1000.0  2003/10/29 19:38:58  gouriano
  319.  * PRODUCTION: IMPORTED [ORIGINAL] Dev-tree R1.1
  320.  *
  321.  * Revision 1.1  2003/10/24 15:07:25  dicuccio
  322.  * Initial revision
  323.  *
  324.  * ===========================================================================
  325.  */
  326. #endif  // __MODELS__HPP