← Back to branch summary

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

« back to all changes in this revision

Viewing changes to include/OpenMS/FILTERING/DATAREDUCTION/SILACAnalyzer.h

  • 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:
include/OpenMS/FILTERING/DATAREDUCTION/SILACAnalyzer.h
 
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: Lars Nilse $
 
32
// $Authors: Lars Nilse, Steffen Sass, Holger Plattfaut, Bastian Blank $
 
33
// --------------------------------------------------------------------------
 
34
 
 
35
#ifndef OPENMS_FILTERING_DATAREDUCTION_SILACANALYZER_H
 
36
#define OPENMS_FILTERING_DATAREDUCTION_SILACANALYZER_H
 
37
 
 
38
//OpenMS includes
 
39
#include <OpenMS/config.h>
 
40
// #include <OpenMS/APPLICATIONS/TOPPBase.h>
 
41
#include <OpenMS/CONCEPT/ProgressLogger.h>
 
42
#include <OpenMS/DATASTRUCTURES/DBoundingBox.h>
 
43
#include <OpenMS/KERNEL/MSExperiment.h>
 
44
#include <OpenMS/KERNEL/StandardTypes.h>
 
45
#include <OpenMS/KERNEL/ConsensusMap.h>
 
46
#include <OpenMS/KERNEL/FeatureMap.h>
 
47
#include <OpenMS/FORMAT/MzMLFile.h>
 
48
#include <OpenMS/FORMAT/ConsensusXMLFile.h>
 
49
#include <OpenMS/FORMAT/FeatureXMLFile.h>
 
50
#include <OpenMS/MATH/STATISTICS/LinearRegression.h>
 
51
#include <OpenMS/KERNEL/RangeUtils.h>
 
52
#include <OpenMS/KERNEL/ChromatogramTools.h>
 
53
#include <OpenMS/FORMAT/MzQuantMLFile.h>
 
54
#include <OpenMS/METADATA/MSQuantifications.h>
 
55
 
 
56
#include <OpenMS/FILTERING/DATAREDUCTION/SILACFilter.h>
 
57
#include <OpenMS/FILTERING/DATAREDUCTION/SILACFiltering.h>
 
58
#include <OpenMS/COMPARISON/CLUSTERING/SILACClustering.h>
 
59
#include <OpenMS/TRANSFORMATIONS/FEATUREFINDER/PeakWidthEstimator.h>
 
60
 
 
61
//Contrib includes
 
62
#include <boost/algorithm/string/split.hpp>
 
63
#include <boost/algorithm/string/classification.hpp>
 
64
 
 
65
//std includes
 
66
#include <cmath>
 
67
#include <vector>
 
68
#include <algorithm>
 
69
#include <fstream>
 
70
#include <limits>
 
71
#include <locale>
 
72
#include <iomanip>
 
73
 
 
74
namespace OpenMS
 
75
{
 
76
  /**
 
77
   * @brief Algorithm to use for SILACAnalysis
 
78
   *
 
79
   * Please initialize before usage using the initialize function.
 
80
   * Next, one can estimate the peak width before filtering the data. 
 
81
   *
 
82
   * In the final step, clusterData will generate the output data.
 
83
   *
 
84
   * @see SILACFiltering
 
85
   */
 
86
  class OPENMS_DLLAPI SILACAnalyzer :
 
87
    public ProgressLogger
 
88
  {
 
89
  private:
 
90
 
 
91
    // input and output files
 
92
    String in;
 
93
    String out;
 
94
    String out_clusters;
 
95
    String out_features;
 
96
    String out_mzq;
 
97
 
 
98
    String out_filters;
 
99
    String in_filters;
 
100
    String out_debug;
 
101
 
 
102
    // section "sample"
 
103
    String selected_labels;
 
104
    UInt charge_min;
 
105
    UInt charge_max;
 
106
    Int missed_cleavages;
 
107
    UInt isotopes_per_peptide_min;
 
108
    UInt isotopes_per_peptide_max;
 
109
 
 
110
    // section "algorithm"
 
111
    DoubleReal rt_threshold;
 
112
    DoubleReal rt_min;
 
113
    DoubleReal intensity_cutoff;
 
114
    DoubleReal intensity_correlation;
 
115
    DoubleReal model_deviation;
 
116
    bool allow_missing_peaks;
 
117
 
 
118
    // section "labels"
 
119
    std::vector<std::vector<String> > SILAClabels;    // list of SILAC labels, e.g. selected_labels="[Lys4,Arg6][Lys8,Arg10]" => SILAClabels[0][1]="Arg6"
 
120
    std::vector<std::vector<DoubleReal> > massShifts; // list of mass shifts
 
121
 
 
122
    typedef SILACClustering Clustering;
 
123
 
 
124
    MSQuantifications msq;
 
125
 
 
126
  public:
 
127
    SILACAnalyzer() :
 
128
      allow_missing_peaks(true)
 
129
    {
 
130
    }
 
131
 
 
132
    /**
 
133
     * @brief Initializes the algorithm with parameters
 
134
     *
 
135
     * @param selected_labels_ Labels used for labelling the sample. For example, [Lys4,Arg6][Lys8,Arg10] describes a mixtures of three samples. One of them unlabelled, one labelled with Lys4 and Arg6 and a third one with Lys8 and Arg10. The used labels must be described in the label_identifiers parameter
 
136
     * @param label_identifiers a map of labels and corresponding mass shits e.g. "Arg6" => 6.0201290268
 
137
     
 
138
     */
 
139
    void initialize(
 
140
      // section "sample"
 
141
      String selected_labels_,
 
142
      UInt charge_min_,
 
143
      UInt charge_max_,
 
144
      Int missed_cleavages_,
 
145
      UInt isotopes_per_peptide_min_,
 
146
      UInt isotopes_per_peptide_max_, 
 
147
      // section "algorithm"
 
148
      DoubleReal rt_threshold_,
 
149
      DoubleReal rt_min_,
 
150
      DoubleReal intensity_cutoff_,
 
151
      DoubleReal intensity_correlation_,
 
152
      DoubleReal model_deviation_,
 
153
      bool allow_missing_peaks_,
 
154
      // labels part
 
155
      std::map<String, DoubleReal> label_identifiers)
 
156
    {
 
157
      selected_labels          = selected_labels_;
 
158
      charge_min               = charge_min_;
 
159
      charge_max               = charge_max_;
 
160
      missed_cleavages         = missed_cleavages_;
 
161
      isotopes_per_peptide_min = isotopes_per_peptide_min_;
 
162
      isotopes_per_peptide_max = isotopes_per_peptide_max_;
 
163
 
 
164
      rt_threshold           = rt_threshold_;
 
165
      rt_min                 = rt_min_;
 
166
      intensity_cutoff       = intensity_cutoff_;
 
167
      intensity_correlation  = intensity_correlation_;
 
168
      model_deviation        = model_deviation_;
 
169
      allow_missing_peaks    = allow_missing_peaks_;
 
170
 
 
171
      calculateLabelsAndMassShifts(label_identifiers);
 
172
    }
 
173
 
 
174
    /**
 
175
     * @brief Calculate the internal massShift and label datastructures from a
 
176
     * map of identifiers and mass shifts (e.g. "Arg6" => 6.0201290268).
 
177
     *
 
178
     * Note that this part is necessary for the algorithm to work and this
 
179
     * function is called from the initialize section.  The algorithm has to be
 
180
     * initialized _first_ with the list of actually used labels
 
181
     * (selected_labels) using the initialize_sample call.
 
182
     */
 
183
    void calculateLabelsAndMassShifts(std::map<String, DoubleReal> label_identifiers);
 
184
 
 
185
    void run_all(MSExperiment<Peak1D> & exp, ConsensusMap & out_map)
 
186
    {
 
187
      PeakWidthEstimator::Result peak_width;
 
188
      std::vector<std::vector<SILACPattern> > data;
 
189
      MSQuantifications msq;
 
190
      std::vector<Clustering *> cluster_data;
 
191
 
 
192
      peak_width = estimatePeakWidth(exp);
 
193
      filterData(exp, peak_width, data); 
 
194
      clusterData(exp, peak_width, cluster_data, data);
 
195
 
 
196
      // write output to consensus map
 
197
      for (std::vector<Clustering *>::const_iterator it = cluster_data.begin(); it != cluster_data.end(); ++it)
 
198
      {
 
199
        generateClusterConsensusByCluster(out_map, **it);
 
200
      }
 
201
    }
 
202
 
 
203
    /**
 
204
     * @brief Peak width estimation
 
205
     */
 
206
    PeakWidthEstimator::Result estimatePeakWidth(const MSExperiment<Peak1D> & exp);
 
207
 
 
208
    /**
 
209
     * @brief Filtering
 
210
     */
 
211
    void filterData(MSExperiment<Peak1D> & exp, const PeakWidthEstimator::Result & peak_width, std::vector<std::vector<SILACPattern> > & data);
 
212
 
 
213
    /**
 
214
     * @brief Clustering
 
215
     */
 
216
    void clusterData(const MSExperiment<> &, const PeakWidthEstimator::Result &, std::vector<Clustering *> &, std::vector<std::vector<SILACPattern> > & data);
 
217
 
 
218
 
 
219
    /**
 
220
     * @brief Returns the list of SILAC labels, e.g.
 
221
     * selected_labels="[Lys4,Arg6][Lys8,Arg10]" => SILAClabels[0][1]="Arg6"
 
222
     */
 
223
    std::vector<std::vector<String> > & getSILAClabels()
 
224
    {
 
225
      return SILAClabels;
 
226
    }
 
227
 
 
228
    /**
 
229
     * @brief Returns the list of mass shifts calcualted
 
230
     */
 
231
    std::vector<std::vector<DoubleReal> > & getMassShifts()
 
232
    {
 
233
      return massShifts;
 
234
    }
 
235
 
 
236
  public:
 
237
 
 
238
    /**
 
239
     * @brief Generate ConsensusMap from clustering result
 
240
     */
 
241
    void generateClusterConsensusByCluster(ConsensusMap &, const Clustering &) const;
 
242
 
 
243
    /**
 
244
     * @brief Generate ConsensusMap from clustering result, one consensus per pattern
 
245
     */
 
246
    void generateClusterConsensusByPattern(ConsensusMap &, const Clustering &, UInt & cluster_id) const;
 
247
 
 
248
    /**
 
249
     * @brief Generate debug output from clustering result
 
250
     */
 
251
    void generateClusterDebug(std::ostream & out, const Clustering & clustering, UInt & cluster_id) const;
 
252
 
 
253
  public:
 
254
    /**
 
255
     * @brief Generate ConsensusMap from filter result
 
256
     */
 
257
    void generateFilterConsensusByPattern(ConsensusMap &, const std::vector<SILACPattern> &) const;
 
258
 
 
259
  private:
 
260
 
 
261
    /**
 
262
     * @brief Generate a consensus entry from a pattern
 
263
     */
 
264
    ConsensusFeature generateSingleConsensusByPattern(const SILACPattern &) const;
 
265
 
 
266
 
 
267
  public:
 
268
 
 
269
    /**
 
270
     * @brief Generate FeatureMap from clustering result
 
271
     */
 
272
    void generateClusterFeatureByCluster(FeatureMap<> &, const Clustering &) const;
 
273
 
 
274
    /**
 
275
     * @brief Read filter result from ConsensusMap
 
276
     */
 
277
    void readFilterConsensusByPattern(ConsensusMap &, std::vector<std::vector<SILACPattern> > &);
 
278
 
 
279
    static const String & selectColor(UInt nr);
 
280
 
 
281
    /**
 
282
     * @brief Read consensusXML from file to ConsensusMap
 
283
     */
 
284
    void readConsensus(const String & filename, ConsensusMap & in) const
 
285
    {
 
286
      ConsensusXMLFile c_file;
 
287
      c_file.load(filename, in);
 
288
    }
 
289
 
 
290
    /**
 
291
     * @brief Write consensusXML from ConsensusMap to file
 
292
     */
 
293
    void writeConsensus(const String & filename, ConsensusMap & out) const
 
294
    {
 
295
      out.sortByPosition();
 
296
      out.applyMemberFunction(&UniqueIdInterface::setUniqueId);
 
297
      out.setExperimentType("silac");
 
298
 
 
299
      ConsensusXMLFile c_file;
 
300
      c_file.store(filename, out);
 
301
    }
 
302
 
 
303
    /**
 
304
     * @brief Write MzQuantML from ConsensusMap to file
 
305
     */
 
306
    void writeMzQuantML(const String & filename, MSQuantifications & msq) const
 
307
    {
 
308
      //~ TODO apply above to ConsensusMap befor putting into Msq
 
309
      //~ out.sortByPosition();
 
310
      //~ out.applyMemberFunction(&UniqueIdInterface::setUniqueId);
 
311
      //~ out.setExperimentType("SILAC");
 
312
 
 
313
      MzQuantMLFile file;
 
314
      file.store(filename, msq);
 
315
    }
 
316
 
 
317
    /**
 
318
     * @brief Write featureXML from FeatureMap to file
 
319
     */
 
320
    void writeFeatures(const String & filename, FeatureMap<> & out) const
 
321
    {
 
322
      out.sortByPosition();
 
323
      out.applyMemberFunction(&UniqueIdInterface::setUniqueId);
 
324
 
 
325
      FeatureXMLFile f_file;
 
326
      f_file.store(filename, out);
 
327
    }
 
328
 
 
329
  };
 
330
}
 
331
 
 
332
#endif /* SILACANALYZER_H_ */