← Back to branch summary

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

« back to all changes in this revision

Viewing changes to source/ANALYSIS/QUANTITATION/IsobaricNormalizer.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/QUANTITATION/IsobaricNormalizer.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: Stephan Aiche $
 
32
// $Authors: Stephan Aiche, Chris Bielow $
 
33
// --------------------------------------------------------------------------
 
34
 
 
35
#include <OpenMS/ANALYSIS/QUANTITATION/IsobaricNormalizer.h>
 
36
 
 
37
namespace OpenMS
 
38
{
 
39
 
 
40
  IsobaricNormalizer::IsobaricNormalizer(const IsobaricQuantitationMethod* const quant_method) :
 
41
    quant_meth_(quant_method)
 
42
  {
 
43
    reference_channel_name_ = quant_meth_->getChannelInformation()[quant_meth_->getReferenceChannel()].name;
 
44
  }
 
45
 
 
46
  IsobaricNormalizer::IsobaricNormalizer(const IsobaricNormalizer& other) :
 
47
    quant_meth_(other.quant_meth_),
 
48
    reference_channel_name_(other.reference_channel_name_)
 
49
  {
 
50
  }
 
51
 
 
52
  IsobaricNormalizer& IsobaricNormalizer::operator=(const IsobaricNormalizer& rhs)
 
53
  {
 
54
    if (this == &rhs)
 
55
      return *this;
 
56
 
 
57
    quant_meth_ = rhs.quant_meth_;
 
58
    reference_channel_name_ = rhs.reference_channel_name_;
 
59
 
 
60
    return *this;
 
61
  }
 
62
 
 
63
  ConsensusFeature::HandleSetType::iterator IsobaricNormalizer::findReferenceChannel_(ConsensusFeature& cf, const ConsensusMap& consensus_map) const
 
64
  {
 
65
    for (ConsensusFeature::HandleSetType::iterator it_elements = cf.begin();
 
66
         it_elements != cf.end();
 
67
         ++it_elements)
 
68
    {
 
69
      if ((Int) consensus_map.getFileDescriptions()[it_elements->getMapIndex()].getMetaValue("channel_name") == reference_channel_name_)
 
70
      {
 
71
        return it_elements;
 
72
      }
 
73
    }
 
74
 
 
75
    return cf.end();
 
76
  }
 
77
 
 
78
  void IsobaricNormalizer::buildVectorIndex_(const ConsensusMap& consensus_map)
 
79
  {
 
80
    // clear old values
 
81
    ref_map_id_ = 0;
 
82
    map_to_vec_index_.clear();
 
83
 
 
84
    Size index = 0;
 
85
    for (ConsensusMap::FileDescriptions::const_iterator file_it = consensus_map.getFileDescriptions().begin();
 
86
         file_it != consensus_map.getFileDescriptions().end();
 
87
         ++file_it)
 
88
    {
 
89
      if ((Int) file_it->second.getMetaValue("channel_name") == reference_channel_name_)
 
90
      {
 
91
        ref_map_id_ = file_it->first;
 
92
      }
 
93
      map_to_vec_index_[file_it->first] = index;
 
94
      ++index;
 
95
    }
 
96
  }
 
97
 
 
98
  void IsobaricNormalizer::collectRatios_(const ConsensusFeature& cf, const Peak2D::IntensityType& ref_intensity)
 
99
  {
 
100
    for (ConsensusFeature::HandleSetType::const_iterator it_elements = cf.begin();
 
101
         it_elements != cf.end();
 
102
         ++it_elements)
 
103
    {
 
104
      if (ref_intensity == 0) //avoid nan's and inf's
 
105
      {
 
106
        if (it_elements->getIntensity() == 0) // 0/0 will give 'nan'
 
107
        {
 
108
          //so leave it out completely (there is no information to be gained)
 
109
        }
 
110
        else // x/0 is 'inf' but std::sort() has problems with that
 
111
        {
 
112
          peptide_ratios_[map_to_vec_index_[it_elements->getMapIndex()]].push_back(std::numeric_limits<Peak2D::IntensityType>::max());
 
113
        }
 
114
      }
 
115
      else // everything seems fine
 
116
      {
 
117
        peptide_ratios_[map_to_vec_index_[it_elements->getMapIndex()]].push_back(it_elements->getIntensity() / ref_intensity);
 
118
      }
 
119
 
 
120
      // control
 
121
      peptide_intensities_[map_to_vec_index_[it_elements->getMapIndex()]].push_back(it_elements->getIntensity());
 
122
    }
 
123
 
 
124
  }
 
125
 
 
126
  void IsobaricNormalizer::computeNormalizationFactors_(std::vector<Peak2D::IntensityType>& normalization_factors)
 
127
  {
 
128
    // ensure that the ref_(ratios|intensities) are sorted
 
129
    std::sort(peptide_ratios_[ref_map_id_].begin(), peptide_ratios_[ref_map_id_].end());
 
130
    std::sort(peptide_intensities_[ref_map_id_].begin(), peptide_intensities_[ref_map_id_].end());
 
131
 
 
132
    // reporting
 
133
    Peak2D::IntensityType max_deviation_from_control = 0;
 
134
 
 
135
    // find MEDIAN of ratios for each channel (store as 0th element in sorted vector)
 
136
    for (Map<Size, Size>::const_iterator it_map = map_to_vec_index_.begin(); it_map != map_to_vec_index_.end(); ++it_map)
 
137
    {
 
138
      // this is solely for readability reasons, the compiler should optimize this anyway
 
139
      const Size vec_idx = it_map->second;
 
140
 
 
141
      // sort vector (partial_sort might improve performance here)
 
142
      std::sort(peptide_ratios_[vec_idx].begin(), peptide_ratios_[vec_idx].end());
 
143
 
 
144
      // save median as first element
 
145
      normalization_factors[vec_idx] = peptide_ratios_[vec_idx][peptide_ratios_[vec_idx].size() / 2];
 
146
 
 
147
      // sort control (intensities)
 
148
      std::sort(peptide_intensities_[vec_idx].begin(), peptide_intensities_[vec_idx].end());
 
149
 
 
150
      // find MEDIAN of control-method (intensities) for each channel
 
151
      peptide_intensities_[vec_idx][0] = peptide_intensities_[vec_idx][peptide_intensities_[vec_idx].size() / 2] /
 
152
                                         peptide_intensities_[ref_map_id_][peptide_intensities_[ref_map_id_].size() / 2];
 
153
 
 
154
      LOG_INFO << "IsobaricNormalizer:  map-id " << (it_map->first) << " has factor " << (normalization_factors[vec_idx]) << " (control: " << (peptide_intensities_[vec_idx][0]) << ")" << std::endl;
 
155
 
 
156
      Peak2D::IntensityType dev = (peptide_ratios_[vec_idx][0] - peptide_intensities_[vec_idx][0]) / normalization_factors[vec_idx];
 
157
      if (fabs(max_deviation_from_control) < fabs(dev))
 
158
      {
 
159
        max_deviation_from_control = dev;
 
160
      }
 
161
    }
 
162
 
 
163
    LOG_INFO << "IsobaricNormalizer: max ratio deviation of alternative method is " << (max_deviation_from_control * 100) << "%\n";
 
164
  }
 
165
 
 
166
  void IsobaricNormalizer::normalize(ConsensusMap& consensus_map)
 
167
  {
 
168
    // determine reference channel as vector index
 
169
    buildVectorIndex_(consensus_map);
 
170
 
 
171
    // build mapping of map_index to ratio_array_index
 
172
    peptide_ratios_.resize(quant_meth_->getNumberOfChannels());
 
173
    peptide_intensities_.resize(quant_meth_->getNumberOfChannels());
 
174
 
 
175
    //build up ratios for each peptide of non-reference channels
 
176
    ConsensusFeature::HandleSetType::iterator ref_it;
 
177
 
 
178
    for (ConsensusMap::Iterator cm_it = consensus_map.begin(); cm_it != consensus_map.end(); ++cm_it)
 
179
    {
 
180
      // find reference index (this is inefficient to do every time,
 
181
      // but the most robust against anyone who tries to change the internals of ConsensusFeature):
 
182
      ref_it = findReferenceChannel_(*cm_it, consensus_map);
 
183
 
 
184
      // reference channel not found in this ConsensusFeature
 
185
      if (ref_it == cm_it->end())
 
186
      {
 
187
        LOG_WARN << "IsobaricNormalizer::normalize() WARNING: ConsensusFeature "
 
188
                 << (cm_it - consensus_map.begin())
 
189
                 << " does not have a reference channel! Skipping"
 
190
                 << std::endl;
 
191
        continue;
 
192
      }
 
193
 
 
194
      collectRatios_(*cm_it, ref_it->getIntensity());
 
195
    } // ! collect ratios
 
196
 
 
197
    // vector to store the channel wise normalization factors
 
198
    std::vector<Peak2D::IntensityType> normalization_factors;
 
199
    normalization_factors.resize(quant_meth_->getNumberOfChannels());
 
200
 
 
201
    // compute the normalization factors based on the medians of the compute ratios
 
202
    computeNormalizationFactors_(normalization_factors);
 
203
 
 
204
    // free memory
 
205
    peptide_intensities_.clear();
 
206
    peptide_ratios_.clear();
 
207
 
 
208
    // adjust intensity ratios
 
209
    for (size_t i = 0; i < consensus_map.size(); ++i)
 
210
    {
 
211
      // find reference index (this is inefficient to do every time,
 
212
      // but the most robust against anyone who tries to change the
 
213
      // internals of ConsensusFeature):
 
214
      ref_it = findReferenceChannel_(consensus_map[i], consensus_map);
 
215
 
 
216
      // reference channel not found in this ConsensusFeature
 
217
      if (ref_it == consensus_map[i].end())
 
218
      {
 
219
        continue;
 
220
      }
 
221
 
 
222
      // now adjust the ratios
 
223
      ConsensusFeature cf = consensus_map[i];
 
224
      cf.clear(); // delete its handles
 
225
      for (ConsensusFeature::HandleSetType::iterator it_elements = consensus_map[i].begin();
 
226
           it_elements != consensus_map[i].end();
 
227
           ++it_elements)
 
228
      {
 
229
        FeatureHandle hd = *it_elements;
 
230
        if (it_elements == ref_it)
 
231
        {
 
232
          hd.setIntensity(1);
 
233
        }
 
234
        else // divide current intensity by normalization factor (which was stored at position 0)
 
235
        {
 
236
          hd.setIntensity(hd.getIntensity() / normalization_factors[map_to_vec_index_[it_elements->getMapIndex()]]);
 
237
        }
 
238
        cf.insert(hd);
 
239
      }
 
240
      // replace consensusFeature with updated intensity
 
241
      consensus_map[i] = cf;
 
242
    } // ! adjust ratios
 
243
  }
 
244
 
 
245
} // namespace