← Back to branch summary

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

« back to all changes in this revision

Viewing changes to source/TRANSFORMATIONS/FEATUREFINDER/ExtendedIsotopeModel.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/TRANSFORMATIONS/FEATUREFINDER/ExtendedIsotopeModel.C
1
 
// -*- mode: C++; tab-width: 2; -*-
2
 
// vi: set ts=2:
3
 
//
4
 
// --------------------------------------------------------------------------
5
 
//                   OpenMS Mass Spectrometry Framework
6
 
// --------------------------------------------------------------------------
7
 
//  Copyright (C) 2003-2011 -- Oliver Kohlbacher, Knut Reinert
8
 
//
9
 
//  This library is free software; you can redistribute it and/or
10
 
//  modify it under the terms of the GNU Lesser General Public
11
 
//  License as published by the Free Software Foundation; either
12
 
//  version 2.1 of the License, or (at your option) any later version.
13
 
//
14
 
//  This library is distributed in the hope that it will be useful,
15
 
//  but WITHOUT ANY WARRANTY; without even the implied warranty of
16
 
//  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
17
 
//  Lesser General Public License for more details.
18
 
//
19
 
//  You should have received a copy of the GNU Lesser General Public
20
 
//  License along with this library; if not, write to the Free Software
21
 
//  Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
 
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.
22
29
//
23
30
// --------------------------------------------------------------------------
24
31
// $Maintainer: Clemens Groepl $
34
41
 
35
42
namespace OpenMS
36
43
{
37
 
        ExtendedIsotopeModel::ExtendedIsotopeModel()
38
 
                : InterpolationModel(),
39
 
                        charge_(0),
40
 
                        monoisotopic_mz_(0.0)
41
 
        {
42
 
                setName(getProductName());
43
 
 
44
 
    defaults_.setValue("averagines:C",0.04443989f,"Number of C atoms per Dalton of mass.", StringList::create("advanced"));
45
 
                defaults_.setValue("averagines:H",0.06981572f,"Number of H atoms per Dalton of mass.", StringList::create("advanced"));
46
 
                defaults_.setValue("averagines:N",0.01221773f,"Number of N atoms per Dalton of mass.", StringList::create("advanced"));
47
 
                defaults_.setValue("averagines:O",0.01329399f,"Number of O atoms per Dalton of mass.", StringList::create("advanced"));
48
 
                defaults_.setValue("averagines:S",0.00037525f,"Number of S atoms per Dalton of mass.", StringList::create("advanced"));
49
 
 
50
 
                defaults_.setValue("isotope:trim_right_cutoff",0.001,"Cutoff in averagine distribution, trailing isotopes below this relative intensity are not considered.", StringList::create("advanced"));
51
 
                defaults_.setValue("isotope:maximum",100,"Maximum isotopic rank to be considered.", StringList::create("advanced"));
52
 
                defaults_.setValue("isotope:distance",1.000495,"Distance between consecutive isotopic peaks.", StringList::create("advanced"));
53
 
                defaults_.setValue("isotope:stdev",0.1,"Standard deviation of gaussian applied to the averagine isotopic pattern to simulate the inaccuracy of the mass spectrometer.", StringList::create("advanced"));
54
 
                defaults_.setValue("charge",1,"Charge state of the model.", StringList::create("advanced"));
55
 
                defaults_.setValue("isotope:monoisotopic_mz",1.0,"Monoisotopic m/z of the model.", StringList::create("advanced"));
56
 
 
57
 
                defaultsToParam_();
58
 
        }
59
 
 
60
 
        ExtendedIsotopeModel::ExtendedIsotopeModel(const ExtendedIsotopeModel& source)
61
 
                : InterpolationModel(source)
62
 
        {
63
 
                setParameters( source.getParameters() );
64
 
                updateMembers_();
65
 
        }
66
 
 
67
 
        ExtendedIsotopeModel::~ExtendedIsotopeModel()
68
 
        {
69
 
        }
70
 
 
71
 
        ExtendedIsotopeModel& ExtendedIsotopeModel::operator = (const ExtendedIsotopeModel& source)
72
 
        {
73
 
                if (&source == this) return *this;
74
 
 
75
 
                InterpolationModel::operator = (source);
76
 
                setParameters( source.getParameters() );
77
 
                updateMembers_();
78
 
 
79
 
                return *this;
80
 
        }
81
 
 
82
 
        void ExtendedIsotopeModel::setSamples()
83
 
        {
84
 
                // MAGIC alert, num stdev for smooth table for normal distribution
85
 
                CoordinateType normal_widening_num_stdev = 4.;
86
 
                // Actual width for values in the smooth table for normal distribution
87
 
                CoordinateType normal_widening_width = isotope_stdev_ * normal_widening_num_stdev;
88
 
 
89
 
                typedef std::vector < double > ContainerType;
90
 
                ContainerType isotopes_exact;
91
 
                CoordinateType mass = monoisotopic_mz_ * charge_;
92
 
 
93
 
                Int C_num = Int( 0.5 + mass * averagine_[C]);
94
 
                Int N_num = Int( 0.5 + mass * averagine_[N]);
95
 
                Int O_num = Int( 0.5 + mass * averagine_[O]);
96
 
                Int H_num = Int( 0.5 + mass * averagine_[H]);
97
 
                Int S_num = Int( 0.5 + mass * averagine_[S]);
98
 
 
99
 
                String form("");
100
 
                if (C_num) form.append("C").append(String(C_num));
101
 
                if (H_num) form.append("H").append(String(H_num));
102
 
                if (N_num) form.append("N").append(String(N_num));
103
 
                if (O_num) form.append("O").append(String(O_num));
104
 
                if (S_num) form.append("S").append(String(S_num));
105
 
 
106
 
                EmpiricalFormula formula(form);
107
 
                IsotopeDistribution isotope_distribution = formula.getIsotopeDistribution(max_isotope_);
108
 
                isotope_distribution.trimRight(trim_right_cutoff_);
109
 
                isotope_distribution.renormalize();
110
 
 
111
 
                // compute the average mass (-offset)
112
 
                for (IsotopeDistribution::iterator iter = isotope_distribution.begin(); iter != isotope_distribution.end(); ++iter)
113
 
                {
114
 
                        isotopes_exact.push_back(iter->second);
115
 
                }
116
 
 
117
 
                // "stretch" the averagine isotope distribution
118
 
                Size isotopes_exact_size = isotopes_exact.size();
119
 
                isotopes_exact.resize(Size( (isotopes_exact_size-1)
120
 
                                                                                                                                *isotope_distance_/interpolation_step_+1.6)); // round up a bit more
121
 
 
122
 
                for ( Size i = isotopes_exact_size-1; i; --i )
123
 
                        // we don't need to move the 0-th entry
124
 
                {
125
 
                        isotopes_exact [Size(CoordinateType(i)*
126
 
                                                                                                         isotope_distance_/interpolation_step_/charge_+0.5)]
127
 
                                =       isotopes_exact [ i ];
128
 
                        isotopes_exact [ i ] = 0;
129
 
                }
130
 
 
131
 
                // compute the normal distribution (to be added for widening the averagine isotope distribution)
132
 
                Math::BasicStatistics<> normal_widening_model;
133
 
                normal_widening_model.setSum  (1);
134
 
                normal_widening_model.setMean (0);
135
 
                normal_widening_model.setVariance(isotope_stdev_*isotope_stdev_);
136
 
                // fill a container with CoordinateType points
137
 
                ContainerType normal_widening_coordinate;
138
 
                for ( DoubleReal coord = -normal_widening_width;
139
 
                                        coord <= normal_widening_width;
140
 
                                        coord += interpolation_step_
141
 
                                )
142
 
                {
143
 
                        normal_widening_coordinate.push_back(coord);
144
 
                }
145
 
                // compute normal approximation at these CoordinateType points
146
 
                ContainerType normal_widening;
147
 
                normal_widening_model.normalApproximation(normal_widening,normal_widening_coordinate );
148
 
 
149
 
                // fill linear interpolation
150
 
                const ContainerType& left = isotopes_exact;
151
 
                const ContainerType& right = normal_widening;
152
 
                ContainerType& result = interpolation_.getData();
153
 
                result.clear ();
154
 
 
155
 
                Int rMax = std::min ( Int( left.size() + right.size() - 1 ), Int(2*normal_widening_width/interpolation_step_*max_isotope_+1) );
156
 
                result.resize ( rMax, 0 );
157
 
 
158
 
                // we loop backwards because then the small products tend to come first
159
 
                // (for better numerics)
160
 
                for ( SignedSize i = left.size() - 1; i >= 0; --i )
161
 
                {
162
 
      if (left[i]==0) continue;
163
 
                        for ( SignedSize j = std::min<SignedSize> ( rMax - i, right.size() ) - 1; j >= 0 ; --j )
164
 
                        {
165
 
                                result[i+j] += left[i] * right[j];
166
 
                        }
167
 
                }
168
 
 
169
 
                // set interpolation
170
 
                interpolation_.setMapping(interpolation_step_, normal_widening_width / interpolation_step_, monoisotopic_mz_ )  ;
171
 
 
172
 
                // scale data so that integral over distribution equals one
173
 
                // multiply sum by interpolation_step_ -> rectangular approximation of integral
174
 
                IntensityType factor = scaling_ / interpolation_step_ /
175
 
                        std::accumulate ( result.begin(), result.end(), IntensityType(0) );
176
 
 
177
 
                for ( ContainerType::iterator iter = result.begin(); iter != result.end(); ++iter)
178
 
                {
179
 
                        *iter *= factor;
180
 
                }
181
 
 
182
 
        }
183
 
 
184
 
        void ExtendedIsotopeModel::setOffset(CoordinateType offset)
185
 
        {
186
 
                DoubleReal diff = offset - getInterpolation().getOffset();
187
 
                monoisotopic_mz_ += diff;
188
 
 
189
 
                InterpolationModel::setOffset(offset);
190
 
 
191
 
                param_.setValue("isotope:monoisotopic_mz", monoisotopic_mz_);
192
 
        }
193
 
 
194
 
        ExtendedIsotopeModel::CoordinateType ExtendedIsotopeModel::getOffset()
195
 
        {
196
 
                return getInterpolation().getOffset();
197
 
        }
198
 
 
199
 
        UInt ExtendedIsotopeModel::getCharge()
200
 
        {
201
 
                return charge_;
202
 
        }
203
 
 
204
 
        ExtendedIsotopeModel::CoordinateType ExtendedIsotopeModel::getCenter() const
205
 
        {
206
 
                return monoisotopic_mz_;
207
 
        }
208
 
 
209
 
        void ExtendedIsotopeModel::updateMembers_()
210
 
        {
211
 
                InterpolationModel::updateMembers_();
212
 
 
213
 
                charge_ = param_.getValue("charge");
214
 
                isotope_stdev_ = param_.getValue("isotope:stdev");
215
 
                monoisotopic_mz_ = param_.getValue("isotope:monoisotopic_mz");
216
 
                max_isotope_ = param_.getValue("isotope:maximum");
217
 
                trim_right_cutoff_ = param_.getValue("isotope:trim_right_cutoff");
218
 
                isotope_distance_ = param_.getValue("isotope:distance");
219
 
 
220
 
                averagine_[C] = param_.getValue("averagines:C");
221
 
                averagine_[H] = param_.getValue("averagines:H");
222
 
                averagine_[N] = param_.getValue("averagines:N");
223
 
                averagine_[O] = param_.getValue("averagines:O");
224
 
                averagine_[S] = param_.getValue("averagines:S");
225
 
 
226
 
                setSamples();
227
 
        }
 
44
  ExtendedIsotopeModel::ExtendedIsotopeModel() :
 
45
    InterpolationModel(),
 
46
    charge_(0),
 
47
    monoisotopic_mz_(0.0)
 
48
  {
 
49
    setName(getProductName());
 
50
 
 
51
    defaults_.setValue("averagines:C", 0.04443989f, "Number of C atoms per Dalton of mass.", StringList::create("advanced"));
 
52
    defaults_.setValue("averagines:H", 0.06981572f, "Number of H atoms per Dalton of mass.", StringList::create("advanced"));
 
53
    defaults_.setValue("averagines:N", 0.01221773f, "Number of N atoms per Dalton of mass.", StringList::create("advanced"));
 
54
    defaults_.setValue("averagines:O", 0.01329399f, "Number of O atoms per Dalton of mass.", StringList::create("advanced"));
 
55
    defaults_.setValue("averagines:S", 0.00037525f, "Number of S atoms per Dalton of mass.", StringList::create("advanced"));
 
56
 
 
57
    defaults_.setValue("isotope:trim_right_cutoff", 0.001, "Cutoff in averagine distribution, trailing isotopes below this relative intensity are not considered.", StringList::create("advanced"));
 
58
    defaults_.setValue("isotope:maximum", 100, "Maximum isotopic rank to be considered.", StringList::create("advanced"));
 
59
    defaults_.setValue("isotope:distance", 1.000495, "Distance between consecutive isotopic peaks.", StringList::create("advanced"));
 
60
    defaults_.setValue("isotope:stdev", 0.1, "Standard deviation of gaussian applied to the averagine isotopic pattern to simulate the inaccuracy of the mass spectrometer.", StringList::create("advanced"));
 
61
    defaults_.setValue("charge", 1, "Charge state of the model.", StringList::create("advanced"));
 
62
    defaults_.setValue("isotope:monoisotopic_mz", 1.0, "Monoisotopic m/z of the model.", StringList::create("advanced"));
 
63
 
 
64
    defaultsToParam_();
 
65
  }
 
66
 
 
67
  ExtendedIsotopeModel::ExtendedIsotopeModel(const ExtendedIsotopeModel & source) :
 
68
    InterpolationModel(source)
 
69
  {
 
70
    setParameters(source.getParameters());
 
71
    updateMembers_();
 
72
  }
 
73
 
 
74
  ExtendedIsotopeModel::~ExtendedIsotopeModel()
 
75
  {
 
76
  }
 
77
 
 
78
  ExtendedIsotopeModel & ExtendedIsotopeModel::operator=(const ExtendedIsotopeModel & source)
 
79
  {
 
80
    if (&source == this)
 
81
      return *this;
 
82
 
 
83
    InterpolationModel::operator=(source);
 
84
    setParameters(source.getParameters());
 
85
    updateMembers_();
 
86
 
 
87
    return *this;
 
88
  }
 
89
 
 
90
  void ExtendedIsotopeModel::setSamples()
 
91
  {
 
92
    // MAGIC alert, num stdev for smooth table for normal distribution
 
93
    CoordinateType normal_widening_num_stdev = 4.;
 
94
    // Actual width for values in the smooth table for normal distribution
 
95
    CoordinateType normal_widening_width = isotope_stdev_ * normal_widening_num_stdev;
 
96
 
 
97
    typedef std::vector<double> ContainerType;
 
98
    ContainerType isotopes_exact;
 
99
    CoordinateType mass = monoisotopic_mz_ * charge_;
 
100
 
 
101
    Int C_num = Int(0.5 + mass * averagine_[C]);
 
102
    Int N_num = Int(0.5 + mass * averagine_[N]);
 
103
    Int O_num = Int(0.5 + mass * averagine_[O]);
 
104
    Int H_num = Int(0.5 + mass * averagine_[H]);
 
105
    Int S_num = Int(0.5 + mass * averagine_[S]);
 
106
 
 
107
    String form("");
 
108
    if (C_num)
 
109
      form.append("C").append(String(C_num));
 
110
    if (H_num)
 
111
      form.append("H").append(String(H_num));
 
112
    if (N_num)
 
113
      form.append("N").append(String(N_num));
 
114
    if (O_num)
 
115
      form.append("O").append(String(O_num));
 
116
    if (S_num)
 
117
      form.append("S").append(String(S_num));
 
118
 
 
119
    EmpiricalFormula formula(form);
 
120
    IsotopeDistribution isotope_distribution = formula.getIsotopeDistribution(max_isotope_);
 
121
    isotope_distribution.trimRight(trim_right_cutoff_);
 
122
    isotope_distribution.renormalize();
 
123
 
 
124
    // compute the average mass (-offset)
 
125
    for (IsotopeDistribution::iterator iter = isotope_distribution.begin(); iter != isotope_distribution.end(); ++iter)
 
126
    {
 
127
      isotopes_exact.push_back(iter->second);
 
128
    }
 
129
 
 
130
    // "stretch" the averagine isotope distribution
 
131
    Size isotopes_exact_size = isotopes_exact.size();
 
132
    isotopes_exact.resize(Size((isotopes_exact_size - 1)
 
133
                               * isotope_distance_ / interpolation_step_ + 1.6));                             // round up a bit more
 
134
 
 
135
    for (Size i = isotopes_exact_size - 1; i; --i)
 
136
    {
 
137
      // we don't need to move the 0-th entry
 
138
      isotopes_exact[Size(CoordinateType(i) *
 
139
                          isotope_distance_ / interpolation_step_ / charge_ + 0.5)]
 
140
        =   isotopes_exact[i];
 
141
      isotopes_exact[i] = 0;
 
142
    }
 
143
 
 
144
    // compute the normal distribution (to be added for widening the averagine isotope distribution)
 
145
    Math::BasicStatistics<> normal_widening_model;
 
146
    normal_widening_model.setSum(1);
 
147
    normal_widening_model.setMean(0);
 
148
    normal_widening_model.setVariance(isotope_stdev_ * isotope_stdev_);
 
149
    // fill a container with CoordinateType points
 
150
    ContainerType normal_widening_coordinate;
 
151
    for (DoubleReal coord = -normal_widening_width;
 
152
         coord <= normal_widening_width;
 
153
         coord += interpolation_step_
 
154
         )
 
155
    {
 
156
      normal_widening_coordinate.push_back(coord);
 
157
    }
 
158
    // compute normal approximation at these CoordinateType points
 
159
    ContainerType normal_widening;
 
160
    normal_widening_model.normalApproximation(normal_widening, normal_widening_coordinate);
 
161
 
 
162
    // fill linear interpolation
 
163
    const ContainerType & left = isotopes_exact;
 
164
    const ContainerType & right = normal_widening;
 
165
    ContainerType & result = interpolation_.getData();
 
166
    result.clear();
 
167
 
 
168
    Int rMax = std::min(Int(left.size() + right.size() - 1), Int(2 * normal_widening_width / interpolation_step_ * max_isotope_ + 1));
 
169
    result.resize(rMax, 0);
 
170
 
 
171
    // we loop backwards because then the small products tend to come first
 
172
    // (for better numerics)
 
173
    for (SignedSize i = left.size() - 1; i >= 0; --i)
 
174
    {
 
175
      if (left[i] == 0)
 
176
        continue;
 
177
      for (SignedSize j = std::min<SignedSize>(rMax - i, right.size()) - 1; j >= 0; --j)
 
178
      {
 
179
        result[i + j] += left[i] * right[j];
 
180
      }
 
181
    }
 
182
 
 
183
    // set interpolation
 
184
    interpolation_.setMapping(interpolation_step_, normal_widening_width / interpolation_step_, monoisotopic_mz_);
 
185
 
 
186
    // scale data so that integral over distribution equals one
 
187
    // multiply sum by interpolation_step_ -> rectangular approximation of integral
 
188
    IntensityType factor = scaling_ / interpolation_step_ /
 
189
                           std::accumulate(result.begin(), result.end(), IntensityType(0));
 
190
 
 
191
    for (ContainerType::iterator iter = result.begin(); iter != result.end(); ++iter)
 
192
    {
 
193
      *iter *= factor;
 
194
    }
 
195
 
 
196
  }
 
197
 
 
198
  void ExtendedIsotopeModel::setOffset(CoordinateType offset)
 
199
  {
 
200
    DoubleReal diff = offset - getInterpolation().getOffset();
 
201
    monoisotopic_mz_ += diff;
 
202
 
 
203
    InterpolationModel::setOffset(offset);
 
204
 
 
205
    param_.setValue("isotope:monoisotopic_mz", monoisotopic_mz_);
 
206
  }
 
207
 
 
208
  ExtendedIsotopeModel::CoordinateType ExtendedIsotopeModel::getOffset()
 
209
  {
 
210
    return getInterpolation().getOffset();
 
211
  }
 
212
 
 
213
  UInt ExtendedIsotopeModel::getCharge()
 
214
  {
 
215
    return charge_;
 
216
  }
 
217
 
 
218
  ExtendedIsotopeModel::CoordinateType ExtendedIsotopeModel::getCenter() const
 
219
  {
 
220
    return monoisotopic_mz_;
 
221
  }
 
222
 
 
223
  void ExtendedIsotopeModel::updateMembers_()
 
224
  {
 
225
    InterpolationModel::updateMembers_();
 
226
 
 
227
    charge_ = param_.getValue("charge");
 
228
    isotope_stdev_ = param_.getValue("isotope:stdev");
 
229
    monoisotopic_mz_ = param_.getValue("isotope:monoisotopic_mz");
 
230
    max_isotope_ = param_.getValue("isotope:maximum");
 
231
    trim_right_cutoff_ = param_.getValue("isotope:trim_right_cutoff");
 
232
    isotope_distance_ = param_.getValue("isotope:distance");
 
233
 
 
234
    averagine_[C] = param_.getValue("averagines:C");
 
235
    averagine_[H] = param_.getValue("averagines:H");
 
236
    averagine_[N] = param_.getValue("averagines:N");
 
237
    averagine_[O] = param_.getValue("averagines:O");
 
238
    averagine_[S] = param_.getValue("averagines:S");
 
239
 
 
240
    setSamples();
 
241
  }
 
242
 
228
243
}