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

生物技术

开发平台:

C/C++

  1. /*
  2.  * ===========================================================================
  3.  * PRODUCTION $Log: msms.cpp,v $
  4.  * PRODUCTION Revision 1000.2  2004/06/01 18:09:02  gouriano
  5.  * PRODUCTION PRODUCTION: UPGRADED [GCC34_MSVC7] Dev-tree R1.6
  6.  * PRODUCTION
  7.  * ===========================================================================
  8.  */
  9. /* $Id: msms.cpp,v 1000.2 2004/06/01 18:09:02 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 authors in any work or product based on this material.
  31.  *
  32.  * ===========================================================================
  33.  *
  34.  * Authors:  Lewis Y. Geer
  35.  *
  36.  * File Description:
  37.  *    Helper classes for ms search algorithms
  38.  *
  39.  * ===========================================================================
  40.  */
  41. #include <ncbi_pch.hpp>
  42. #include <corelib/ncbistd.hpp>
  43. #include <fstream>
  44. #include "msms.hpp"
  45. #include "Mod.hpp"
  46. USING_NCBI_SCOPE;
  47. USING_SCOPE(objects);
  48. USING_SCOPE(omssa);
  49. /////////////////////////////////////////////////////////////////////////////
  50. //
  51. //  CAA::
  52. //
  53. //  lookup table for AA index
  54. //
  55. CAA::CAA(void) 
  56. {
  57.     int i;
  58.     for(i = 0; i < 256; i++) AAMap[i] = kNumUniqueAA;
  59.     for(i = 0; i < kNumUniqueAA; i++) {
  60. AAMap[UniqueAA[i]] = i;
  61. // deal with blast readdb encoding.
  62. AAMap[i] = i;
  63.     }
  64. }
  65. /////////////////////////////////////////////////////////////////////////////
  66. //
  67. //  CCleave::
  68. //
  69. //  Classes for cleaving sequences quickly and computing masses
  70. //
  71. CCleave::CCleave(void): CleaveAt(0), kCleave(0)
  72. {
  73.     ProtonMass = static_cast <int> (kProton*MSSCALE);
  74.     TermMass = static_cast <int> 
  75. ((kTermMass[kCIon] + kTermMass[kXIon])*MSSCALE);
  76.     Reverse = ReverseAA.GetMap();
  77. }
  78. // char based replacement for find_first_of()
  79. int CCleave::findfirst(char* Seq, int Pos, int SeqLen)
  80. {
  81.     int i, j;
  82.     for (i = Pos; i < SeqLen; i++)
  83. for(j = 0; j < kCleave; j++) 
  84.     if(Seq[i] == CleaveAt[j]) return i;
  85.     return i;
  86. }
  87. void CTrypsin::Cleave(char *Seq, int SeqLen, TCleave& Positions)
  88. {
  89.     int Pos = 0;
  90.     Positions.clear();
  91.     Pos = findfirst(Seq, Pos, SeqLen);
  92.     Positions.push_back(0);  // beginning of sequence
  93.     while(Pos < SeqLen - 1) {
  94. if(Seq[Pos+1] == 'P' ||
  95.    Seq[Pos+1] == 'x0e' ) { // not before proline
  96.     Pos = findfirst(Seq, Pos+1, SeqLen);
  97.     continue;
  98. }
  99. Positions.push_back(Pos);
  100. Pos = findfirst(Seq, Pos+1, SeqLen);
  101.     }
  102. }
  103. // - cuts trypsin and calculates mass using integer arithmetic
  104. // - needs to take variable mods into account (i.e. different masses)
  105. // - MassArray contains corrected masses for the fixed mods that
  106. // are not position specific
  107. // - FixedMods contain position specific mods.
  108. //   - open question on how to deal with C terminal fixed mods.  E.g. how do
  109. //     you know you are at the Cterm K?  Doesn't matter for mass, as ANY k will
  110. //     increase total mass, but does matter for ladder.
  111. // - charge is an array of ints for indeterminate charge states
  112. // - the api will also eventually create or extend CLadders for large 
  113. //   search datasets, as most peptides will be examined.
  114. // dealing with missed cleavages:
  115. // - starts with existing mass array
  116. // - extends and shifts existing ladders
  117. //
  118. // - needs to be made into a general method
  119. //
  120. // returns true on end of sequence
  121. // note that the coordinates are inclusive, i.e. [start, end]
  122. bool CTrypsin::CalcAndCut(const char *SeqStart, 
  123.   const char *SeqEnd,  // the end, not beyond the end
  124.   const char **PepStart,  // return value
  125.   int *Masses,
  126.   int& NumMod,
  127.   int MaxNumMod,
  128.   int *EndMasses,
  129.   CMSMod &VariableMods,
  130.   const char **Site,
  131.   int *DeltaMass,
  132.   const int *IntCalcMass  // array of int AA masses
  133.   )
  134. {
  135.     char SeqChar;
  136.     // iterator thru mods
  137.     CMSRequest::TVariable::iterator iMods;
  138.     // iterate thru mods characters
  139.     int iChar;
  140.     // iterate through sequence
  141.     // note that this loop doesn't check at the end of the sequence
  142.     for(; *PepStart < SeqEnd; (*PepStart)++) {
  143. SeqChar = **PepStart;
  144. // check for mods that are type AA only
  145. for(iMods = VariableMods.GetAAMods(eModAA).begin();
  146.     iMods !=  VariableMods.GetAAMods(eModAA).end(); iMods++) {
  147.     for(iChar = 0; iChar < NumModChars[*iMods]; iChar++) {
  148. if (SeqChar == ModChar[iChar][*iMods] && NumMod < MaxNumMod) {
  149.     Site[NumMod] = *PepStart;
  150.     DeltaMass[NumMod] = ModMass[*iMods];
  151.     NumMod++; 
  152. }
  153.     }
  154. }
  155. // if((SeqChar == 'M' || SeqChar == 'x0c') && NumMod < MaxNumMod) {
  156. //     NumMod++;
  157. //     Masses[NumCleave][NumMod-1] = Masses[NumCleave][NumMod-2] +
  158. // 16*MSSCALE;
  159. // }
  160. CalcMass(SeqChar, Masses, IntCalcMass);
  161. // check for cleavage point
  162. if(SeqChar == CleaveAt[0] ||  
  163.    SeqChar == CleaveAt[1] ) { 
  164.     if(*(*PepStart+1) == 'x0e' )  continue;  // not before proline
  165.     EndMass(EndMasses);
  166.     return false;
  167. }
  168.     }
  169.     // todo: deal with mods on the end
  170.     CalcMass(**PepStart, Masses, IntCalcMass);
  171.     EndMass(EndMasses);
  172.     return true;  // end of sequence
  173. }
  174. void CCNBr::Cleave(char *Seq, int SeqLen, TCleave& Positions)
  175. {
  176.     int Pos = 0;
  177.     Positions.clear();
  178.     Positions.push_back(0);  // beginning of sequence
  179.     Pos = findfirst(Seq, Pos, SeqLen);
  180.     while(Pos < SeqLen) {
  181. Positions.push_back(Pos);
  182. Pos = findfirst(Seq, Pos+1, SeqLen);
  183.     }
  184. }
  185. void CFormicAcid::Cleave(char *Seq, int SeqLen, TCleave& Positions)
  186. {
  187.     int Pos = 0;
  188.     Positions.clear();
  189.     Pos = findfirst(Seq, Pos, SeqLen);
  190.     Positions.push_back(0);  // beginning of sequence
  191.     while(Pos < SeqLen - 1) {
  192. if(Seq[Pos+1] != 'P' &&
  193.    Seq[Pos+1] != 'x0e' ) { //  before proline
  194.     Pos = findfirst(Seq, Pos+1, SeqLen);
  195.     continue;
  196. }
  197. Positions.push_back(Pos);
  198. Pos = findfirst(Seq, Pos+1, SeqLen);
  199.     }
  200. }
  201. /////////////////////////////////////////////////////////////////////////////
  202. //
  203. //  CMassArray::
  204. //
  205. //  Holds AA indexed mass array
  206. //
  207. void CMassArray::Init(const CMSRequest::TSearchtype &SearchType)
  208. {
  209.     x_Init(SearchType);
  210. }
  211. void CMassArray::x_Init(const CMSRequest::TSearchtype &SearchType)
  212. {
  213.     int i;
  214.     if(SearchType == eMSSearchType_average) {
  215. for(i = 0; i < kNumUniqueAA; i++ ) {
  216.     CalcMass[i] = AverageMass[i];
  217.     IntCalcMass[i] = static_cast <int> 
  218. (AverageMass[i]*MSSCALE);
  219. }
  220.     }
  221.     else if(SearchType == eMSSearchType_monoisotopic) {
  222. for(i = 0; i < kNumUniqueAA; i++ ) {
  223.     CalcMass[i] = MonoMass[i];
  224.     IntCalcMass[i] = static_cast <int>
  225. (MonoMass[i]*MSSCALE);
  226. }
  227.     }
  228. }
  229. // set up the mass array with fixed mods
  230. void CMassArray::Init(const CMSRequest::TFixed &Mods, 
  231.       const CMSRequest::TSearchtype &SearchType)
  232. {
  233.     x_Init(SearchType);
  234.     CMSRequest::TFixed::const_iterator i;  // iterate thru fixed mods
  235.     int j; // the number of characters affected by the fixed mod
  236.     for(i = Mods.begin(); i != Mods.end(); i++) {
  237. for(j = 0; j < NumModChars[*i]; j++) {
  238.     CalcMass[ModChar[j][*i]] +=  ModMass[*i]/(double)MSSCALE;
  239.     IntCalcMass[ModChar[j][*i]] += ModMass[*i];
  240. }
  241.     }
  242. }
  243. /*
  244.   $Log: msms.cpp,v $
  245.   Revision 1000.2  2004/06/01 18:09:02  gouriano
  246.   PRODUCTION: UPGRADED [GCC34_MSVC7] Dev-tree R1.6
  247.   Revision 1.6  2004/05/21 21:41:03  gorelenk
  248.   Added PCH ncbi_pch.hpp
  249.   Revision 1.5  2004/03/30 19:36:59  lewisg
  250.   multiple mod code
  251.   Revision 1.4  2004/03/16 20:18:54  gorelenk
  252.   Changed includes of private headers.
  253.   Revision 1.3  2004/03/01 18:24:07  lewisg
  254.   better mod handling
  255.   Revision 1.2  2003/10/21 21:12:16  lewisg
  256.   reorder headers
  257.   Revision 1.1  2003/10/20 21:32:13  lewisg
  258.   ommsa toolkit version
  259.   Revision 1.11  2003/10/07 18:02:28  lewisg
  260.   prep for toolkit
  261.   Revision 1.10  2003/08/14 23:49:22  lewisg
  262.   first pass at variable mod
  263.   Revision 1.9  2003/07/17 18:45:49  lewisg
  264.   multi dta support
  265.   Revision 1.8  2003/03/21 21:14:40  lewisg
  266.   merge ming's code, other stuff
  267.   Revision 1.7  2003/02/10 19:37:55  lewisg
  268.   perf and web page cleanup
  269.   Revision 1.6  2003/01/21 21:55:51  lewisg
  270.   fixes
  271.   Revision 1.5  2003/01/21 21:46:13  lewisg
  272.   *** empty log message ***
  273.   Revision 1.4  2002/11/26 00:41:57  lewisg
  274.   changes for msfilter
  275.   Revision 1.3  2002/09/20 20:19:34  lewisg
  276.   msms search update
  277.   Revision 1.2  2002/07/16 13:26:23  lewisg
  278.   *** empty log message ***
  279.   Revision 1.1.1.1  2002/02/14 02:14:02  lewisg
  280. */