开通博客赚积分
发布资源赚积分
分类

源码开发语言/平台
当前位置：首页 > 源码/资料 > 数值算法/人工智能 > 生物技术 > 查看源码
mspeak.cpp
资源名称：ncbi_cxx_win--Jul_08_2004.zip　[点击查看]
上传用户：yhdzpy8989
上传日期：2007-06-13
资源大小：13604k
文件大小：26k
源码类别：

生物技术

开发平台：

C/C++

mspeak.cpp：源码内容
  1. /*
  2.  * ===========================================================================
  3.  * PRODUCTION $Log: mspeak.cpp,v $
  4.  * PRODUCTION Revision 1000.4  2004/06/01 18:09:04  gouriano
  5.  * PRODUCTION PRODUCTION: UPGRADED [GCC34_MSVC7] Dev-tree R1.13
  6.  * PRODUCTION
  7.  * ===========================================================================
  8.  */
  10. /* $Id: mspeak.cpp,v 1000.4 2004/06/01 18:09:04 gouriano Exp $
  11.  * ===========================================================================
  12.  *
  13.  *                            PUBLIC DOMAIN NOTICE
  14.  *               National Center for Biotechnology Information
  15.  *
  16.  *  This software/database is a "United States Government Work" under the
  17.  *  terms of the United States Copyright Act.  It was written as part of
  18.  *  the author's official duties as a United States Government employee and
  19.  *  thus cannot be copyrighted.  This software/database is freely available
  20.  *  to the public for use. The National Library of Medicine and the U.S.
  21.  *  Government have not placed any restriction on its use or reproduction.
  22.  *
  23.  *  Although all reasonable efforts have been taken to ensure the accuracy
  24.  *  and reliability of the software and data, the NLM and the U.S.
  25.  *  Government do not and cannot warrant the performance or results that
  26.  *  may be obtained by using this software or data. The NLM and the U.S.
  27.  *  Government disclaim all warranties, express or implied, including
  28.  *  warranties of performance, merchantability or fitness for any particular
  29.  *  purpose.
  30.  *
  31.  *  Please cite the authors in any work or product based on this material.
  32.  *
  33.  * ===========================================================================
  34.  *
  35.  * Author:  Lewis Y. Geer
  36.  *  
  37.  *
  38.  * File Description:
  39.  *    code to deal with spectra and m/z ladders
  40.  *
  41.  * ===========================================================================
  42.  */
  44. #include <ncbi_pch.hpp>
  45. #include <corelib/ncbistd.hpp>
  46. #include <corelib/ncbi_limits.h>
  48. #include <algorithm>
  50. #include "mspeak.hpp"
  53. USING_NCBI_SCOPE;
  54. USING_SCOPE(objects);
  55. USING_SCOPE(omssa);
  57. /////////////////////////////////////////////////////////////////////////////
  58. //
  59. //  CMSHit::
  60. //
  61. //  Used by CMSPeak class to hold hits
  62. //
  65. // helper function for RecordHits that scans thru a single ladder
  66. void CMSHit::RecordMatchesScan(CLadder& Ladder, int& iHitInfo, CMSPeak *Peaks,
  67.        int Which)
  68. {
  69.     try {
  70. TIntensity Intensity(new unsigned [Ladder.size()]);
  71. Peaks->CompareSorted(Ladder, Which, &Intensity);
  73. // examine hits array
  74. unsigned i;
  75. for(i = 0; i < Ladder.size(); i++) {
  76.     // if hit, add to hitlist
  77.     if(Ladder.GetHit()[i] > 0) {
  78. SetHitInfo(iHitInfo).SetCharge() = (char) Ladder.GetCharge();
  79. SetHitInfo(iHitInfo).SetIon() = (char) Ladder.GetType();
  80. SetHitInfo(iHitInfo).SetNumber() = (short) i;
  81. SetHitInfo(iHitInfo).SetIntensity() = *(Intensity.get() + i);
  82. iHitInfo++;
  83.     }
  84. }
  85.     } catch (NCBI_NS_STD::exception& e) {
  86. ERR_POST(Info << "Exception caught in CMSHit::RecordMatchesScan: " << e.what());
  87. throw;
  88.     }
  90. }
  92. // make a record of the ions hit
  93. void CMSHit::RecordMatches(CLadder& BLadder, CLadder& YLadder,
  94.    CLadder& B2Ladder,
  95.    CLadder& Y2Ladder, CMSPeak *Peaks)
  96. {
  97.     // create hitlist.  note that this is deleted in the copy operator
  98.     HitInfo.reset(new CMSHitInfo[Hits]);
  99.     // increment thru hithist
  100.     int iHitInfo(0); 
  101.     int Which = Peaks->GetWhich(Charge);
  103.     // scan thru each ladder
  104.     if(Charge >= kConsiderMult) {
  105. RecordMatchesScan(BLadder, iHitInfo, Peaks, Which);
  106. RecordMatchesScan(YLadder, iHitInfo, Peaks, Which);
  107. RecordMatchesScan(B2Ladder, iHitInfo, Peaks, Which);
  108. RecordMatchesScan(Y2Ladder, iHitInfo, Peaks, Which);
  109.     }
  110.     else {
  111. RecordMatchesScan(BLadder, iHitInfo, Peaks, Which);
  112. RecordMatchesScan(YLadder, iHitInfo, Peaks, Which);
  113.     }
  115. }
  117. // return number of hits above threshold
  118. int CMSHit::GetHits(double Threshold, int MaxI)
  119. {
  120.     int i, retval(0);
  122.     for(i = 0; i < Hits; i++)
  123.         if(SetHitInfo(i).GetIntensity() > MaxI*Threshold)
  124.             retval++;
  125.     return retval;
  126. }
  130. /////////////////////////////////////////////////////////////////////////////
  131. //
  132. //  Comparison functions used in CMSPeak for sorting
  133. //
  135. // compares m/z.  Lower m/z first in sort.
  136. struct CMZICompare {
  137.     bool operator() (CMZI x, CMZI y)
  138.     {
  139. if(x.MZ < y.MZ) return true;
  140. return false;
  141.     }
  142. };
  144. // compares m/z.  High m/z first in sort.
  145. struct CMZICompareHigh {
  146.     bool operator() (CMZI x, CMZI y)
  147.     {
  148. if(x.MZ > y.MZ) return true;
  149. return false;
  150.     }
  151. };
  153. // compares intensity.  Higher intensities first in sort.
  154. struct CMZICompareIntensity {
  155.     bool operator() (CMZI x, CMZI y)
  156.     {
  157. if(x.Intensity > y.Intensity) return true;
  158. return false;
  159.     }
  160. };
  164. /////////////////////////////////////////////////////////////////////////////
  165. //
  166. //  CMSPeak::
  167. //
  168. //  Class used to hold spectra and convert to mass ladders
  169. //
  172. void CMSPeak::xCMSPeak(void)
  173. {
  174.     AAMap = AA.GetMap();
  175.     Sorted[0] = Sorted[1] = false;
  176.     MZI[MSCULLED1] = 0;
  177.     MZI[MSCULLED2] = 0;
  178.     MZI[MSORIGINAL] = 0;
  179.     MZI[MSTOPHITS] = 0;
  180.     Used[MSCULLED1] = 0;
  181.     Used[MSCULLED2] = 0;
  182.     Used[MSORIGINAL] = 0;    
  183.     Used[MSTOPHITS] = 0;    
  184.     PlusOne = 0.8L;
  185.     ComputedCharge = eChargeUnknown;
  186.     Error = eMSHitError_none;
  187. }
  189. CMSPeak::CMSPeak(void)
  190. { 
  191.     xCMSPeak(); 
  192. }
  194. CMSPeak::CMSPeak(int HitListSizeIn): HitListSize(HitListSizeIn)
  195. {  
  196.     xCMSPeak();
  197. }
  199. CMSPeak::~CMSPeak(void)
  200. {  
  201.     delete [] MZI[MSCULLED1];
  202.     delete [] Used[MSCULLED1];
  203.     delete [] MZI[MSCULLED2];
  204.     delete [] Used[MSCULLED2];
  205.     delete [] MZI[MSORIGINAL];
  206.     delete [] Used[MSTOPHITS];
  207.     delete [] MZI[MSTOPHITS];
  208.     int iCharges;
  209.     for(iCharges = 0; iCharges < GetNumCharges(); iCharges++)
  210. delete [] HitList[iCharges];
  211. }
  213. // add hit to hitlist.  returns true if successful
  214. bool CMSPeak::AddHit(CMSHit& in, CMSHit *& out)
  215. {
  216.     out = 0;
  217.     // initialize index using charge
  218.     int Index = in.GetCharge() - Charges[0];
  219.     // check to see if hitlist is full
  220.     if(HitListIndex[Index] >= HitListSize) {
  221. // if less or equal hits than recorded min, don't bother
  222. if(in.GetHits() <= LastHitNum[Index]) return false;  
  223. int i, min(HitList[Index][0].GetHits()); //the minimum number of hits in the list
  224. int minpos(0);     // the position of min
  225. // find the minimum in the hitslist
  226. for(i = 1; i < HitListSize; i++) {
  227.     if(min > HitList[Index][i].GetHits()) {
  228. min = HitList[Index][i].GetHits();
  229. minpos = i;
  230.     }
  231. }
  232. // keep record of min
  233. LastHitNum[Index] = min;
  234. // replace in list
  235. HitList[Index][minpos] = in;
  236. out =  & (HitList[Index][minpos]);
  237. return true;
  238.     } 
  239.     else {
  240. // add to end of list
  241. HitList[Index][HitListIndex[Index]] = in;
  242. out = & (HitList[Index][HitListIndex[Index]]);
  243. HitListIndex[Index]++;
  244. return true;
  245.     }
  246. }
  249. void CMSPeak::AddTotalMass(int massin, int tolin)
  250. {
  251.     TotalMass = massin;
  252.     tol = tolin;
  253. }
  256. void CMSPeak::Sort(int Which)
  257. {
  258.     if(Which < MSORIGINAL || Which >=  MSNUMDATA) return;
  259.     sort(MZI[Which], MZI[Which] + Num[Which], CMZICompare());
  260.     Sorted[Which] = true;
  261. }
  263. bool CMSPeak::Contains(int value, int Which)
  264. {
  265.     CMZI precursor;
  266.     precursor.MZ = value -  tol; // /2;
  267.     CMZI *p = lower_bound(MZI[Which], MZI[Which] + Num[Which], precursor, CMZICompare());
  268.     if(p == MZI[Which] + Num[Which]) return false;
  269.     if(p->MZ < value + tol/*/2*/) return true;
  270.     return false;
  271. }
  274. bool CMSPeak::ContainsFast(int value, int Which)
  275. {
  277.     int x, l, r;
  278.     
  279.     l = 0;
  280.     r = Num[Which] - 1;
  282.     while(l <= r) {
  283.         x = (l + r)/2;
  284.         if (MZI[Which][x].MZ < value - tol) 
  285.     l = x + 1;
  286.         else if (MZI[Which][x].MZ > value + tol)
  287.     r = x - 1;
  288. else return true;
  289.     } 
  290.     
  291.     if (x < Num[Which] - 1 && 
  292. MZI[Which][x+1].MZ < value + tol && MZI[Which][x+1].MZ > value - tol) 
  293. return true;
  294.     return false;
  295. }
  297. // compare assuming all lists are sorted
  298. // the intensity array holds the intensity if there is a match to the ladder
  299. // returns total number of matches, which may be more than is recorded in the ladder due to overlap
  300. int CMSPeak::CompareSorted(CLadder& Ladder, int Which, TIntensity* Intensity)
  301. {
  302.     unsigned i(0), j(0);
  303.     int retval(0);
  304.     if(Ladder.size() == 0 ||  Num[Which] == 0) return 0;
  306.     do {
  307. if (MZI[Which][j].MZ < Ladder[i] - tol) {
  308.     j++;
  309.     if(j >= Num[Which]) break;
  310.     continue;
  311. }
  312. else if(MZI[Which][j].MZ > Ladder[i] + tol) {
  313.     i++;
  314.     if(i >= Ladder.size()) break;
  315.     continue;
  316. }
  317. else {
  318.     if(Used[Which][j] == 0) {
  319. Used[Which][j] = 1;
  320. Ladder.GetHit()[i] = Ladder.GetHit()[i] + 1;
  321.     }
  322.     retval++;
  323.     // record the intensity if requested, used for auto adjust
  324.     if(Intensity) {
  325. *(Intensity->get() + i) = MZI[Which][j].Intensity;
  326.     }
  327.     j++;
  328.     if(j >= Num[Which]) break;
  329.     i++;
  330.     if(i >= Ladder.size()) break;
  331. }
  332.     } while(true);
  333.     return retval;
  334. }
  336. int CMSPeak::Compare(CLadder& Ladder, int Which)
  337. {
  338.     unsigned i;
  339.     int retval(0);
  340.     for(i = 0; i < Ladder.size(); i++) {
  341. if(ContainsFast(Ladder[i], Which)) {
  342.     retval++;
  343.     Ladder.GetHit()[i] = Ladder.GetHit()[i] + 1;
  344. }
  345.     }
  346.     return retval;
  347. }
  349. // compares the top peaks to the ladder.  returns immediately when finding a hit
  350. bool CMSPeak::CompareTop(CLadder& Ladder)
  351. {
  352.     unsigned i;
  353.     for(i = 0; i < Num[MSTOPHITS]; i++) {
  354. if(Ladder.ContainsFast(MZI[MSTOPHITS][i].MZ, tol)) return true;
  355.     }
  356.     return false;
  357. }
  360. // read in an asn.1 spectrum and initialize peak values
  361. int CMSPeak::Read(CMSSpectrum& Spectrum, double MSMSTolerance)
  362. {
  363.     try {
  364. int Scale = Spectrum.GetScale();
  365. TotalMass = Spectrum.GetPrecursormz()*MSSCALE/Scale;
  366. SetTolerance(MSMSTolerance);
  367. Charge = *(Spectrum.GetCharge().begin());
  368. Num[MSORIGINAL] = 0;   
  369.     
  370. const CMSSpectrum::TMz& Mz = Spectrum.GetMz();
  371. const CMSSpectrum::TAbundance& Abundance = Spectrum.GetAbundance();
  372. MZI[MSORIGINAL] = new CMZI [Mz.size()];
  373. Num[MSORIGINAL] = Mz.size();
  375. int i;
  376. for(i = 0; i < Num[MSORIGINAL]; i++) {
  377.     MZI[MSORIGINAL][i].MZ = Mz[i]*MSSCALE/Scale;
  378.     MZI[MSORIGINAL][i].Intensity = Abundance[i]/Scale;
  379. }
  380. Sort(MSORIGINAL);
  381.     } catch (NCBI_NS_STD::exception& e) {
  382. ERR_POST(Info << "Exception in CMSPeak::Read: " << e.what());
  383. throw;
  384.     }
  386.     return 0;
  387. }
  389. // return the number of peaks in a range
  390. int CMSPeak::CountRange(double StartFraction, double StopFraction)
  391. {
  392.     CMZI Start, Stop;
  393.     int Precursor = static_cast <int> (TotalMass/Charge + tol/2.0);
  394.     Start.MZ = static_cast <int> (StartFraction * Precursor);
  395.     Stop.MZ = static_cast <int> (StopFraction * Precursor);
  396.     CMZI *LoHit = lower_bound(MZI[MSORIGINAL], MZI[MSORIGINAL] +
  397.       Num[MSORIGINAL], Start, CMZICompare());
  398.     CMZI *HiHit = upper_bound(MZI[MSORIGINAL], MZI[MSORIGINAL] + 
  399.       Num[MSORIGINAL], Stop, CMZICompare());
  401.     return HiHit - LoHit;
  402. }
  404. // return the percentage of peaks below and at the precursor
  405. int CMSPeak::PercentBelow(void)
  406. {
  407.     CMZI precursor;
  408.     precursor.MZ = static_cast <int> (TotalMass/Charge + tol/2.0);
  409.     CMZI *Hit = upper_bound(MZI[MSORIGINAL], MZI[MSORIGINAL] + Num[MSORIGINAL],
  410.     precursor, CMZICompare());
  411.     Hit++;  // go above the peak
  412.     if(Hit >= MZI[MSORIGINAL] + Num[MSORIGINAL]) return Num[MSORIGINAL];
  413.     return Hit-MZI[MSORIGINAL];
  414. }
  416. // is the data charge +1?
  417. bool CMSPeak::IsPlus1(double PercentBelowIn)
  418. {
  419.     if(PercentBelowIn/Num[MSORIGINAL] > PlusOne) return true;
  420.     return false;
  421. }
  423. // takes the ratio, low/high, of two ranges in the spectrum
  424. double CMSPeak::RangeRatio(double Start, double Middle, double Stop)
  425. {
  426.     int Lo = CountRange(Start, Middle);
  427.     int Hi = CountRange(Middle, Stop);
  428.     return (double)Lo/Hi;
  429. }
  431. // calculate the charge
  432. void CMSPeak::SetComputedCharge(int MaxCharge)
  433. {
  434.     if(IsPlus1(PercentBelow())) {
  435. ComputedCharge = eCharge1;
  436. Charges[0] = 1;
  437. NumCharges = 1;
  438.     }
  439.     else {
  440. #if 0
  441. NumCharges = 1;
  442. if(RangeRatio(0.5, 1.0, 1.5) > 1.25 ) {
  443.     Charges[0] = 2;
  444.     ComputedCharge = eCharge2;
  445. }
  446. else if(RangeRatio(1.0, 1.5, 2.0) > 1.25) {
  447.     Charges[0] = 3;
  448.     ComputedCharge = eCharge3;
  449. }
  450. else {  // assume +4
  451.     Charges[0] = 4;
  452.     ComputedCharge = eCharge4;
  453. }
  454. #endif
  455. ComputedCharge = eChargeNot1; 
  456. Charges[0] = 2;
  457. Charges[1] = 3;
  458. NumCharges = 2;
  459.     }
  460. }
  462. // initializes arrays used to track hits
  463. void CMSPeak::InitHitList(void)
  464. {
  465.     int iCharges;
  466.     for(iCharges = 0; iCharges < GetNumCharges(); iCharges++) {
  467. LastHitNum[iCharges] = MSHITMIN - 1;
  468. HitListIndex[iCharges] = 0;
  469. HitList[iCharges] = new CMSHit[HitListSize];
  470. PeptidesExamined[iCharges] = 0;
  471.     }
  472. }
  474. void CMSPeak::xWrite(std::ostream& FileOut, CMZI *Temp, int Num)
  475. {
  476.     FileOut << (double)TotalMass/MSSCALE << " " << Charge << endl;
  478.     int i;
  479.     unsigned Intensity;
  480.     for(i = 0; i < Num; i++) {
  481.         Intensity = Temp[i].Intensity;
  482.         if(Intensity == 0.0) Intensity = 1;  // make Mascot happy
  483.         FileOut << (double)Temp[i].MZ/MSSCALE << " " << 
  484.     Intensity << endl;
  485.     }
  486. }
  488. void CMSPeak::Write(std::ostream& FileOut, CMSPeak::EFileType FileType, int Which)
  489. {
  490.     if(!FileOut || FileType != CMSPeak::eDTA) return;
  491.     xWrite(FileOut, MZI[Which], Num[Which]);
  492. }
  495. // counts the number of peaks above % of maximum peak
  497. int CMSPeak::AboveThresh(double Threshold, int Which)
  498. {
  499.     CMZI *MZISort = new CMZI[Num[Which]];
  500.     unsigned i;
  501.     for(i = 0; i < Num[Which]; i++) MZISort[i] = MZI[Which][i];
  503.     // sort so higher intensities first
  504.     sort(MZISort, MZISort+Num[Which], CMZICompareIntensity());
  506.     for(i = 1; i < Num[Which] &&
  507.     (double)MZISort[i].Intensity/MZISort[0].Intensity > Threshold; i++);
  508.     delete [] MZISort;
  509.     return i;
  510. }
  513. // Truncate the at the precursor if plus one and set charge to 1
  514. // if charge is erronously set to 1, set it to 2
  516. void CMSPeak::TruncatePlus1(void)
  517. {
  518.     int PercentBelowVal = PercentBelow();
  519.     if(IsPlus1(PercentBelowVal)) {
  520.         Num[MSORIGINAL] = PercentBelowVal;
  521.         Charge = 1;
  522.     }
  523.     else if(Charge == 1) Charge = 2;
  524. }
  526. // functions used in SmartCull
  527. // take out original peaks below a threshold.  assumes original peaks sorted by intensity
  528. void CMSPeak::CullBaseLine(double Threshold, CMZI *Temp, int& TempLen)
  529. {
  530.     unsigned iMZI;
  531.     for(iMZI = 0; iMZI < Num[MSORIGINAL] && MZI[MSORIGINAL][iMZI].Intensity > Threshold * MZI[MSORIGINAL][0].Intensity; iMZI++);
  532.     copy(&MZI[MSORIGINAL][0], &MZI[MSORIGINAL][iMZI], Temp);
  533.     TempLen = iMZI;
  534. }
  537. // cull precursors
  538. void CMSPeak::CullPrecursor(CMZI *Temp, int& TempLen, double Precursor)
  539. {
  540.     // chop out precursors
  541.     int iTemp(0), iMZI;
  542.     
  543.     for(iMZI = 0; iMZI < TempLen; iMZI++) { 
  544. if(Temp[iMZI].MZ > Precursor - tol && 
  545.     Temp[iMZI].MZ < Precursor + tol) continue;
  547. Temp[iTemp] = Temp[iMZI];
  548. iTemp++;
  549.     }
  550.     TempLen = iTemp;
  551. }
  553. // iterate thru peaks, deleting ones that pass the test
  554. void CMSPeak::CullIterate(CMZI *Temp, int& TempLen, TMZIbool FCN)
  555. {
  556.     if(!FCN) return;
  557.     int iTemp(0), iMZI, jMZI;
  558.     set <int> Deleted;
  559.     
  560.     // scan for isotopes, starting at highest peak
  561.     for(iMZI = 0; iMZI < TempLen - 1; iMZI++) { 
  562. if(Deleted.count(iMZI) != 0) continue;
  563. for(jMZI = iMZI + 1; jMZI < TempLen; jMZI++) { 
  564.     if(Deleted.count(jMZI) != 0) continue;
  565.     if((*FCN)(Temp[iMZI], Temp[jMZI], tol)) {
  566. Deleted.insert(jMZI);
  567. // Temp[iMZI].Intensity += Temp[jMZI].Intensity; // keep the intensity
  568.     }
  569. }
  570.     }
  571.     
  572.     for(iMZI = 0; iMZI < TempLen; iMZI++) {
  573. if(Deleted.count(iMZI) == 0) {
  574.     Temp[iTemp] = Temp[iMZI];
  575.     iTemp++;
  576. }
  577.     }
  578.     TempLen = iTemp;
  579. }
  581. // object for looking for isotopes.  true if isotope
  582. static bool FCompareIsotope(const CMZI& x, const CMZI& y, int tol)
  583. {
  584.     if(y.MZ < x.MZ + 2*MSSCALE + tol && y.MZ > x.MZ - tol) return true;
  585.     return false;
  586. }
  588. // cull isotopes using the Markey Method
  590. void CMSPeak::CullIsotope(CMZI *Temp, int& TempLen)
  591. {
  592.     CullIterate(Temp, TempLen, &FCompareIsotope);
  593. }
  596. // function for looking for H2O or NH3.  true if is.
  597. static bool FCompareH2ONH3(const CMZI& x, const CMZI& y, int tol)
  598. {
  599.     if((y.MZ > x.MZ - 18*MSSCALE - tol && y.MZ < x.MZ - 18*MSSCALE + tol ) ||
  600.        (y.MZ > x.MZ - 17*MSSCALE - tol && y.MZ < x.MZ - 17*MSSCALE + tol))
  601. return true;
  602.     return false;
  603. }
  605. // cull peaks that are water or ammonia loss
  606. // note that this only culls the water or ammonia loss if these peaks have a lesser
  607. // less intensity
  608. void CMSPeak::CullH20NH3(CMZI *Temp, int& TempLen)
  609. {
  610.     // need to start at top of m/z ladder, work way down.  sort
  611.     CullIterate(Temp, TempLen, &FCompareH2ONH3);
  612. }
  615. // use smartcull on all charges
  616. void CMSPeak::CullAll(double Threshold, int SingleWindow,
  617.       int DoubleWindow, int SingleHit, int DoubleHit)
  618. {    
  619.     int iCharges;
  620.     for(iCharges = 0; iCharges < GetNumCharges(); iCharges++)
  621. SmartCull(Threshold, GetCharges()[iCharges], SingleWindow,
  622.   DoubleWindow, SingleHit, DoubleHit);
  624.     // make the high intensity list
  625.     iCharges = GetNumCharges() - 1;
  626.     int Which = GetWhich(GetCharges()[iCharges]);
  627.     CMZI *Temp = new CMZI [Num[Which]]; // temporary holder
  628.     copy(MZI[Which], MZI[Which] + Num[Which], Temp);
  629.     sort(Temp, Temp + Num[Which], CMZICompareIntensity());
  631.     if(Num[Which] > MSNUMTOP)  Num[MSTOPHITS] = MSNUMTOP;
  632.     else  Num[MSTOPHITS] = Num[Which];
  634.     MZI[MSTOPHITS] = new CMZI [Num[MSTOPHITS]]; // holds highest hits
  635.     Used[MSTOPHITS] = new char [Num[MSTOPHITS]];
  636.     ClearUsed(MSTOPHITS);
  637.     Sorted[MSTOPHITS] = false;
  638.     copy(Temp, Temp + Num[MSTOPHITS], MZI[MSTOPHITS]);
  639.     Sort(MSTOPHITS);
  641.     delete [] Temp;
  642. }
  644. // check to see if TestMZ is Diff away from BigMZ
  645. bool CMSPeak::IsAtMZ(int BigMZ, int TestMZ, int Diff, int tol)
  646. {
  647.     if(TestMZ < BigMZ - Diff*MSSCALE + tol && 
  648.        TestMZ > BigMZ - Diff*MSSCALE - tol)
  649. return true;
  650.     return false;
  651. }
  653. // see if TestMZ can be associated with BigMZ, e.g. water loss, etc.
  654. bool CMSPeak::IsMajorPeak(int BigMZ, int TestMZ, int tol)
  655. {
  656.     if (IsAtMZ(BigMZ, TestMZ, 1, tol) || 
  657. IsAtMZ(BigMZ, TestMZ, 16, tol) ||
  658. IsAtMZ(BigMZ, TestMZ, 17, tol) ||
  659. IsAtMZ(BigMZ, TestMZ, 18, tol)) return true;
  660.     return false;
  661. }
  663. // recursively culls the peaks
  664. void CMSPeak::SmartCull(double Threshold, int Charge, int SingleWindow,
  665. int DoubleWindow, int SingleNum, int DoubleNum)
  666. {
  667.     //       sort(MZI[MSORIGINAL], MZI[MSORIGINAL] + Num[MSORIGINAL], CMZICompare());
  668.     sort(MZI[MSORIGINAL], MZI[MSORIGINAL] + Num[MSORIGINAL], CMZICompareIntensity());
  669.     Sorted[MSORIGINAL] = false;
  670.     double Precursor;
  671.     int Which = GetWhich(Charge);
  672.     Precursor = GetMass()/(double)Charge;
  674.     int iMZI = 0;  // starting point
  676.     int TempLen(0);
  677.     CMZI *Temp = new CMZI [Num[MSORIGINAL]]; // temporary holder
  679.     // prep the data
  680.     CullBaseLine(Threshold, Temp, TempLen);
  681. #ifdef DEBUG_PEAKS1
  682.     {
  683. sort(Temp, Temp+TempLen , CMZICompare());
  684. ofstream FileOut("afterbase.dta");
  685. xWrite(FileOut, Temp, TempLen);
  686. sort(Temp, Temp+TempLen , CMZICompareIntensity());
  687.     }
  688. #endif
  689.     CullPrecursor(Temp, TempLen, Precursor);
  690. #ifdef DEBUG_PEAKS1
  691.     {
  692. sort(Temp, Temp+TempLen , CMZICompare());
  693. ofstream FileOut("afterprecurse.dta");
  694. xWrite(FileOut, Temp, TempLen);
  695. sort(Temp, Temp+TempLen , CMZICompareIntensity());
  696.     }
  697. #endif
  698.     CullIsotope(Temp, TempLen);
  699. #ifdef DEBUG_PEAKS1
  700.     {
  701. sort(Temp, Temp+TempLen , CMZICompare());
  702. ofstream FileOut("afteriso.dta");
  703. xWrite(FileOut, Temp, TempLen);
  704. sort(Temp, Temp+TempLen , CMZICompareIntensity());
  705.     }
  706. #endif
  707.     // h20 and nh3 loss ignored until this cut can be studied
  708. #if 0
  709.     sort(Temp, Temp + TempLen, CMZICompareHigh());
  710.     CullH20NH3(Temp, TempLen);
  711.     sort(Temp, Temp + TempLen, CMZICompareIntensity());
  712. #endif
  714.     // now do the recursion
  715.     int iTemp(0), jMZI;
  716.     set <int> Deleted;
  717.     map <int, int> MajorPeak; // maps a peak to it's "major peak"
  718.     int HitsAllowed; // the number of hits allowed in a window
  719.     int Window; // the m/z range used to limit the number of peaks
  720.     int HitCount;  // number of nearby peaks
  721.    
  722.     // scan for isotopes, starting at highest peak
  723.     for(iMZI = 0; iMZI < TempLen - 1; iMZI++) { 
  724. if(Deleted.count(iMZI) != 0) continue;
  725. HitCount = 0;
  726. if(Charge <  kConsiderMult || Temp[iMZI].MZ > Precursor) {
  727.     // if charge 1 region, allow fewer peaks
  728.     Window = SingleWindow; //27;
  729.     HitsAllowed = SingleNum;
  730. }
  731. else {
  732.     // otherwise allow a greater number
  733.     Window = DoubleWindow; //14;
  734.     HitsAllowed = DoubleNum;
  735. }
  736. // go over smaller peaks
  737. set <int> Considered;
  738. for(jMZI = iMZI + 1; jMZI < TempLen; jMZI++) { 
  739.     if(Deleted.count(jMZI) != 0) continue;
  740.     if(Temp[jMZI].MZ < Temp[iMZI].MZ + Window*MSSCALE + tol &&
  741.        Temp[jMZI].MZ > Temp[iMZI].MZ - Window*MSSCALE - tol) {
  742. if(IsMajorPeak(Temp[iMZI].MZ, Temp[jMZI].MZ, tol)) {
  743.     // link little peak to big peak
  744.     MajorPeak[jMZI] = iMZI;
  745.     // ignore for deletion
  746.     continue;
  747. }
  748. // if linked to a major peak, skip
  749. if(MajorPeak.find(jMZI) != MajorPeak.end()) {
  750.     if(Considered.find(jMZI) != Considered.end())
  751. continue;
  752.     Considered.insert(jMZI);
  753.     HitCount++;
  754.     if(HitCount <= HitsAllowed)  continue;
  755. }
  756. // if the first peak, skip
  757. HitCount++;
  758. if(HitCount <= HitsAllowed)  continue;
  759. Deleted.insert(jMZI);
  760.     }
  761. }
  762.     }
  763.     
  764.     // delete the culled peaks
  765.     for(iMZI = 0; iMZI < TempLen; iMZI++) {
  766. if(Deleted.count(iMZI) == 0) {
  767.     Temp[iTemp] = Temp[iMZI];
  768.     iTemp++;
  769. }
  770.     }
  771.     TempLen = iTemp;
  774.     // make the array of culled peaks
  775.     if(MZI[Which]) delete [] MZI[Which];
  776.     if(Used[Which]) delete [] Used[Which];
  777.     Num[Which] = TempLen;
  778.     MZI[Which] = new CMZI [TempLen];
  779.     Used[Which] = new char [TempLen];
  780.     ClearUsed(Which);
  781.     copy(Temp, Temp+TempLen, MZI[Which]);
  782.     Sort(Which);
  783. #ifdef DEBUG_PEAKS1
  784.     {
  785. ofstream FileOut("aftercull.dta");
  786. xWrite(FileOut, MZI[Which], Num[Which]);
  787.     }
  788. #endif
  790.     delete [] Temp;
  791. }
  793. // return the lowest culled peak and the highest culled peak less than the
  794. // +1 precursor mass
  795. void CMSPeak::HighLow(int& High, int& Low, int& NumPeaks, int PrecursorMass,
  796.       int Charge, double Threshold, int& NumLo, int& NumHi)
  797. {
  798.     int Which = GetWhich(Charge);
  800.     if(!Sorted[Which]) Sort(Which);
  801.     if(Num[Which] < 2) {
  802. High = Low = -1;
  803. NumPeaks = NumLo = NumHi = 0;
  804. return;
  805.     }
  807.     Low = PrecursorMass;
  808.     High = 0;
  809.     NumPeaks = 0;
  810.     NumLo = 0;
  811.     NumHi = 0;
  813.     int MaxI = GetMaxI(Which);
  815.     unsigned iMZI;
  816.     for(iMZI = 0; iMZI < Num[Which]; iMZI++) {
  817. if(MZI[Which][iMZI].Intensity > Threshold*MaxI &&
  818.    MZI[Which][iMZI].MZ <= PrecursorMass) {
  819.     if(MZI[Which][iMZI].MZ > High) {
  820. High = MZI[Which][iMZI].MZ;
  821.     }
  822.     if(MZI[Which][iMZI].MZ < Low) {
  823. Low = MZI[Which][iMZI].MZ;
  824.     }
  825.     NumPeaks++;
  826.     if(MZI[Which][iMZI].MZ < PrecursorMass/2.0) NumLo++;
  827.     else NumHi++;
  828. }
  829.     }
  830. }
  832. // count the number of putative AA intervals.
  833. // Nodup true means don't count the peaks twice
  834. int CMSPeak::CountAAIntervals(CMassArray& MassArray, bool Nodup, int Which)
  835. {
  836.     unsigned ipeaks, low;
  837.     int i;
  838.     int PeakCount(0);
  839.     if(!Sorted[Which]) return -1;
  840.     const int *IntMassArray = MassArray.GetIntMass();
  842.     // list <int> intensity;
  843.     // int maxintensity = 0;
  845.     for(ipeaks = 0 ; ipeaks < Num[Which] - 1; ipeaks++) {
  846. // if(maxintensity < MZI[ipeaks].Intensity) maxintensity = MZI[ipeaks].Intensity;
  847. low = ipeaks;
  848. low++;
  849. // if(maxintensity < MZI[low].Intensity) maxintensity = MZI[low].Intensity;
  850. for(; low < Num[Which]; low++) {
  851.     for(i = 0; i < kNumUniqueAA; i++) {
  852. if(IntMassArray[i] == 0) continue;  // skip gaps, etc.
  853. if(MZI[Which][low].MZ- MZI[Which][ipeaks].MZ < IntMassArray[i] + tol/2.0 &&
  854.    MZI[Which][low].MZ - MZI[Which][ipeaks].MZ > IntMassArray[i] - tol/2.0 ) {   
  855.     PeakCount++;
  856.     // intensity.push_back(MZI[ipeaks].Intensity);
  857.     if(Nodup) goto newpeak;
  858.     else goto oldpeak;
  859. }
  860.     }
  861. oldpeak:
  862.     i = i;
  863. }
  864.     newpeak:
  865. i = i;
  866.     }
  868.     return PeakCount;
  869. }
  872. // return maximum intensity value
  873. int CMSPeak::GetMaxI(int Which)
  874. {
  875.     unsigned Intensity(0), i;
  876.     for(i = 0; i < Num[Which]; i++) {
  877.         if(Intensity < MZI[Which][i].Intensity)
  878.     Intensity = MZI[Which][i].Intensity;
  879.     }
  880.     return Intensity;
  881. }
  885. /////////////////////////////////////////////////////////////////////////////
  886. //
  887. //  CMSPeakSet::
  888. //
  889. //  Class used to hold sets of CMSPeak and access them quickly
  890. //
  893. CMSPeakSet::~CMSPeakSet() 
  894. {
  895.     while(!PeakSet.empty()) {
  896. delete *PeakSet.begin();
  897. PeakSet.pop_front();
  898.     }
  899. }
  901. // compares m/z.  Lower m/z first in sort.
  902. struct CMassPeakCompareHi {
  903.     bool operator() (TMassPeak x, TMassPeak y)
  904.     {
  905. if(x.Mass + x.Peptol < y.Mass + y.Peptol) return true;
  906. return false;
  907.     }
  908. };
  910. void CMSPeakSet::SortPeaks(int Peptol)
  911. {
  912.     int iCharges;
  913.     CMSPeak* Peaks;
  914.     TPeakSet::iterator iPeakSet;
  915.     int CalcMass; // the calculated mass
  916.     TMassPeak temp;
  917.     int ptol; // charge corrected mass tolerance
  919.     // first sort
  920.     for(iPeakSet = GetPeaks().begin();
  921. iPeakSet != GetPeaks().end();
  922. iPeakSet++ ) {
  923. Peaks = *iPeakSet;
  924. // skip empty spectra
  925. if(Peaks->GetError() == eMSHitError_notenuffpeaks) continue;
  927. // loop thru possible charges
  928. for(iCharges = 0; iCharges < Peaks->GetNumCharges(); iCharges++) {
  929.     // correction for incorrect charge determination.
  930.     // see 12/13/02 notebook, pg. 135
  931.     ptol = Peaks->GetCharges()[iCharges] * Peptol;
  932.     CalcMass = static_cast <int> ((Peaks->GetMass() +
  933.    Peaks->GetCharge()*kProton*MSSCALE) * 
  934.   Peaks->GetCharges()[iCharges]/(double)(Peaks->GetCharge()) - 
  935.   Peaks->GetCharges()[iCharges]*kProton*MSSCALE);
  936.     temp.Mass = CalcMass;
  937.     temp.Peptol = ptol;
  938.     temp.Charge = Peaks->GetCharges()[iCharges];
  939.     temp.Peak = Peaks;
  940.     // order by upper bound
  941.     MassMap.insert(pair <const int, TMassPeak>(temp.Mass + temp.Peptol, temp)); 
  942. }
  943.     } 
  945.     // then create static array
  947.     ArraySize = MassMap.size();
  948.     MassPeak.reset(new TMassPeak[ArraySize]);
  950.     TMassPeakMap::iterator iMassMap;
  951.     int i(0);
  952.     for(iMassMap = MassMap.begin(); iMassMap != MassMap.end(); iMassMap++, i++) {
  953. GetMassPeak(i).Mass =
  954.  iMassMap->second.Mass;
  955. GetMassPeak(i).Peptol = iMassMap->second.Peptol;
  956. GetMassPeak(i).Peak = iMassMap->second.Peak;
  957. GetMassPeak(i).Charge = iMassMap->second.Charge;
  958.     }
  960.     MassMap.clear();
  962. }
  965. // Get the first index into the sorted array where the mass
  966. // + tolerance is >= the given calculated mass. 
  967. TMassPeak *CMSPeakSet::GetIndexLo(int Mass)
  968. {
  969.     TMassPeak temp;
  970.     TMassPeak *retval;
  971.     temp.Mass = Mass;
  972.     temp.Peptol = 0;
  973.     // look for first spectrum whose upper mass value + tolerance exceeds calculated mass
  974.     retval = lower_bound(MassPeak.get(), MassPeak.get() + ArraySize, temp, 
  975.  CMassPeakCompareHi());
  976.     return retval;
  977. }
  979. // get peak for sorted list by index into list
  980. CMSPeak *CMSPeakSet::GetPeak(int Index)
  981. {
  982.     if(Index < 0 || Index >= ArraySize) return 0;
  983.     return GetMassPeak(Index).Peak;
  984. }