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

生物技术

开发平台:

C/C++

  1. /*
  2.  * ===========================================================================
  3.  * PRODUCTION $Log: polya.hpp,v $
  4.  * PRODUCTION Revision 1000.2  2004/06/01 18:04:31  gouriano
  5.  * PRODUCTION PRODUCTION: UPGRADED [GCC34_MSVC7] Dev-tree R1.6
  6.  * PRODUCTION
  7.  * ===========================================================================
  8.  */
  9. /*  $Id: polya.hpp,v 1000.2 2004/06/01 18:04:31 gouriano Exp $
  10. * ===========================================================================
  11. *
  12. *                            PUBLIC DOMAIN NOTICE
  13. *               National Center for Biotechnology Information
  14. *
  15. *  This software/database is a "United States Government Work" under the
  16. *  terms of the United States Copyright Act.  It was written as part of
  17. *  the author's official duties as a United States Government employee and
  18. *  thus cannot be copyrighted.  This software/database is freely available
  19. *  to the public for use. The National Library of Medicine and the U.S.
  20. *  Government have not placed any restriction on its use or reproduction.
  21. *
  22. *  Although all reasonable efforts have been taken to ensure the accuracy
  23. *  and reliability of the software and data, the NLM and the U.S.
  24. *  Government do not and cannot warrant the performance or results that
  25. *  may be obtained by using this software or data. The NLM and the U.S.
  26. *  Government disclaim all warranties, express or implied, including
  27. *  warranties of performance, merchantability or fitness for any particular
  28. *  purpose.
  29. *
  30. *  Please cite the author in any work or product based on this material.
  31. *
  32. * ===========================================================================
  33. *
  34. * Author: Philip Johnson
  35. *
  36. * File Description: finds mRNA 3' modification (poly-A tails)
  37. *
  38. * ---------------------------------------------------------------------------
  39. */
  40. #ifndef ALGO_SEQUENCE___POLYA__HPP
  41. #define ALGO_SEQUENCE___POLYA__HPP
  42. #include <corelib/ncbistd.hpp>
  43. BEGIN_NCBI_SCOPE
  44. enum EPolyTail {
  45.     ePolyTail_None,
  46.     ePolyTail_A3, //> 3' poly-A tail
  47.     ePolyTail_T5  //> 5' poly-T head (submitted to db reversed?)
  48. };
  49. ///////////////////////////////////////////////////////////////////////////////
  50. /// PRE : two random access iterators pointing to sequence data [begin,
  51. /// end)
  52. /// POST: poly-A tail cleavage site, if any (-1 if not)
  53. template <typename Iterator>
  54. TSignedSeqPos
  55. FindPolyA(Iterator begin, Iterator end);
  56. ///////////////////////////////////////////////////////////////////////////////
  57. /// PRE : two random access iterators pointing to sequence data [begin,
  58. /// end)
  59. /// POST: cleavageSite (if any) and whether we found a poly-A tail, a poly-T
  60. /// head, or neither
  61. template <typename Iterator>
  62. EPolyTail
  63. FindPolyTail(Iterator begin, Iterator end, TSignedSeqPos &cleavageSite);
  64. ///////////////////////////////////////////////////////////////////////////////
  65. /// Implementation [in header because of templates]
  66. template<typename Iterator>
  67. class CRevComp_It {
  68. public:
  69.     CRevComp_It(void) {}
  70.     CRevComp_It(const Iterator &it) {
  71.         m_Base = it;
  72.     }
  73.     char operator*(void) const {
  74.         Iterator tmp = m_Base;
  75.         switch (*--tmp) {
  76.         case 'A': return 'T';
  77.         case 'T': return 'A';
  78.         case 'C': return 'G';
  79.         case 'G': return 'C';
  80.         default: return *tmp;
  81.         }
  82.     }
  83.     CRevComp_It& operator++(void) {
  84.         --m_Base;
  85.         return *this;
  86.     }
  87.     CRevComp_It operator++(int) {
  88.         CRevComp_It it = m_Base;
  89.         --m_Base;
  90.         return it;
  91.     }
  92.     CRevComp_It& operator+=(int i) {
  93.         m_Base -= i;
  94.         return *this;
  95.     }
  96.     CRevComp_It& operator-=(int i) {
  97.         m_Base += i;
  98.         return *this;
  99.     }
  100.     CRevComp_It operator+ (int i) const {
  101.         CRevComp_It it(m_Base);
  102.         it += i;
  103.         return it;
  104.     }
  105.     CRevComp_It operator- (int i) const {
  106.         CRevComp_It it(m_Base);
  107.         it -= i;
  108.         return it;
  109.     }
  110.     int operator- (const CRevComp_It &it) const {
  111.         return it.m_Base - m_Base;
  112.     }
  113.     
  114.     //booleans
  115.     bool operator>= (const CRevComp_It &it) const {
  116.         return m_Base <= it.m_Base;
  117.     }
  118.     bool operator>  (const CRevComp_It &it) const {
  119.         return m_Base < it.m_Base;
  120.     }
  121.     bool operator<= (const CRevComp_It &it) const {
  122.         return m_Base >= it.m_Base;
  123.     }
  124.     bool operator<  (const CRevComp_It &it) const {
  125.         return m_Base > it.m_Base;
  126.     }
  127.     bool operator== (const CRevComp_It &it) const {
  128.         return m_Base == it.m_Base;
  129.     }
  130.     bool operator!= (const CRevComp_It &it) const {
  131.         return m_Base != it.m_Base;
  132.     }
  133. private:
  134.     Iterator m_Base;
  135. };
  136. ///////////////////////////////////////////////////////////////////////////////
  137. // PRE : same conditions as STL 'search', but iterators must have ptrdiff_t
  138. // difference type
  139. // POST: same as STL 'search'
  140. template <typename ForwardIterator1, typename ForwardIterator2>
  141. ForwardIterator1 ItrSearch(ForwardIterator1 first1, ForwardIterator1 last1,
  142.                            ForwardIterator2 first2, ForwardIterator2 last2)
  143. {
  144.     ptrdiff_t d1 = last1 - first1;
  145.     ptrdiff_t d2 = last2 - first2;
  146.     if (d1 < d2) {
  147.         return last1;
  148.     }
  149.     ForwardIterator1 current1 = first1;
  150.     ForwardIterator2 current2 = first2;
  151.     while (current2 != last2) {
  152.         if (!(*current1 == *current2)) {
  153.             if (d1-- == d2) {
  154.                 return last1;
  155.             } else {
  156.                 current1 = ++first1;
  157.                 current2 = first2;
  158.             }
  159.         } else {
  160.             ++current1;
  161.             ++current2;
  162.         }
  163.     }
  164.     return (current2 == last2) ? first1 : last1;
  165. }
  166. ///////////////////////////////////////////////////////////////////////////////
  167. // PRE : two random access iterators pointing to sequence data [begin,
  168. // end)
  169. // POST: poly-A tail cleavage site, if any (-1 if not)
  170. template <typename Iterator>
  171. TSignedSeqPos FindPolyA(Iterator begin, Iterator end)
  172. {
  173.     string motif1("AATAAA");
  174.     string motif2("ATTAAA");
  175.     Iterator pos = max(begin, end-250);
  176.     Iterator uStrmMotif = pos;
  177.     while (uStrmMotif != end) {
  178.         pos = uStrmMotif;
  179.         uStrmMotif = ItrSearch(pos, end, motif1.begin(), motif1.end());
  180.         if (uStrmMotif == end) {
  181.             uStrmMotif = ItrSearch(pos, end, motif2.begin(), motif2.end());
  182.         }
  183.         if (uStrmMotif != end) {
  184.             uStrmMotif += 6; // skip over upstream motif
  185.             pos = uStrmMotif;
  186.             if (end - pos < 10) { //min skip of 10
  187.                 break;
  188.             }
  189.             Iterator maxCleavage = (end - pos < 30) ? end : pos + 30;
  190.             unsigned int aRun = 0;
  191.             for (pos += 10;  pos < maxCleavage  &&  aRun < 3;  ++pos) {
  192.                 if (*pos == 'A') {
  193.                     ++aRun;
  194.                 } else {
  195.                     aRun = 0;
  196.                 }
  197.             }
  198.             
  199.             if (aRun) {
  200.                 pos -= aRun;
  201.             }
  202.             TSignedSeqPos cleavageSite = pos - begin;
  203.             //now let's look for poly-adenylated tail..
  204.             unsigned int numA = 0, numOther = 0;
  205.             while (pos < end) {
  206.                 if (*pos == 'A') {
  207.                     ++numA;
  208.                 } else {
  209.                     ++numOther;
  210.                 }
  211.                 ++pos;
  212.             }
  213.             if (numOther + numA > 0  &&
  214.                 ((double) numA / (numA+numOther)) > 0.95) {
  215.                 return cleavageSite;
  216.             }
  217.         }
  218.     }
  219.     return -1;
  220. }
  221. ///////////////////////////////////////////////////////////////////////////////
  222. // PRE : two random access iterators pointing to sequence data [begin,
  223. // end)
  224. // POST: cleavageSite (if any) and whether we found a poly-A tail, a poly-T
  225. // head, or neither
  226. template<typename Iterator>
  227. EPolyTail FindPolyTail(Iterator begin, Iterator end,
  228.                        TSignedSeqPos &cleavageSite)
  229. {
  230.     cleavageSite = FindPolyA(begin, end);
  231.     if (cleavageSite >= 0) {
  232.         return ePolyTail_A3;
  233.     } else {
  234.         int seqLen = end - begin;
  235.         cleavageSite = FindPolyA(CRevComp_It<Iterator>(end),
  236.                                  CRevComp_It<Iterator>(begin));
  237.         if (cleavageSite >= 0) {
  238.             cleavageSite = seqLen - cleavageSite - 1;
  239.             return ePolyTail_T5;
  240.         }
  241.     }
  242.     return ePolyTail_None;
  243. }
  244. END_NCBI_SCOPE
  245. /*
  246. * ===========================================================================
  247. * $Log: polya.hpp,v $
  248. * Revision 1000.2  2004/06/01 18:04:31  gouriano
  249. * PRODUCTION: UPGRADED [GCC34_MSVC7] Dev-tree R1.6
  250. *
  251. * Revision 1.6  2004/05/12 19:13:36  johnson
  252. * make sure no iterator ever goes past end
  253. *
  254. * Revision 1.5  2004/05/11 18:36:56  johnson
  255. * made reverse iterator more standard (== base-1); avoid add a value to the
  256. * iterator that could potentially take it past the end
  257. *
  258. * Revision 1.4  2004/04/28 15:19:46  johnson
  259. * Uses templated iterators; no longer accepts user suggestion for cleavage
  260. * site
  261. *
  262. * Revision 1.3  2003/12/31 20:41:39  johnson
  263. * FindPolySite takes cleavage prompt for both 3' poly-A and 5' poly-T
  264. *
  265. * Revision 1.2  2003/12/30 21:28:31  johnson
  266. * added msvc export specifiers
  267. *
  268. * ===========================================================================
  269. */
  270. #endif /*ALGO_SEQUENCE___POLYA__HPP*/