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

生物技术

开发平台:

C/C++

  1. /*
  2.  * ===========================================================================
  3.  * PRODUCTION $Log: nuc_prop.cpp,v $
  4.  * PRODUCTION Revision 1000.1  2004/06/01 18:10:37  gouriano
  5.  * PRODUCTION PRODUCTION: UPGRADED [GCC34_MSVC7] Dev-tree R1.6
  6.  * PRODUCTION
  7.  * ===========================================================================
  8.  */
  9. /*  $Id: nuc_prop.cpp,v 1000.1 2004/06/01 18:10:37 gouriano Exp $
  10.  * ===========================================================================
  11.  *
  12.  *                            PUBLIC DOMAIN NOTICE
  13.  *               National Center for Biotechnology Information
  14.  *
  15.  *  This software/database is a "United States Government Work" under the
  16.  *  terms of the United States Copyright Act.  It was written as part of
  17.  *  the author's official duties as a United States Government employee and
  18.  *  thus cannot be copyrighted.  This software/database is freely available
  19.  *  to the public for use. The National Library of Medicine and the U.S.
  20.  *  Government have not placed any restriction on its use or reproduction.
  21.  *
  22.  *  Although all reasonable efforts have been taken to ensure the accuracy
  23.  *  and reliability of the software and data, the NLM and the U.S.
  24.  *  Government do not and cannot warrant the performance or results that
  25.  *  may be obtained by using this software or data. The NLM and the U.S.
  26.  *  Government disclaim all warranties, express or implied, including
  27.  *  warranties of performance, merchantability or fitness for any particular
  28.  *  purpose.
  29.  *
  30.  *  Please cite the author in any work or product based on this material.
  31.  *
  32.  * ===========================================================================
  33.  *
  34.  * Authors:  Josh Cherry
  35.  *
  36.  * File Description:
  37.  *
  38.  */
  39. #include <ncbi_pch.hpp>
  40. #include <objects/seq/seqport_util.hpp>
  41. #include <algo/sequence/nuc_prop.hpp>
  42. BEGIN_NCBI_SCOPE
  43. BEGIN_SCOPE(objects)
  44. void CNucProp::CountNmers(CSeqVector& seqvec, int n, vector<int>& table)
  45. {
  46.     TSeqPos len = seqvec.size();
  47.     table.resize(NumberOfNmers(n));
  48.     // clear table
  49.     for (int i = 0;  i < NumberOfNmers(n);  i++) {
  50.         table[i] = 0;
  51.     }
  52.     string seq_string;
  53.     seqvec.GetSeqData(0, len, seq_string);
  54.     const char *seq;
  55.     seq = seq_string.data();
  56.     for (TSeqPos i = 0;  i <= len-n;  ++i) {
  57.         int nmerint = Nmer2Int(seq+i, n);
  58.         if (nmerint >= 0) {   // if no ambiguity chars
  59.             table[nmerint]++;
  60.         }
  61.     }
  62. }
  63. // convert nmer of length n, pointed to by seq,
  64. // to an integer representation
  65. // if there's a character other than AGCT, return -1
  66. // n better be small enough that the result fits into an int
  67. int CNucProp::Nmer2Int(const char *seq, int n)
  68. {
  69.     int rv = 0;
  70.     for (int i = 0;  i<n;  i++) {
  71.         rv <<= 2;
  72.         int nucint = Nuc2Nybble(seq[i]);
  73.         if(nucint < 0) {
  74.             return -1;
  75.         }
  76.         rv |= nucint;
  77.     }
  78.     return rv;
  79. }
  80. // convert int from Nmer2Int back to a string
  81. void CNucProp::Int2Nmer(int nmer_int, int nmer_size, string& out)
  82. {
  83.     out.resize(nmer_size);
  84.     for (int i = nmer_size-1;  i >= 0; i--) {
  85.         out[i] = Nybble2Nuc(nmer_int & 3);   // analyze two low-order bits
  86.         nmer_int >>= 2;
  87.     }
  88. }
  89. // there are 4^n n-mers
  90. int CNucProp::NumberOfNmers(int n)
  91. {
  92.     return 1 << (2*n);
  93. }
  94. int CNucProp::Nuc2Nybble(char nuc)
  95. {
  96.     switch (nuc) {
  97.     case 'G':
  98.         return 0;
  99.     case 'A':
  100.         return 1;
  101.     case 'T':
  102.         return 2;
  103.     case 'C':
  104.         return 3;
  105.     default:     // other than GATC
  106.         return -1;
  107.     }
  108. }
  109. char CNucProp::Nybble2Nuc(int n)
  110. {
  111.     switch (n) {
  112.     case 0:
  113.         return 'G';
  114.     case 1:
  115.         return 'A';
  116.     case 2:
  117.         return 'T';
  118.     case 3:
  119.         return 'C';
  120.     default:
  121.         // this should never happen
  122.         return '?';
  123.     }
  124. }
  125. int CNucProp::GetPercentGC(const CSeqVector& seqvec)
  126. {
  127.     TSeqPos gc_count = 0;
  128.     TSeqPos len = seqvec.size();
  129.     for (TSeqPos i = 0;  i < len;  ++i) {
  130.         switch (seqvec[i]) {
  131.         case 'C':
  132.         case 'G':
  133.         case 'S':
  134.             ++gc_count;
  135.             break;
  136.         default:
  137.             break;
  138.         }
  139.     }
  140.     return (int) ((gc_count * 100.0) / len + 0.5);
  141. }
  142. END_SCOPE(objects)
  143. END_NCBI_SCOPE
  144. /*
  145.  * ===========================================================================
  146.  * $Log: nuc_prop.cpp,v $
  147.  * Revision 1000.1  2004/06/01 18:10:37  gouriano
  148.  * PRODUCTION: UPGRADED [GCC34_MSVC7] Dev-tree R1.6
  149.  *
  150.  * Revision 1.6  2004/05/21 21:41:04  gorelenk
  151.  * Added PCH ncbi_pch.hpp
  152.  *
  153.  * Revision 1.5  2003/08/18 17:35:29  jcherry
  154.  * Fixed CountNmers to alter result vectors size, not just its capacity.
  155.  * Made Int2Nmer build its string from scratch, and do so properly.
  156.  *
  157.  * Revision 1.4  2003/07/28 20:41:01  jcherry
  158.  * Changed GetPercentGC() to round properly
  159.  *
  160.  * Revision 1.3  2003/07/28 11:54:34  dicuccio
  161.  * Changed Int2Nmer to use std::string instead of char*
  162.  *
  163.  * Revision 1.2  2003/07/01 19:01:13  ucko
  164.  * Fix scope use
  165.  *
  166.  * Revision 1.1  2003/07/01 15:10:40  jcherry
  167.  * Initial versions of nuc_prop and prot_prop
  168.  *
  169.  * ===========================================================================
  170.  */