C++ PeakMap::RTBegin方法代码示例

  void FeatureFinder::run(const String& algorithm_name, PeakMap& input_map, FeatureMap& features, const Param& param, const FeatureMap& seeds)
  {
    // Nothing to do if there is no data
    if ((algorithm_name != "mrm" && input_map.empty()) || (algorithm_name == "mrm" && input_map.getChromatograms().empty()))
    {
      features.clear(true);
      return;
    }

    // check input
    {
      // We need updated ranges => check number of peaks
      if (algorithm_name != "mrm" && input_map.getSize() == 0)
      {
        throw Exception::IllegalArgument(__FILE__, __LINE__, OPENMS_PRETTY_FUNCTION, "FeatureFinder needs updated ranges on input map. Aborting.");
      }

      // We need MS1 data only => check levels
      if (algorithm_name != "mrm" && (input_map.getMSLevels().size() != 1 || input_map.getMSLevels()[0] != 1))
      {
        throw Exception::IllegalArgument(__FILE__, __LINE__, OPENMS_PRETTY_FUNCTION, "FeatureFinder can only operate on MS level 1 data. Please do not use MS/MS data. Aborting.");
      }

      //Check if the peaks are sorted according to m/z
      if (!input_map.isSorted(true))
      {
        LOG_WARN << "Input map is not sorted by RT and m/z! This is done now, before applying the algorithm!" << std::endl;
        input_map.sortSpectra(true);
        input_map.sortChromatograms(true);
      }
      for (Size s = 0; s < input_map.size(); ++s)
      {
        if (input_map[s].empty())
          continue;
        if (input_map[s][0].getMZ() < 0)
        {
          throw Exception::IllegalArgument(__FILE__, __LINE__, OPENMS_PRETTY_FUNCTION, "FeatureFinder can only operate on spectra that contain peaks with positive m/z values. Filter the data accordingly beforehand! Aborting.");
        }
      }
    }

    // initialize
    if (algorithm_name != "mrm" && algorithm_name != "centroided")
    {
      // Resize peak flag vector
      flags_.resize(input_map.size());
      for (Size i = 0; i < input_map.size(); ++i)
      {
        flags_[i].assign(input_map[i].size(), UNUSED);
      }
    }

    // do the work
    if (algorithm_name != "none")
    {
      FeatureFinderAlgorithm* algorithm = Factory<FeatureFinderAlgorithm>::create(algorithm_name);
      algorithm->setParameters(param);
      algorithm->setData(input_map, features, *this);
      algorithm->setSeeds(seeds);
      algorithm->run();
      delete(algorithm);
    }

    if (algorithm_name != "mrm") // mrm  works on chromatograms; the next section is only for conventional data
    {
      //report RT apex spectrum index and native ID for each feature
      for (Size i = 0; i < features.size(); ++i)
      {
        //index
        Size spectrum_index = input_map.RTBegin(features[i].getRT()) - input_map.begin();
        features[i].setMetaValue("spectrum_index", spectrum_index);
        //native id
        if (spectrum_index < input_map.size())
        {
          String native_id = input_map[spectrum_index].getNativeID();
          features[i].setMetaValue("spectrum_native_id", native_id);
        }
        else
        {
          /// @todo that happens sometimes using IsotopeWaveletFeatureFinder (Rene, Marc, Andreas, Clemens)
          std::cerr << "FeatureFinderAlgorithm_impl, line=" << __LINE__ << "; FixMe this cannot be, but happens" << std::endl;
        }
      }
    }
  }

    set<Size> correctToNearestMS1Peak(PeakMap & exp, double mz_tolerance, bool ppm, vector<double> & deltaMZs, vector<double> & mzs, vector<double> & rts)
    {
      set<Size> corrected_precursors;
      // load experiment and extract precursors
      vector<Precursor> precursors;  // precursor
      vector<double> precursors_rt;  // RT of precursor MS2 spectrum
      vector<Size> precursor_scan_index;
      getPrecursors_(exp, precursors, precursors_rt, precursor_scan_index);

      for (Size i = 0; i != precursors_rt.size(); ++i)
      {
        // get precursor rt
        double rt = precursors_rt[i];

        // get precursor MZ
        double mz = precursors[i].getMZ();

        //cout << rt << " " << mz << endl;

        // get precursor spectrum
        MSExperiment<Peak1D>::ConstIterator rt_it = exp.RTBegin(rt - 1e-8);

        // store index of MS2 spectrum
        UInt precursor_spectrum_idx = rt_it - exp.begin();

        // get parent (MS1) of precursor spectrum
        rt_it = exp.getPrecursorSpectrum(rt_it);

        if (rt_it->getMSLevel() != 1)
        {
          LOG_WARN << "Error: no MS1 spectrum for this precursor" << endl;
        }

        //cout << rt_it->getRT() << " " << rt_it->size() << endl;

        // find peak (index) closest to expected position
        Size nearest_peak_idx = rt_it->findNearest(mz);

        // get actual position of closest peak
        double nearest_peak_mz = (*rt_it)[nearest_peak_idx].getMZ();

        // calculate error between expected and actual position
        double nearestPeakError = ppm ? abs(nearest_peak_mz - mz)/mz * 1e6 : abs(nearest_peak_mz - mz);

        // check if error is small enough
        if (nearestPeakError < mz_tolerance)
        {
          // sanity check: do we really have the same precursor in the original and the picked spectrum
          if (fabs(exp[precursor_spectrum_idx].getPrecursors()[0].getMZ() - mz) > 0.0001)
          {
            LOG_WARN << "Error: index is referencing different precursors in original and picked spectrum." << endl;
          }

          // cout << mz << " -> " << nearest_peak_mz << endl;
          double deltaMZ = nearest_peak_mz - mz;
          deltaMZs.push_back(deltaMZ);
          mzs.push_back(mz);
          rts.push_back(rt);
          // correct entries
          Precursor corrected_prec = precursors[i];
          corrected_prec.setMZ(nearest_peak_mz);
          exp[precursor_spectrum_idx].getPrecursors()[0] = corrected_prec;
          corrected_precursors.insert(precursor_spectrum_idx);
        }
      }
      return corrected_precursors;
    }

  void correct(PeakMap & exp, vector<DoubleReal> & deltaMZs, vector<DoubleReal> & mzs, vector<DoubleReal> & rts)
  {
    // load experiment and extract precursors
    vector<Precursor> precursors;  // precursor
    vector<double> precursors_rt;  // RT of precursor MS2 spectrum
    getPrecursors_(exp, precursors, precursors_rt);

    for (Size i = 0; i != precursors_rt.size(); ++i)
    {
      // get precursor rt
      DoubleReal rt = precursors_rt[i];

      // get precursor MZ
      DoubleReal mz = precursors[i].getMZ();

      //cout << rt << " " << mz << endl;

      // get precursor spectrum
      MSExperiment<Peak1D>::ConstIterator rt_it = exp.RTBegin(rt);

      // store index of MS2 spectrum
      UInt precursor_spectrum_idx = rt_it - exp.begin();

      // get parent (MS1) of precursor spectrum
      rt_it = exp.getPrecursorSpectrum(rt_it);

      if (rt_it->getMSLevel() != 1)
      {
        cout << "Error: no MS1 spectrum for this precursor" << endl;
      }

      //cout << rt_it->getRT() << " " << rt_it->size() << endl;

      // find peak (index) closest to expected position
      Size nearest_peak_idx = rt_it->findNearest(mz);

      // get actual position of closest peak
      DoubleReal nearest_peak_mz = (*rt_it)[nearest_peak_idx].getMZ();

      // calculate error between expected and actual position
      DoubleReal nearestPeakError = abs(nearest_peak_mz - mz);

      // check if error is small enough
      if (nearestPeakError < 0.1)
      {
        // sanity check: do we really have the same precursor in the original and the picked spectrum
        if (fabs(exp[precursor_spectrum_idx].getPrecursors()[0].getMZ() - mz) > 0.0001)
        {
          cout << "Error: index is referencing different precursors in original and picked spectrum." << endl;
        }

        // cout << mz << " -> " << nearest_peak_mz << endl;
        DoubleReal deltaMZ = nearest_peak_mz - mz;
        deltaMZs.push_back(deltaMZ);
        mzs.push_back(mz);
        rts.push_back(rt);
        // correct entries
        Precursor corrected_prec = precursors[i];
        corrected_prec.setMZ(nearest_peak_mz);
        exp[precursor_spectrum_idx].getPrecursors()[0] = corrected_prec;
      }
    }
  }