OpenMS: PSLPFormulation.h Source File

Go to the documentation of this file.
00001 // --------------------------------------------------------------------------
00002 //                   OpenMS -- Open-Source Mass Spectrometry
00003 // --------------------------------------------------------------------------
00004 // Copyright The OpenMS Team -- Eberhard Karls University Tuebingen,
00005 // ETH Zurich, and Freie Universitaet Berlin 2002-2012.
00006 //
00007 // This software is released under a three-clause BSD license:
00008 //  * Redistributions of source code must retain the above copyright
00009 //    notice, this list of conditions and the following disclaimer.
00010 //  * Redistributions in binary form must reproduce the above copyright
00011 //    notice, this list of conditions and the following disclaimer in the
00012 //    documentation and/or other materials provided with the distribution.
00013 //  * Neither the name of any author or any participating institution
00014 //    may be used to endorse or promote products derived from this software
00015 //    without specific prior written permission.
00016 // For a full list of authors, refer to the file AUTHORS.
00017 // --------------------------------------------------------------------------
00018 // THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
00019 // AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
00020 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
00021 // ARE DISCLAIMED. IN NO EVENT SHALL ANY OF THE AUTHORS OR THE CONTRIBUTING
00022 // INSTITUTIONS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
00023 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
00024 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS;
00025 // OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY,
00026 // WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR
00027 // OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF
00028 // ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
00029 //
00030 // --------------------------------------------------------------------------
00031 // $Maintainer: Alexandra Zerck $
00032 // $Authors: $
00033 // --------------------------------------------------------------------------
00034 
00035 #ifndef OPENMS_ANALYSIS_TARGETED_PSLPFORMULATION_H
00036 #define OPENMS_ANALYSIS_TARGETED_PSLPFORMULATION_H
00037 #include <OpenMS/KERNEL/FeatureMap.h>
00038 #include <OpenMS/KERNEL/MSExperiment.h>
00039 
00040 #include <OpenMS/DATASTRUCTURES/LPWrapper.h>
00041 
00042 namespace OpenMS
00043 {
00044   class PrecursorIonSelectionPreprocessing;
00045   class PSProteinInference;
00051   class OPENMS_DLLAPI PSLPFormulation :
00052     public DefaultParamHandler
00053   {
00054 
00055 
00056 
00057 public:
00058 
00059     PSLPFormulation();
00060 
00061     virtual ~PSLPFormulation();
00062 
00066     struct IndexTriple
00067     {
00068       Size feature;
00069       Int scan;
00070       Size variable;
00071       DoubleReal rt_probability;
00072       DoubleReal signal_weight;
00073       String prot_acc;
00074     };
00075 
00076 
00088     template <typename InputPeakType>
00089     void createAndSolveILPForKnownLCMSMapFeatureBased(const FeatureMap<> & features,
00090                                                       const MSExperiment<InputPeakType> & experiment,
00091                                                       std::vector<IndexTriple> & variable_indices,
00092                                                       std::vector<std::vector<std::pair<Size, Size> > > & mass_ranges,
00093                                                       std::set<Int> & charges_set, UInt ms2_spectra_per_rt_bin,
00094                                                       std::vector<int> & solution_indices);
00095 
00100     void createAndSolveILPForInclusionListCreation(PrecursorIonSelectionPreprocessing & preprocessing,
00101                                                    UInt ms2_spectra_per_rt_bin, UInt max_list_size,
00102                                                    FeatureMap<> & precursors,
00103                                                    bool solve_ILP = true);
00104 
00105     template <typename InputPeakType>
00106     void createAndSolveCombinedLPForKnownLCMSMapFeatureBased(const FeatureMap<> & features,
00107                                                              const MSExperiment<InputPeakType> & experiment,
00108                                                              std::vector<IndexTriple> & variable_indices,
00109                                                              std::vector<int> & solution_indices,
00110                                                              std::vector<std::vector<std::pair<Size, Size> > > & mass_ranges,
00111                                                              std::set<Int> & charges_set, UInt ms2_spectra_per_rt_bin,
00112                                                              Size step_size = 0, bool sequential_order = false);
00113 
00114     void updateStepSizeConstraint(Size iteration, UInt step_size);
00115     void updateFeatureILPVariables(FeatureMap<> & new_features, std::vector<IndexTriple> & variable_indices, std::map<Size, std::vector<String> > & feature_constraints_map);
00116     void updateRTConstraintsForSequentialILP(Size & rt_index, UInt ms2_spectra_per_rt_bin, Size max_rt_index);
00117     void updateCombinedILP(FeatureMap<> & features, PrecursorIonSelectionPreprocessing & preprocessed_db, std::vector<IndexTriple> & variable_indices,
00118                            std::vector<String> & new_protein_accs, std::vector<String> & protein_accs, PSProteinInference & prot_inference, Size & variable_counter,
00119                            std::map<String, std::vector<Size> > & protein_feature_map, Feature & new_feature, std::map<String, Size> & protein_variable_index_map,
00120                            std::map<String, std::set<String> > & prot_id_counter);
00121 
00122 
00126     void solveILP(std::vector<int> & solution_indices);
00127 
00128     void setLPSolver(LPWrapper::SOLVER solver)
00129     {
00130       solver_ = solver;
00131     }
00132 
00133     LPWrapper::SOLVER getLPSolver()
00134     {
00135       return solver_;
00136     }
00137 
00138     struct IndexLess :
00139       std::binary_function<IndexTriple, IndexTriple, bool>
00140     {
00141       inline bool operator()(IndexTriple  const & left,
00142                              IndexTriple const & right) const
00143       {
00144         return left.feature < right.feature;
00145       }
00146 
00147     };
00148 
00149 
00150     struct ScanLess :
00151       std::binary_function<IndexTriple, IndexTriple, bool>
00152     {
00153       inline bool operator()(IndexTriple  const & left,
00154                              IndexTriple  const & right) const
00155       {
00156         return left.scan < right.scan;
00157       }
00158 
00159     };
00160 
00161     struct VariableIndexLess :
00162       std::binary_function<IndexTriple, IndexTriple, bool>
00163     {
00164       inline bool operator()(IndexTriple  const & left,
00165                              IndexTriple  const & right) const
00166       {
00167         return left.variable < right.variable;
00168       }
00169 
00170     };
00171 
00172 protected:
00173 
00174     template <typename InputPeakType>
00175     void getXIC_(const std::vector<std::pair<Size, Size> > & end_points,
00176                  std::vector<DoubleReal> & weights,
00177                  const MSExperiment<InputPeakType> & experiment,
00178                  const bool normalize);
00179 
00183     template <typename InputPeakType>
00184     void calculateXICs_(std::vector<std::vector<DoubleReal> > & xics,
00185                         const FeatureMap<> & features,
00186                         const MSExperiment<InputPeakType> & experiment,
00187                         const std::vector<std::vector<std::pair<Size, Size> > > & mass_ranges,
00188                         const bool normalize);
00189 
00193     void createAndSolveILP_(const FeatureMap<> & features, std::vector<std::vector<DoubleReal> > & intensity_weights,
00194                             std::set<Int> & charges_set, std::vector<std::vector<std::pair<Size, Size> > > & mass_ranges,
00195                             std::vector<IndexTriple> & variable_indices, std::vector<int> & solution_indices,
00196                             UInt ms2_spectra_per_rt_bin, Size number_of_scans);
00197 
00198     void createAndSolveCombinedLPFeatureBased_(const FeatureMap<> & features, std::vector<std::vector<DoubleReal> > & intensity_weights,
00199                                                std::set<Int> & charges_set, std::vector<std::vector<std::pair<Size, Size> > > & mass_ranges,
00200                                                std::vector<IndexTriple> & variable_indices, std::vector<Int> & solution_indices,
00201                                                UInt ms2_spectra_per_rt_bin, Size number_of_scans, Size step_size = 0, bool sequential_order = false);
00202 
00203     void addProteinToILP_(PrecursorIonSelectionPreprocessing & preprocessing,
00204                           std::map<String, std::vector<DoubleReal> >::const_iterator map_iter,
00205                           Size & counter, Size & pep_counter, Size & feature_counter,
00206                           std::vector<IndexTriple> & variable_indices,
00207                           std::map<String, Size> & protein_penalty_index_map, FeatureMap<> & precursors);
00208 
00209     void addPrecursorAcquisitionNumberConstraint_(std::vector<IndexTriple> & variable_indices, Size number_of_features, UInt number_of_msms_per_precursor);
00210 
00211     void addMaxInclusionListSizeConstraints_(std::vector<IndexTriple> & variable_indices, /*Size number_of_features,*/ UInt max_list_size);
00212 
00213     void addRTBinCapacityConstraint_(std::vector<IndexTriple> & variable_indices,
00214                                      Size max_rt_index, UInt ms2_spectra_per_rt_bin, bool sequential_order = false);
00215 
00216     void addProteinCoverageConstraint_(std::vector<IndexTriple> & variable_indices,
00217                                        PrecursorIonSelectionPreprocessing & preprocessing,
00218                                        std::map<String, Size> protein_variable_index_map);
00219 
00220     void addStepSizeConstraint_(std::vector<IndexTriple> & variable_indices, UInt step_size);
00221 
00222 
00223     void assembleInclusionListForProteinBasedLP_(std::vector<IndexTriple> & variable_indices, FeatureMap<> & precursors, std::vector<int> & solution_indices, PrecursorIonSelectionPreprocessing & preprocessing);
00224 
00225     void updateObjFunction_(String acc, FeatureMap<> & features, PrecursorIonSelectionPreprocessing & preprocessed_db, std::vector<IndexTriple> & variable_indices);
00226 
00227 
00228     Int getNumberOfPrecsInSpectrum_(Int constr_idx);
00229 
00230     LPWrapper * model_;
00231     LPWrapper::SOLVER solver_;
00232   };
00233 
00234   template <typename InputPeakType>
00235   void PSLPFormulation::getXIC_(const std::vector<std::pair<Size, Size> > & end_points,
00236                                 std::vector<DoubleReal> & weights,
00237                                 const MSExperiment<InputPeakType> & experiment,
00238                                 const bool normalize)
00239   {
00240     DoubleReal max_weight = 0.;
00241     weights.clear();
00242     for (Size i = 0; i < end_points.size(); i += 2)
00243     {
00244       DoubleReal weight = 0.;
00245       for (Size j = end_points[i].second; j <= end_points[i + 1].second; ++j)
00246       {
00247         weight += experiment[end_points[i].first][j].getIntensity();
00248         // std::cout << " add "<<experiment[end_points[i].first][j].getIntensity()<<std::endl;
00249       }
00250       if (weight > max_weight)
00251         max_weight = weight;
00252 
00253       weights.push_back(weight);
00254     }
00255 
00256     if (normalize)
00257     {
00258       // normalize weights
00259       for (Size i = 0; i < weights.size(); ++i)
00260       {
00261 #ifdef DEBUG_OPS
00262         if (end_points.size() >= i)
00263         {
00264           std::cout << "scan " << end_points[i].first << " " << weights[i] << " " << max_weight
00265                     << " " << weights[i] / max_weight << std::endl;
00266         }
00267 #endif
00268         weights[i] /= max_weight;
00269       }
00270     }
00271   }
00272 
00273   template <typename InputPeakType>
00274   void PSLPFormulation::calculateXICs_(std::vector<std::vector<DoubleReal> > & xics,
00275                                        const FeatureMap<> & features,
00276                                        const MSExperiment<InputPeakType> & experiment,
00277                                        const std::vector<std::vector<std::pair<Size, Size> > > & mass_ranges,
00278                                        const bool normalize)
00279   {
00280     xics.clear();
00281     xics.resize(features.size());
00282     for (Size i = 0; i < features.size(); ++i)
00283     {
00284       getXIC_(mass_ranges[i], xics[i], experiment, normalize);
00285     }
00286   }
00287 
00288   template <typename InputPeakType>
00289   void PSLPFormulation::createAndSolveILPForKnownLCMSMapFeatureBased(const FeatureMap<> & features,
00290                                                                      const MSExperiment<InputPeakType> & experiment,
00291                                                                      std::vector<IndexTriple> & variable_indices,
00292                                                                      std::vector<std::vector<std::pair<Size, Size> > > & mass_ranges,
00293                                                                      std::set<Int> & charges_set, UInt ms2_spectra_per_rt_bin,
00294                                                                      std::vector<int> & solution_indices)
00295   {
00296 
00297     std::vector<std::vector<DoubleReal> > intensity_weights;
00298     calculateXICs_(intensity_weights, features, experiment, mass_ranges, true);
00299 #ifdef DEBUG_OPS
00300     std::cout << "got xics" << std::endl;
00301 #endif
00302 
00303     createAndSolveILP_(features, intensity_weights, charges_set, mass_ranges, variable_indices, solution_indices,
00304                        ms2_spectra_per_rt_bin, experiment.size());
00305   }
00306 
00307   inline OPENMS_DLLAPI std::ostream & operator<<(std::ostream & os, const PSLPFormulation::IndexTriple & triple)
00308   {
00309     os << "feature: " << triple.feature << " scan: " << triple.scan << " variable: " << triple.variable << " prot_acc: " << triple.prot_acc;
00310     return os;
00311   }
00312 
00313   template <typename InputPeakType>
00314   void PSLPFormulation::createAndSolveCombinedLPForKnownLCMSMapFeatureBased(const FeatureMap<> & features,
00315                                                                             const MSExperiment<InputPeakType> & experiment,
00316                                                                             std::vector<IndexTriple> & variable_indices,
00317                                                                             std::vector<Int> & solution_indices,
00318                                                                             std::vector<std::vector<std::pair<Size, Size> > > & mass_ranges,
00319                                                                             std::set<Int> & charges_set, UInt ms2_spectra_per_rt_bin,
00320                                                                             Size step_size, bool sequential_order)
00321   {
00322 
00323     std::vector<std::vector<DoubleReal> > intensity_weights;
00324     calculateXICs_(intensity_weights, features, experiment, mass_ranges, true);
00325 #ifdef DEBUG_OPS
00326     std::cout << "got xics" << std::endl;
00327 #endif
00328 
00329     createAndSolveCombinedLPFeatureBased_(features, intensity_weights, charges_set, mass_ranges, variable_indices, solution_indices, ms2_spectra_per_rt_bin,
00330                                           experiment.size(), step_size, sequential_order);
00331   }
00332 
00333 } // namespace
00334 
00335 #endif // OPENMS_ANALYSIS_ID_PSLPFORMULATION_H