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

生物技术

开发平台:

C/C++

  1. /*
  2.  * ===========================================================================
  3.  * PRODUCTION $Log: prot_prop.cpp,v $
  4.  * PRODUCTION Revision 1000.1  2004/06/01 18:10:44  gouriano
  5.  * PRODUCTION PRODUCTION: UPGRADED [GCC34_MSVC7] Dev-tree R1.4
  6.  * PRODUCTION
  7.  * ===========================================================================
  8.  */
  9. /*  $Id: prot_prop.cpp,v 1000.1 2004/06/01 18:10:44 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 <algo/sequence/prot_prop.hpp>
  41. BEGIN_NCBI_SCOPE
  42. BEGIN_SCOPE(objects)
  43. // simply count the occurrences of the different aa (ncbistdaa) types
  44. // puts counts in aacount, returns total count
  45. TSeqPos CProtProp::AACount(CSeqVector& v, vector<TSeqPos>& aacount)
  46. {
  47.     v.SetCoding(CSeq_data::e_Ncbistdaa);
  48.     
  49.     TSeqPos size = v.size();
  50.     aacount.resize(26);
  51.     for (int i = 0;  i < 26;  i++) {
  52.         aacount[i] = 0;
  53.     }
  54.     for (CSeqVector_CI vit(v);  vit.GetPos() < size;  ++vit) {
  55.         CSeqVector::TResidue res = *vit;
  56.         aacount[res]++;
  57.     }
  58.     return size;
  59. }
  60. double CProtProp::GetProteinPI(CSeqVector& v)
  61. {
  62.     // amino acid count
  63.     vector<TSeqPos> aacount;
  64.     CProtProp::AACount(v, aacount);
  65.     // first and last residues
  66.     CSeqVector::TResidue nter = v[0];
  67.     CSeqVector::TResidue cter = v[v.size() - 1];
  68.     // use these to get pI
  69.     
  70. #define PH_MAX 14
  71. #define PH_MIN 0
  72. #define MAXLOOP 2000  // max. number of iterations for root-finding
  73. #define EPSI 0.0001   // epsilon (desired precision)
  74.     
  75.     double phMin = PH_MIN;
  76.     double phMax = PH_MAX;
  77.     double phMid;
  78.     for (int i = 0;  i < MAXLOOP && (phMax - phMin) > EPSI;  i++) {
  79.         phMid = phMin + (phMax - phMin) / 2;
  80.         double charge = GetProteinCharge(aacount, nter, cter, phMid);
  81.         if (charge > 0) {
  82.             phMin = phMid;
  83.         } else {
  84.             phMax = phMid;
  85.         }
  86.     }
  87.     return phMid;
  88. }
  89. //  pKas according to Expasy, by NCBIstdaa:
  90. //  A  B  C  D  E  F  G  H  I  K  L  M  N  P  Q  R  S  T  V  W  X  Y  Z  U  *
  91. //  note: we don't have values for U; use most common for N and C, ignore sidechain
  92. static const double pKaN[26] =
  93. {0, 7.59, 7.5, 7.5, 7.5, 7.7, 7.5, 7.5, 7.5, 7.5, 7.5, 7.5, 7, 7.5, 8.36, 7.5, 7.5, 6.93, 6.82, 7.44, 7.5, 7.5, 7.5, 7.5, 7.5, 0};
  94. static const double pKaC[26] =
  95. {0, 3.55, 3.55, 3.55, 4.55, 4.75, 3.55, 3.55, 3.55, 3.55, 3.55, 3.55, 3.55, 3.55, 3.55, 3.55, 3.55, 3.55, 3.55, 3.55, 3.55, 3.55, 3.55, 3.55, 3.55, 0};
  96. static const double pKaSide[26] =
  97. {0, 0, 0, 9, 4.05, 4.45, 0, 0, 5.98, 0, 10, 0, 0, 0, 0, 0, 12.0, 0, 0, 0, 0, 0, 10, 0, 0, 0};
  98. static const size_t maxRes = sizeof(pKaN) / sizeof(*pKaN) - 1;
  99. double CProtProp::GetProteinCharge(const vector<TSeqPos>& aacount,
  100.                         CSeqVector::TResidue nter,
  101.                         CSeqVector::TResidue cter,
  102.                         double pH)
  103. {
  104.     // N- and C-termini
  105.     double cnter = exp10(-pH) / (exp10(-pKaN[nter]) + exp10(-pH));
  106.     double ccter = exp10(-pKaC[cter]) / (exp10(-pKaC[cter]) + exp10(-pH));
  107.     // basic aas
  108.     double carg = aacount[16] * exp10(-pH) / (exp10(-pKaSide[16]) + exp10(-pH));
  109.     double clys = aacount[10] * exp10(-pH) / (exp10(-pKaSide[10]) + exp10(-pH));
  110.     double chis = aacount[8] * exp10(-pH) / (exp10(-pKaSide[8]) + exp10(-pH));
  111.     // classic acidic aas
  112.     double casp = aacount[4] * exp10(-pKaSide[4]) / (exp10(-pKaSide[4]) + exp10(-pH));
  113.     double cglu = aacount[5] * exp10(-pKaSide[5]) / (exp10(-pKaSide[5]) + exp10(-pH));
  114.     // mildly acidic aas
  115.     double ccys = aacount[3] * exp10(-pKaSide[3]) / (exp10(-pKaSide[3]) + exp10(-pH));
  116.     double ctyr = aacount[22] * exp10(-pKaSide[22]) / (exp10(-pKaSide[22]) + exp10(-pH));
  117.     return carg + clys + chis + cnter -
  118.         (casp + cglu + ctyr + ccys + ccter);
  119. }
  120. END_SCOPE(objects)
  121. END_NCBI_SCOPE
  122. /*
  123.  * ===========================================================================
  124.  * $Log: prot_prop.cpp,v $
  125.  * Revision 1000.1  2004/06/01 18:10:44  gouriano
  126.  * PRODUCTION: UPGRADED [GCC34_MSVC7] Dev-tree R1.4
  127.  *
  128.  * Revision 1.4  2004/05/21 21:41:04  gorelenk
  129.  * Added PCH ncbi_pch.hpp
  130.  *
  131.  * Revision 1.3  2004/05/14 01:24:15  jcherry
  132.  * Pass vector by reference
  133.  *
  134.  * Revision 1.2  2003/07/01 19:01:13  ucko
  135.  * Fix scope use
  136.  *
  137.  * Revision 1.1  2003/07/01 15:10:40  jcherry
  138.  * Initial versions of nuc_prop and prot_prop
  139.  *
  140.  * ===========================================================================
  141.  */