← Back to branch summary

~ubuntu-branches/ubuntu/wily/openms/wily

« back to all changes in this revision

Viewing changes to source/ANALYSIS/OPENSWATH/OPENSWATHALGO/ALGO/MRMScoring.C

  • Committer: Package Import Robot
  • Author(s): Filippo Rusconi
  • Date: 2013-12-20 11:30:16 UTC
  • mfrom: (5.1.2 sid)
  • Revision ID: package-import@ubuntu.com-20131220113016-wre5g9bteeheq6he
Tags: 1.11.1-3
* remove version number from libbost development package names;
* ensure that AUTHORS is correctly shipped in all packages.

Show diffs side-by-side

added added

removed removed

Lines of Context:
source/ANALYSIS/OPENSWATH/OPENSWATHALGO/ALGO/MRMScoring.C
 
1
// --------------------------------------------------------------------------
 
2
//                   OpenMS -- Open-Source Mass Spectrometry
 
3
// --------------------------------------------------------------------------
 
4
// Copyright The OpenMS Team -- Eberhard Karls University Tuebingen,
 
5
// ETH Zurich, and Freie Universitaet Berlin 2002-2013.
 
6
//
 
7
// This software is released under a three-clause BSD license:
 
8
//  * Redistributions of source code must retain the above copyright
 
9
//    notice, this list of conditions and the following disclaimer.
 
10
//  * Redistributions in binary form must reproduce the above copyright
 
11
//    notice, this list of conditions and the following disclaimer in the
 
12
//    documentation and/or other materials provided with the distribution.
 
13
//  * Neither the name of any author or any participating institution
 
14
//    may be used to endorse or promote products derived from this software
 
15
//    without specific prior written permission.
 
16
// For a full list of authors, refer to the file AUTHORS.
 
17
// --------------------------------------------------------------------------
 
18
// THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
 
19
// AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
 
20
// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
 
21
// ARE DISCLAIMED. IN NO EVENT SHALL ANY OF THE AUTHORS OR THE CONTRIBUTING
 
22
// INSTITUTIONS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
 
23
// EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
 
24
// PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS;
 
25
// OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY,
 
26
// WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR
 
27
// OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF
 
28
// ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
 
29
//
 
30
// --------------------------------------------------------------------------
 
31
// $Maintainer: Hannes Roest $
 
32
// $Authors: Hannes Roest $
 
33
// --------------------------------------------------------------------------
 
34
 
 
35
//#define MRMSCORING_TESTING
 
36
#include <algorithm>
 
37
#include <algorithm>
 
38
#include <iterator>
 
39
#include <iostream>
 
40
 
 
41
#include "OpenMS/ANALYSIS/OPENSWATH/OPENSWATHALGO/ALGO/MRMScoring.h"
 
42
#include "OpenMS/ANALYSIS/OPENSWATH/OPENSWATHALGO/ALGO/Scoring.h"
 
43
#include "OpenMS/ANALYSIS/OPENSWATH/OPENSWATHALGO/ALGO/StatsHelpers.h"
 
44
 
 
45
#ifdef OPENMS_ASSERTIONS
 
46
#define OPENMS_PRECONDITION(condition, message)\
 
47
        if (!(condition))\
 
48
    { throw std::runtime_error(message); }
 
49
#else
 
50
#define OPENMS_PRECONDITION(condition, message)
 
51
#endif
 
52
 
 
53
namespace OpenSwath
 
54
{
 
55
 
 
56
  const MRMScoring::XCorrMatrixType & MRMScoring::getXCorrMatrix() const
 
57
  {
 
58
    return xcorr_matrix_;
 
59
  }
 
60
 
 
61
  void MRMScoring::initializeXCorrMatrix(OpenSwath::IMRMFeature * mrmfeature, OpenSwath::ITransitionGroup * transition_group, bool normalize)
 
62
  {
 
63
    std::vector<double> intensityi, intensityj;
 
64
    xcorr_matrix_.resize(transition_group->size());
 
65
    for (std::size_t i = 0; i < transition_group->size(); i++)
 
66
    {
 
67
      String native_id = transition_group->getNativeIDs()[i];
 
68
      FeatureType fi = mrmfeature->getFeature(native_id);
 
69
      xcorr_matrix_[i].resize(transition_group->size());
 
70
      intensityi.clear();
 
71
      fi->getIntensity(intensityi);
 
72
      for (std::size_t j = i; j < transition_group->size(); j++)
 
73
      {
 
74
        String native_id2 = transition_group->getNativeIDs()[j];
 
75
        FeatureType fj = mrmfeature->getFeature(native_id2);
 
76
        intensityj.clear();
 
77
        fj->getIntensity(intensityj);
 
78
        //std::cout << " = Computing crosscorrelation " << i << " / "  << j << " or " << native_id << " vs " << native_id2 << std::endl;
 
79
        if (normalize)
 
80
        {
 
81
          xcorr_matrix_[i][j] = Scoring::normalizedCalcxcorr(intensityi, intensityj, boost::numeric_cast<int>(intensityi.size()), 1);
 
82
        }
 
83
        else
 
84
        {
 
85
          throw "not implemented";
 
86
        }
 
87
 
 
88
      }
 
89
    }
 
90
  }
 
91
 
 
92
  // see /IMSB/users/reiterl/bin/code/biognosys/trunk/libs/mrm_libs/MRM_pgroup.pm
 
93
  // _calc_xcorr_coelution_score
 
94
  //
 
95
  //   for each i,j get xcorr_matrix array => find max of the crosscorrelation
 
96
  //   store the delta to the retention time
 
97
  // return $deltascore_mean + $deltascore_stdev
 
98
  double MRMScoring::calcXcorrCoelutionScore()
 
99
  {
 
100
    OPENMS_PRECONDITION(xcorr_matrix_.size() > 1, "Expect cross-correlation matrix of at least 2x2");
 
101
 
 
102
    std::vector<int> deltas;
 
103
    for (std::size_t i = 0; i < xcorr_matrix_.size(); i++)
 
104
    {
 
105
      for (std::size_t  j = i; j < xcorr_matrix_.size(); j++)
 
106
      {
 
107
        // first is the X value (RT), should be an int
 
108
        deltas.push_back(std::abs(Scoring::xcorrArrayGetMaxPeak(xcorr_matrix_[i][j])->first));
 
109
#ifdef MRMSCORING_TESTING
 
110
        std::cout << "&&_xcoel append " << std::abs(Scoring::xcorrArrayGetMaxPeak(xcorr_matrix_[i][j])->first) << std::endl;
 
111
#endif
 
112
      }
 
113
    }
 
114
 
 
115
    OpenSwath::mean_and_stddev msc;
 
116
    msc = std::for_each(deltas.begin(), deltas.end(), msc);
 
117
    double deltas_mean = msc.mean();
 
118
    double deltas_stdv = msc.sample_stddev();
 
119
    //double deltas_mean = gsl_stats_int_mean(&deltas[0], 1, deltas.size());
 
120
    //double deltas_stdv = gsl_stats_int_sd(&deltas[0], 1, deltas.size());
 
121
 
 
122
    double xcorr_coelution_score = deltas_mean + deltas_stdv;
 
123
    return xcorr_coelution_score;
 
124
  }
 
125
 
 
126
  double MRMScoring::calcXcorrCoelutionScore_weighted(
 
127
    const std::vector<double> & normalized_library_intensity)
 
128
  {
 
129
    OPENMS_PRECONDITION(xcorr_matrix_.size() > 1, "Expect cross-correlation matrix of at least 2x2");
 
130
 
 
131
#ifdef MRMSCORING_TESTING
 
132
    double weights = 0;
 
133
#endif
 
134
    std::vector<double> deltas;
 
135
    for (std::size_t i = 0; i < xcorr_matrix_.size(); i++)
 
136
    {
 
137
      deltas.push_back(
 
138
        std::abs(Scoring::xcorrArrayGetMaxPeak(xcorr_matrix_[i][i])->first)
 
139
        * normalized_library_intensity[i]
 
140
        * normalized_library_intensity[i]);
 
141
#ifdef MRMSCORING_TESTING
 
142
      std::cout << "_xcoel_weighted " << i << " " << i << " " << Scoring::xcorrArrayGetMaxPeak(xcorr_matrix_[i][i])->first << " weight " <<
 
143
      normalized_library_intensity[i] * normalized_library_intensity[i] << std::endl;
 
144
      weights += normalized_library_intensity[i] * normalized_library_intensity[i];
 
145
#endif
 
146
      for (std::size_t j = i + 1; j < xcorr_matrix_.size(); j++)
 
147
      {
 
148
        // first is the X value (RT), should be an int
 
149
        deltas.push_back(
 
150
          std::abs(Scoring::xcorrArrayGetMaxPeak(xcorr_matrix_[i][j])->first)
 
151
          * normalized_library_intensity[i]
 
152
          * normalized_library_intensity[j] * 2);
 
153
#ifdef MRMSCORING_TESTING
 
154
        std::cout << "_xcoel_weighted " << i << " " << j << " " << Scoring::xcorrArrayGetMaxPeak(xcorr_matrix_[i][j])->first << " weight " <<
 
155
        normalized_library_intensity[i] * normalized_library_intensity[j] * 2 << std::endl;
 
156
        weights += normalized_library_intensity[i] * normalized_library_intensity[j];
 
157
#endif
 
158
 
 
159
      }
 
160
    }
 
161
 
 
162
#ifdef MRMSCORING_TESTING
 
163
    std::cout << " all weights sum " << weights << std::endl;
 
164
#endif
 
165
 
 
166
    /*
 
167
     double deltas_mean = gsl_stats_int_mean(&deltas[0],1,deltas.size() );
 
168
     double deltas_stdv = gsl_stats_int_sd(  &deltas[0],1,deltas.size() );
 
169
 
 
170
     double xcorr_coelution_score = deltas_mean + deltas_stdv;
 
171
     return xcorr_coelution_score;
 
172
     */
 
173
    return std::accumulate(deltas.begin(), deltas.end(), 0.0);
 
174
  }
 
175
 
 
176
  // see /IMSB/users/reiterl/bin/code/biognosys/trunk/libs/mrm_libs/MRM_pgroup.pm
 
177
  // _calc_xcorr_shape_score
 
178
  //
 
179
  //   for each i,j get xcorr_matrix array => find max of the crosscorrelation
 
180
  //   calculate whether the maximal crosscorrelation coincides with the maximal intensity
 
181
  ///
 
182
  double MRMScoring::calcXcorrShape_score()
 
183
  {
 
184
    OPENMS_PRECONDITION(xcorr_matrix_.size() > 1, "Expect cross-correlation matrix of at least 2x2");
 
185
 
 
186
    std::vector<double> intensities;
 
187
    for (std::size_t i = 0; i < xcorr_matrix_.size(); i++)
 
188
    {
 
189
      for (std::size_t j = i; j < xcorr_matrix_.size(); j++)
 
190
      {
 
191
        // second is the Y value (intensity)
 
192
        intensities.push_back(Scoring::xcorrArrayGetMaxPeak(xcorr_matrix_[i][j])->second);
 
193
      }
 
194
    }
 
195
    OpenSwath::mean_and_stddev msc;
 
196
    msc = std::for_each(intensities.begin(), intensities.end(), msc);
 
197
    return msc.mean();
 
198
  }
 
199
 
 
200
  double MRMScoring::calcXcorrShape_score_weighted(
 
201
    const std::vector<double> & normalized_library_intensity)
 
202
  {
 
203
    OPENMS_PRECONDITION(xcorr_matrix_.size() > 1, "Expect cross-correlation matrix of at least 2x2");
 
204
 
 
205
    // TODO (hroest) : check implementation
 
206
    //         see _calc_weighted_xcorr_shape_score in MRM_pgroup.pm
 
207
    //         -- they only multiply up the intensity once
 
208
    std::vector<double> intensities;
 
209
    for (std::size_t i = 0; i < xcorr_matrix_.size(); i++)
 
210
    {
 
211
      intensities.push_back(
 
212
        Scoring::xcorrArrayGetMaxPeak(xcorr_matrix_[i][i])->second
 
213
        * normalized_library_intensity[i]
 
214
        * normalized_library_intensity[i]);
 
215
#ifdef MRMSCORING_TESTING
 
216
      std::cout << "_xcorr_weighted " << i << " " << i << " " << Scoring::xcorrArrayGetMaxPeak(xcorr_matrix_[i][i])->second << " weight " <<
 
217
      normalized_library_intensity[i] * normalized_library_intensity[i] << std::endl;
 
218
#endif
 
219
      for (std::size_t j = i + 1; j < xcorr_matrix_.size(); j++)
 
220
      {
 
221
        intensities.push_back(
 
222
          Scoring::xcorrArrayGetMaxPeak(xcorr_matrix_[i][j])->second
 
223
          * normalized_library_intensity[i]
 
224
          * normalized_library_intensity[j] * 2);
 
225
#ifdef MRMSCORING_TESTING
 
226
        std::cout << "_xcorr_weighted " << i << " " << j << " " << Scoring::xcorrArrayGetMaxPeak(xcorr_matrix_[i][j])->second << " weight " <<
 
227
        normalized_library_intensity[i] * normalized_library_intensity[j] * 2 << std::endl;
 
228
#endif
 
229
      }
 
230
    }
 
231
    return std::accumulate(intensities.begin(), intensities.end(), 0.0);
 
232
  }
 
233
 
 
234
  void MRMScoring::calcLibraryScore(OpenSwath::IMRMFeature * mrmfeature, const std::vector<TransitionType> & transitions,
 
235
                                           double & correlation, double & rmsd, double & manhattan, double & dotprod)
 
236
  {
 
237
    std::vector<double> library_intensity;
 
238
    std::vector<double> experimental_intensity;
 
239
    String native_id;
 
240
 
 
241
    for (std::size_t k = 0; k < transitions.size(); k++)
 
242
    {
 
243
      native_id = transitions[k].getNativeID();
 
244
      double intensity = transitions[k].getLibraryIntensity();
 
245
      // the library intensity should never be below zero
 
246
      if (intensity < 0.0)
 
247
      {
 
248
        intensity = 0.0;
 
249
      }
 
250
      experimental_intensity.push_back(mrmfeature->getFeature(native_id)->getIntensity());
 
251
      library_intensity.push_back(intensity);
 
252
    }
 
253
 
 
254
    OPENMS_PRECONDITION(library_intensity.size() == experimental_intensity.size(), "Both vectors need to have the same size");
 
255
 
 
256
#ifdef MRMSCORING_TESTING
 
257
    for (std::size_t k = 0; k < transitions.size(); k++)
 
258
    {
 
259
      native_id = transitions[k].getNativeID();
 
260
      std::cout << native_id << " Lib vs exp " << library_intensity[k] << " " << experimental_intensity[k] << std::endl;
 
261
    }
 
262
#endif
 
263
 
 
264
    manhattan = OpenSwath::manhattanScoring(experimental_intensity, library_intensity);
 
265
    dotprod = OpenSwath::dotprodScoring(experimental_intensity, library_intensity);
 
266
 
 
267
 
 
268
    Scoring::normalize_sum(&experimental_intensity[0], boost::numeric_cast<unsigned int>(transitions.size()));
 
269
    Scoring::normalize_sum(&library_intensity[0], boost::numeric_cast<unsigned int>(transitions.size()));
 
270
 
 
271
    rmsd = Scoring::RMSD(&experimental_intensity[0], &library_intensity[0], boost::numeric_cast<unsigned int>(transitions.size()));
 
272
    correlation = OpenSwath::cor_pearson(experimental_intensity.begin(), experimental_intensity.end(), library_intensity.begin());
 
273
 
 
274
 
 
275
    //double c0, c1, cov00, cov01, cov11, sumsq;
 
276
    //int ret = gsl_fit_linear(normx, 1, normy, 1, transitions.size(), &c0, &c1, &cov00, &cov01, &cov11, &sumsq);
 
277
    if (boost::math::isnan(correlation))
 
278
      correlation = -1.0;
 
279
  }
 
280
 
 
281
  double MRMScoring::calcRTScore(const PeptideType & peptide, double normalized_experimental_rt)
 
282
  {
 
283
    double expected_rt;
 
284
    expected_rt = peptide.rt;
 
285
 
 
286
    if (expected_rt <= -1000)
 
287
    {
 
288
      return 0;
 
289
    }
 
290
 
 
291
    // use the transformed experimental retention time and then take the difference.
 
292
    double rt_score = std::fabs(normalized_experimental_rt - expected_rt);
 
293
    return rt_score;
 
294
  }
 
295
 
 
296
  double MRMScoring::calcSNScore(OpenSwath::IMRMFeature * mrmfeature, std::vector<OpenSwath::ISignalToNoisePtr> & signal_noise_estimators)
 
297
  {
 
298
    OPENMS_PRECONDITION(signal_noise_estimators.size() > 1, "Input S/N estimators needs to be larger than 1");
 
299
 
 
300
    double sn_score = 0;
 
301
    if (signal_noise_estimators.size() == 0) 
 
302
    {
 
303
      return 0;
 
304
    }
 
305
 
 
306
    for (std::size_t k = 0; k < signal_noise_estimators.size(); k++)
 
307
    {
 
308
      sn_score += signal_noise_estimators[k]->getValueAtRT(mrmfeature->getRT());
 
309
    }
 
310
    return sn_score / signal_noise_estimators.size();
 
311
  }
 
312
 
 
313
}