1 /*
2 * ===========================================================================
3 *
4 *                            PUBLIC DOMAIN NOTICE
5 *               National Center for Biotechnology Information
6 *
7 *  This software/database is a "United States Government Work" under the
8 *  terms of the United States Copyright Act.  It was written as part of
9 *  the author's official duties as a United States Government employee and
10 *  thus cannot be copyrighted.  This software/database is freely available
11 *  to the public for use. The National Library of Medicine and the U.S.
12 *  Government have not placed any restriction on its use or reproduction.
13 *
14 *  Although all reasonable efforts have been taken to ensure the accuracy
15 *  and reliability of the software and data, the NLM and the U.S.
16 *  Government do not and cannot warrant the performance or results that
17 *  may be obtained by using this software or data. The NLM and the U.S.
18 *  Government disclaim all warranties, express or implied, including
19 *  warranties of performance, merchantability or fitness for any particular
20 *  purpose.
21 *
22 *  Please cite the authors in any work or product based on this material.
23 *
24 * ===========================================================================
25 *
26 * Author:  Lewis Y. Geer
27 *
28 * File Description:
29 *    code to do mass spec scoring
30 *
31 * ===========================================================================
32 */
33 
34 //
35 //
36 //
37 // Please keep this file clear of NCBI toolkit dependent code as it
38 // is meant to be used externally from the toolkit.
39 //
40 // NCBI dependent code can go into omssascore.cpp
41 //
42 //
43 //
44 
45 // precompiled headers #include found in msscore.cpp
46 #if 0
47     #include <ncbi_pch.hpp>
48 #endif
49 
50 #ifdef HAVE_LIMITS
51     #include <limits>
52 #endif
53 #include <math.h>
54 #include "msscore.hpp"
55 
56 #ifdef USING_OMSSA_SCOPE
57 USING_OMSSA_SCOPE
58 #endif
59 
60 // define log of gamma fcn if not defined already
61 #ifndef MSLNGAMMA
62     #define MSLNGAMMA lgamma
63 #endif
64 
65 // define error function
66 #ifndef MSERF
67     #define MSERF erf
68 #endif
69 
70 
IsForwardSeries(EMSIonSeries Series)71 const int CMSBasicMatchedPeak::IsForwardSeries(EMSIonSeries Series)
72 {
73     if ( Series == eMSIonTypeA ||
74          Series == eMSIonTypeB ||
75          Series == eMSIonTypeC)
76         return 1;
77     if ( Series == eMSIonTypeX ||
78          Series == eMSIonTypeY ||
79          Series == eMSIonTypeZ)
80         return 0;
81     return -1;
82 }
83 
CMSMatchedPeak(void)84 CMSMatchedPeak::CMSMatchedPeak(void):
85 MassTolerance(-1),
86 ExpIons(-1),
87 MatchType(eMSMatchTypeUnknown)
88 {
89 }
90 
Assign(CMSMatchedPeak * in)91 void CMSMatchedPeak::Assign(CMSMatchedPeak *in)
92 {
93     if (!in || in == this) return;
94     SetCharge() = in->GetCharge();
95     SetExpIons() = in->GetExpIons();
96     SetIntensity() = in->GetIntensity();
97     SetIonSeries() = in->GetIonSeries();
98     SetMassTolerance() = in->GetMassTolerance();
99     SetMatchType() = in->GetMatchType();
100     SetMZ() = in->GetMZ();
101     SetNumber() = in->GetNumber();
102 }
103 
Assign(CMSBasicMatchedPeak * in)104 void CMSMatchedPeak::Assign(CMSBasicMatchedPeak *in)
105 {
106     if (!in || in == this) return;
107     SetCharge() = in->GetCharge();
108     SetIntensity() = in->GetIntensity();
109     SetIonSeries() = in->GetIonSeries();
110     SetMZ() = in->GetMZ();
111     SetNumber() = in->GetNumber();
112 
113     SetMassTolerance() = -1;
114     SetExpIons() = -1;
115     SetMatchType() = eMSMatchTypeUnknown;
116 }
117 
118 
CMSMatchedPeakSet(void)119 CMSMatchedPeakSet::CMSMatchedPeakSet(void)
120 {
121 }
122 
~CMSMatchedPeakSet()123 CMSMatchedPeakSet::~CMSMatchedPeakSet()
124 {
125     DeleteMatchedPeakSet();
126 }
127 
CreateMatchedPeakSet(int SizeIn)128 void CMSMatchedPeakSet::CreateMatchedPeakSet(int SizeIn)
129 {
130     DeleteMatchedPeakSet();
131     int i;
132     for (i = 0; i < SizeIn; ++i)
133         SetMatchedPeakSet().push_back(new CMSMatchedPeak);
134 }
135 
DeleteMatchedPeakSet(void)136 void CMSMatchedPeakSet::DeleteMatchedPeakSet(void)
137 {
138     while (!GetMatchedPeakSet().empty()) {
139         delete SetMatchedPeakSet().back();
140         SetMatchedPeakSet().pop_back();
141     }
142 }
143 
144 void
Compare(CMSMatchedPeakSet * Other,bool SameDirection)145 CMSMatchedPeakSet::Compare(CMSMatchedPeakSet *Other, bool SameDirection)
146 {
147     if(!Other ||
148        (GetMatchedPeakSet().size() != Other->GetMatchedPeakSet().size()) ) return;
149 
150     unsigned i, j;
151 
152     for(i = 0; i < GetMatchedPeakSet().size(); ++i) {
153         if(SameDirection) j = i;
154         else j = Other->GetSize() - i - 1;
155         // don't do the comparison if there are no peaks to compare!
156         if(j >= Other->GetMatchedPeakSet().size()) break;
157 
158         if ((GetMatchedPeakSet()[i]->GetMatchType() == eMSMatchTypeIndependent ||
159              GetMatchedPeakSet()[i]->GetMatchType() == eMSMatchTypeSemiIndependent) &&
160             (Other->GetMatchedPeakSet()[j]->GetMatchType() == eMSMatchTypeIndependent ||
161              Other->GetMatchedPeakSet()[j]->GetMatchType() == eMSMatchTypeSemiIndependent)) {
162             Other->GetMatchedPeakSet()[j]->SetMatchType() = eMSMatchTypeDependent;
163         }
164     }
165 }
166 
CMSMatchedPeakSetMap(void)167 CMSMatchedPeakSetMap::CMSMatchedPeakSetMap(void)
168 {
169 }
170 
~CMSMatchedPeakSetMap()171 CMSMatchedPeakSetMap::~CMSMatchedPeakSetMap()
172 {
173     TIonSeriesMatchMap::iterator i;
174     for ( i = SetMatchMap().begin(); i != SetMatchMap().end(); ++i)
175         delete i->second;
176     SetMatchMap().clear();
177 }
178 
179 CMSMatchedPeakSet *
CreateSeries(TMSCharge Charge,TMSIonSeries Series,int Size,int Maxproductions)180     CMSMatchedPeakSetMap::CreateSeries(
181         TMSCharge Charge,
182         TMSIonSeries Series,
183         int Size,
184         int Maxproductions)
185 {
186     CMSMatchedPeakSet *retval;
187     int Key = ChargeSeries2Key(Charge, Series);
188     unsigned Realsize = min(Size, Maxproductions);
189 
190     // find old matches
191     if (GetMatchMap().find(Key) != GetMatchMap().end()) {
192         // if size is same, don't bother reallocating
193         if (SetMatchMap()[Key]->SetMatchedPeakSet().size() == Realsize) {
194             return SetMatchMap()[Key];
195         }
196         // otherwise delete and reallocate
197         else {
198             delete SetMatchMap()[Key];
199             SetMatchMap().erase(Key);
200         }
201     }
202 
203     retval = new CMSMatchedPeakSet;
204     if (retval) {
205         retval->CreateMatchedPeakSet(Realsize);
206         retval->SetSize() = Size;
207         SetMatchMap().insert(std::pair <int, CMSMatchedPeakSet * > (Key, retval));
208     }
209 
210     return retval;
211 }
212 
SetSeries(TMSCharge Charge,TMSIonSeries Series)213 CMSMatchedPeakSet * CMSMatchedPeakSetMap::SetSeries(TMSCharge Charge, TMSIonSeries Series)
214 {
215     CMSMatchedPeakSet *retval;
216     const int Key = ChargeSeries2Key(Charge, Series);
217     if (GetMatchMap().find(Key) != GetMatchMap().end())
218         retval = SetMatchMap()[Key];
219     else
220         retval = 0;
221 
222     return retval;
223 }
224 
225 //typedef std::map <int, CMSMatchedPeakSet *> TIonSeriesMatchMap;
226 
227 const int
ChargeSeries2Key(TMSCharge Charge,TMSIonSeries Series)228 CMSMatchedPeakSetMap::ChargeSeries2Key(TMSCharge Charge, TMSIonSeries Series)
229 {
230     return eMSIonTypeMax * Charge + Series;
231 }
232 
233 const TMSCharge
Key2Charge(int Key)234 CMSMatchedPeakSetMap::Key2Charge(int Key)
235 {
236     return static_cast <TMSCharge> (Key/eMSIonTypeMax);
237 }
238 
239 const TMSIonSeries
Key2Series(int Key)240 CMSMatchedPeakSetMap::Key2Series(int Key)
241 {
242     return static_cast <TMSIonSeries> (Key - Key2Charge(Key)*eMSIonTypeMax);
243 }
244 
245 
246 
CMSSpectrumMatch(void)247 CMSSpectrumMatch::CMSSpectrumMatch(void):
248 HitInfo(0),
249 Hits(0),
250 ExpMass(-1),
251 TheoreticalMass(-1),
252 Sum(0),
253 M(0),
254 N(0)
255 {
256 }
257 
~CMSSpectrumMatch()258 CMSSpectrumMatch::~CMSSpectrumMatch()
259 {
260     delete [] HitInfo;
261 }
262 
CreateHitInfo(void)263 void CMSSpectrumMatch::CreateHitInfo(void)
264 {
265     delete [] HitInfo;
266     HitInfo = 0;
267     if (GetHits() > 0)
268         HitInfo = new CMSBasicMatchedPeak[GetHits()];
269 }
270 
Find(TMSNumber Number,TMSCharge ChargeIn,TMSIonSeries Series)271 CMSBasicMatchedPeak * CMSSpectrumMatch::Find(
272                                             TMSNumber Number,
273                                             TMSCharge ChargeIn,
274                                             TMSIonSeries Series)
275 {
276     int i;
277     for (i = 0; i < GetHits(); ++i) {
278         if (GetHitInfo(i).GetNumber() == Number &&
279             GetHitInfo(i).GetCharge() == ChargeIn &&
280             GetHitInfo(i).GetIonSeries() == Series)
281             return &SetHitInfo(i);
282     }
283     return 0;
284 }
285 
286 
FillMatchedPeaks(TMSCharge ChargeIn,TMSIonSeries Series,unsigned Size,TMSIntensity MinIntensity,bool Skipb1,EMSTerminalBias TerminalIon,int Maxproductions,string & Sequence,bool NoProline)287 void CMSSpectrumMatch::FillMatchedPeaks(
288                                        TMSCharge ChargeIn,
289                                        TMSIonSeries Series,
290                                        unsigned Size,
291                                        TMSIntensity MinIntensity,
292                                        bool Skipb1,
293                                        EMSTerminalBias TerminalIon,
294                                        int Maxproductions
295 //#if 0
296                                        ,
297                                        string &Sequence,
298                                        bool NoProline
299 //#endif
300                                        )
301 {
302     // create a new match set
303     CMSMatchedPeakSet *MatchPeakSet =
304     SetIonSeriesMatchMap().CreateSeries(ChargeIn, Series, Size, Maxproductions);
305 
306     unsigned i;
307 
308     // iterate over series and look for matching hits
309     for (i = 0; i < MatchPeakSet->GetMatchedPeakSet().size(); ++i) {
310         CMSBasicMatchedPeak * FoundPeak = Find(i, ChargeIn, Series);
311         if (FoundPeak && FoundPeak->GetIntensity() > MinIntensity) {
312             // copy hit to match
313             MatchPeakSet->SetMatchedPeakSet()[i]->Assign(FoundPeak);
314             // set as untyped hit
315             MatchPeakSet->SetMatchedPeakSet()[i]->SetMatchType() =
316             eMSMatchTypeNotTyped;
317 
318             // check for terminal match
319             // n-term
320             if (TerminalIon == eMSNTerminalBias || TerminalIon == eMSBothTerminalBias) {
321                 if ( (i==0 && CMSBasicMatchedPeak::IsForwardSeries(static_cast <EMSIonSeries> (Series))) ||
322                      (i==Size-1 && !CMSBasicMatchedPeak::IsForwardSeries(static_cast <EMSIonSeries> (Series)))) {
323                     MatchPeakSet->SetMatchedPeakSet()[i]->SetMatchType() = eMSMatchTypeTerminus;
324                 }
325             }
326             // c-term
327             else if (TerminalIon == eMSCTerminalBias || TerminalIon == eMSBothTerminalBias) {
328                 if ( (i==0 && !CMSBasicMatchedPeak::IsForwardSeries(static_cast <EMSIonSeries> (Series))) ||
329                      (i==Size-1 && CMSBasicMatchedPeak::IsForwardSeries(static_cast <EMSIonSeries> (Series)))) {
330                     MatchPeakSet->SetMatchedPeakSet()[i]->SetMatchType() = eMSMatchTypeTerminus;
331                 }
332             }
333         }
334         else {
335             MatchPeakSet->SetMatchedPeakSet()[i]->SetCharge() = ChargeIn;
336             MatchPeakSet->SetMatchedPeakSet()[i]->SetIonSeries() = Series;
337             MatchPeakSet->SetMatchedPeakSet()[i]->SetMatchType() = eMSMatchTypeUnknown;
338             MatchPeakSet->SetMatchedPeakSet()[i]->SetMZ() = 0;
339             MatchPeakSet->SetMatchedPeakSet()[i]->SetIntensity() = 0;
340             MatchPeakSet->SetMatchedPeakSet()[i]->SetNumber() = i;
341         }
342         // if b1 should not be searched, mark it as so.
343         if (i == 0 && Skipb1 && CMSBasicMatchedPeak::IsForwardSeries(static_cast <EMSIonSeries> (Series)))
344             MatchPeakSet->SetMatchedPeakSet()[i]->SetMatchType() = eMSMatchTypeNoSearch;
345 //#if 0
346         // apply rule that no
347         if (NoProline &&
348             CMSBasicMatchedPeak::IsForwardSeries(static_cast <EMSIonSeries> (Series)) &&
349             i+1 < Sequence.size() &&
350             Sequence[i+1] == 'P') {
351             MatchPeakSet->SetMatchedPeakSet()[i]->SetMatchType() = eMSMatchTypeNoSearch;
352         }
353 
354         if (NoProline &&
355             !CMSBasicMatchedPeak::IsForwardSeries(static_cast <EMSIonSeries> (Series)) &&
356             i < Sequence.size() &&
357             Sequence[Sequence.size() - i - 1] == 'P') {
358             MatchPeakSet->SetMatchedPeakSet()[i]->SetMatchType() = eMSMatchTypeNoSearch;
359         }
360 //#endif
361 
362     }
363 
364 
365     // iterate over array, look for unset m/z.  keep track of index of last set m/z
366     // if unset, iterate until first set, then fill in mz
367 
368     // find first -1
369     unsigned j;
370     for (i = 0; i < MatchPeakSet->GetMatchedPeakSet().size(); ++i) {
371         if (MatchPeakSet->GetMatchedPeakSet()[i]->GetMZ() != 0) break;
372     }
373     // if first filled hit is not the first one, fill in array
374     for (j = 0; j < i && i < MatchPeakSet->GetMatchedPeakSet().size(); ++j) {
375         MatchPeakSet->SetMatchedPeakSet()[j]->SetMZ() =
376         (MatchPeakSet->GetMatchedPeakSet()[i]->GetMZ()/ (i+1)) * (j+1);
377     }
378     // now fill out sections in the center
379     unsigned LastIndex(i);
380     for (; i < MatchPeakSet->GetMatchedPeakSet().size(); ++i) {
381         if (MatchPeakSet->GetMatchedPeakSet()[i]->GetMZ() != 0) {
382             if (i != LastIndex) {
383                 // fill in the blank spaces
384                 for (j = LastIndex+1; j < i; j++) {
385                     MatchPeakSet->SetMatchedPeakSet()[j]->SetMZ() =
386                     MatchPeakSet->GetMatchedPeakSet()[LastIndex]->GetMZ() +
387                     ((MatchPeakSet->GetMatchedPeakSet()[i]->GetMZ() -
388                       MatchPeakSet->GetMatchedPeakSet()[LastIndex]->GetMZ()) /
389                      (i - LastIndex )) * (j - LastIndex) ;
390                 }
391             }
392             LastIndex = i;
393         }
394     }
395 
396     // now fill out last section
397     if (MatchPeakSet->GetMatchedPeakSet()[MatchPeakSet->GetMatchedPeakSet().size()-1]->GetMZ() == 0) {
398 
399         int StartIndex;
400         TMSMZ StartMass;
401         TMSMZ MassIncrement;
402         // if completely empty, force the array to be completely filled out
403         if (MatchPeakSet->GetMatchedPeakSet()[0]->GetMZ() == 0) {
404             MassIncrement = (GetExpMass()/ChargeIn)/(Size+1);
405             StartIndex = 0;
406             StartMass = (GetExpMass()/ChargeIn)/(Size+1);
407         }
408         else {
409             MassIncrement = ((GetExpMass()/ChargeIn) - MatchPeakSet->GetMatchedPeakSet()[LastIndex]->GetMZ()) /
410                             (Size - LastIndex);
411             StartIndex = LastIndex + 1;
412             StartMass = MatchPeakSet->GetMatchedPeakSet()[LastIndex]->GetMZ() +
413                         MassIncrement;
414         }
415 
416         for (j = StartIndex; j < MatchPeakSet->GetMatchedPeakSet().size(); ++j) {
417             MatchPeakSet->SetMatchedPeakSet()[j]->SetMZ() =
418             StartMass + MassIncrement * (j - StartIndex);
419         }
420     }
421 
422 
423     // first determine if the first peak is independent.  This simplifies the loop
424     if (MatchPeakSet->GetMatchedPeakSet()[0]->GetMatchType() == eMSMatchTypeNotTyped)
425         MatchPeakSet->SetMatchedPeakSet()[0]->SetMatchType() = eMSMatchTypeIndependent;
426 
427     // now characterize all of the matches depending on the preceeding match
428     for (i = 1; i < MatchPeakSet->GetMatchedPeakSet().size(); ++i) {
429         // only type matched peaks
430         if (MatchPeakSet->GetMatchedPeakSet()[i]->GetMatchType() == eMSMatchTypeNotTyped) {
431             if (MatchPeakSet->GetMatchedPeakSet()[i-1]->GetMatchType() == eMSMatchTypeIndependent ||
432                 MatchPeakSet->GetMatchedPeakSet()[i-1]->GetMatchType() == eMSMatchTypeSemiIndependent) {
433                 MatchPeakSet->SetMatchedPeakSet()[i]->SetMatchType() = eMSMatchTypeDependent;
434             }
435             else if (MatchPeakSet->GetMatchedPeakSet()[i-1]->GetMatchType() == eMSMatchTypeNoSearch ||
436                      MatchPeakSet->GetMatchedPeakSet()[i-1]->GetMatchType() == eMSMatchTypeUnknown ||
437                      MatchPeakSet->GetMatchedPeakSet()[i-1]->GetMatchType() == eMSMatchTypeTerminus) {
438                 MatchPeakSet->SetMatchedPeakSet()[i]->SetMatchType() = eMSMatchTypeIndependent;
439             }
440             else if (MatchPeakSet->GetMatchedPeakSet()[i-1]->GetMatchType() == eMSMatchTypeDependent) {
441                 MatchPeakSet->SetMatchedPeakSet()[i]->SetMatchType() = eMSMatchTypeSemiIndependent;
442             }
443         }
444     }
445 
446     return;
447 }
448 
449 
450 
451 const double
CalcPoissonMean(double ProbTerminal,int NumTerminalMasses,double ProbDependent,int NumUniqueMasses,double ToleranceAdjust) const452 CMSSpectrumMatch::CalcPoissonMean(double ProbTerminal,
453                                   int NumTerminalMasses,
454                                   double ProbDependent,
455                                   int NumUniqueMasses,
456                                   double ToleranceAdjust
457                                   ) const
458 {
459     double Mean(0.0L);
460     // scan thru entire map
461     TIonSeriesMatchMap::const_iterator i;
462     for (i = GetIonSeriesMatchMap().GetMatchMap().begin();
463         i != GetIonSeriesMatchMap().GetMatchMap().end();
464         ++i) {
465         CMSMatchedPeakSet *PeakSet = i->second;
466         // scan thru each match set
467         TMatchedPeakSet::const_iterator j;
468         for (j = PeakSet->GetMatchedPeakSet().begin();
469             j != PeakSet->GetMatchedPeakSet().end();
470             ++j) {
471             CMSMatchedPeak *Peak = *j;
472 
473             // skip unsearched peaks
474             if (Peak->GetMatchType() == eMSMatchTypeNoSearch) continue;
475 
476             // calculate single poisson
477             // add it to total poisson
478             Mean += 2.0 * Peak->GetMassTolerance() * Peak->GetExpIons() * ToleranceAdjust;
479             _TRACE( "mean=" << Mean << " masstol=" << Peak->GetMassTolerance() << " expions=" <<  Peak->GetExpIons());
480 
481             // if a statistically dependent peak, add in dependent match probability
482             if (Peak->GetMatchType() == eMSMatchTypeDependent && NumUniqueMasses != 0) {
483                 Mean += ProbDependent / NumUniqueMasses;
484                 _TRACE( "mean=" << Mean << " probdep=" << ProbDependent << " numunique=" << NumUniqueMasses);
485             }
486 
487             // if a statistically biased terminal peak, add in biased match probability
488             if (Peak->GetMatchType() == eMSMatchTypeTerminus && NumTerminalMasses != 0) {
489                 Mean += ProbTerminal / NumTerminalMasses;
490             }
491         }
492     }
493     return Mean;
494 }
495 
CalcPoisson(double Mean,int i) const496 const double CMSSpectrumMatch::CalcPoisson(double Mean, int i) const
497 {
498     return exp(-Mean) * pow(Mean, i) / exp(MSLNGAMMA(i+1));
499 }
500 
CalcPoissonTopHit(double Mean,int i,double TopHitProb) const501 const double CMSSpectrumMatch::CalcPoissonTopHit(double Mean, int i, double TopHitProb) const
502 {
503     double retval;
504     retval =  exp(-Mean) * pow(Mean, i) / exp(MSLNGAMMA(i+1));
505     if (TopHitProb != 1.0L) retval *= 1.0L - pow((1.0L-TopHitProb), i);
506     return retval;
507 }
508 
CalcNormalTopHit(double Mean,double TopHitProb) const509 const double CMSSpectrumMatch::CalcNormalTopHit(double Mean, double TopHitProb) const
510 {
511     int i;
512     double retval(0.0L), before(-1.0L), increment;
513 
514     for (i = 1; i < 1000; i++) {
515         increment = CalcPoissonTopHit(Mean, i, TopHitProb);
516         // convergence hack -- on linux (at least) convergence test doesn't work
517         // for gcc release build
518         if (increment <= MSDOUBLELIMIT && i > 10) break;
519         retval += increment;
520         if (retval == before) break;  // convergence
521         before = retval;
522     }
523     return retval;
524 }
525 
CalcPvalue(double Mean,int HitsIn) const526 const double CMSSpectrumMatch::CalcPvalue(double Mean, int HitsIn) const
527 {
528     if (HitsIn <= 0) return 1.0L;
529 
530     int i;
531     double retval(0.0L), increment, beforeincrement(0.0L);
532 
533     for (i = 0; i < HitsIn; i++) {
534         increment = CalcPoisson(Mean, i);
535 #ifdef HAVE_LIMITS
536         // detect limit of increment if numeric_limits available
537         if(increment < retval * std::numeric_limits <double> ::epsilon()) break;
538 #endif
539         retval += increment;
540         if(retval == 1.0L) break;
541         beforeincrement = increment;
542     }
543 
544     // if retval is calculated to the maximum precision of double
545     // back off one increment to give a lower bound of the e-value
546     if(retval == 1.0L)
547         retval -= beforeincrement;
548 
549     retval = 1.0L - retval;
550     if (retval <= 0.0L) retval = MSDOUBLELIMIT;
551     return retval;
552 }
553 
CalcPvalueTopHit(double Mean,int HitsIn,double Normal,double TopHitProb) const554 const double CMSSpectrumMatch::CalcPvalueTopHit(double Mean, int HitsIn, double Normal, double TopHitProb) const
555 {
556     if (HitsIn <= 0) return 1.0L;
557 
558     int i;
559     double retval(0.0L), increment, beforeretval(-1.0), beforeincrement(0.0L);
560 
561     for (i = 1; i < HitsIn; i++) {
562         increment = CalcPoissonTopHit(Mean, i, TopHitProb);
563 #ifdef HAVE_LIMITS
564         // detect limit of increment if numeric_limits available
565         if(increment < retval * std::numeric_limits <double> ::epsilon()) break;
566 #endif
567         retval += increment;
568         // note that in gcc 3.4.2 on linux with -O9 -ffast-math, the following statement doesn't appear to work
569         // trying to fix this with retval - beforeretval < increment also does not work -- it breaks too soon!
570         if(retval == beforeretval) break;
571         beforeretval = retval;
572         beforeincrement = increment;
573     }
574 
575     // if retval and Normal are calculated to the maximum precision of double
576     // back off one increment to give a lower bound of the e-value
577     if(retval == Normal) retval -= beforeincrement;
578 
579     retval /= Normal;
580     retval = 1.0L - retval;
581     if (retval <= 0.0L) retval = MSDOUBLELIMIT;
582     return retval;
583 }
584 
CalcRankProb(void) const585 const double CMSSpectrumMatch::CalcRankProb(void) const
586 {
587     double Average, StdDev, Perf(1.0L);
588     if (GetM() != 0.0) {
589         Average = (GetM() *(GetN()+1))/2.0;
590         StdDev = sqrt(Average * (GetN() - GetM()) / 6.0);
591         if ( StdDev != 0.0)
592             Perf = 0.5*(1.0 + MSERF((GetSum() - Average)/(StdDev*sqrt(2.0))));
593     }
594     return Perf;
595 }
596 
597 
598 const TMSMZ
GetMeanDelta(void) const599 CMSSpectrumMatch::GetMeanDelta(void) const
600 {
601     TMSMZ Mean(0);
602 
603     int i;
604     for (i = 0; i < GetHits(); ++i) {
605         Mean += GetHitInfo(i).GetDelta();
606     }
607     return static_cast <int> (Mean/(static_cast <double> (GetHits())));
608 }
609 
610 
611 const TMSMZ
GetStdDevDelta(void) const612 CMSSpectrumMatch::GetStdDevDelta(void) const
613 {
614     TMSMZ Mean;
615     double StdDev(0.0L);
616 
617     Mean = GetMeanDelta();
618 
619     int i;
620      for (i = 0; i < GetHits(); ++i) {
621          StdDev += pow((double)GetHitInfo(i).GetDelta() - Mean, 2.0);
622      }
623      StdDev = pow((double)StdDev/GetHits(), 0.5);
624      return (TMSMZ)StdDev;
625 }
626 
627 
628 const TMSMZ
GetMaxDelta(void) const629 CMSSpectrumMatch::GetMaxDelta(void) const
630 {
631      TMSMZ Max(0);
632      int i;
633      for (i = 0; i < GetHits(); ++i) {
634          if(abs(GetHitInfo(i).GetDelta()) > Max)
635              Max = abs(GetHitInfo(i).GetDelta());
636      }
637      return Max;
638 }
639 
640