Skip to content

Precomputation for DIAScoring and DIAHelper #99

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
taranehstrunk wants to merge 1 commit into
base: develop
Choose a base branch
from
Open

Precomputation for DIAScoring and DIAHelper #99

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
17 changes: 10 additions & 7 deletions src/openms/include/OpenMS/ANALYSIS/OPENSWATH/DIAHelper.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -36,6 +36,7 @@

#include <OpenMS/CHEMISTRY/AASequence.h>
#include <OpenMS/OPENSWATHALGO/DATAACCESS/DataStructures.h>
#include <OpenMS/FILTERING/DATAREDUCTION/IsotopeDistributionCache.h>
Copy link
Owner

@cbielow cbielow May 28, 2020

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

no need to include the header here.
Just forward declare should be enough (to save on compile time).
class IsotopeDistributionCache; within the OpenMS namespace.


#include <vector>

Expand Down Expand Up @@ -110,14 +111,16 @@ namespace OpenMS
UInt charge = 1u);

/// get averagine distribution given mass
OPENMS_DLLAPI void getAveragineIsotopeDistribution(const double product_mz,
std::vector<std::pair<double, double> >& isotopesSpec,
const double charge = 1.,
const int nr_isotopes = 4,
const double mannmass = 1.00048);
OPENMS_DLLAPI void getAveragineIsotopeDistribution(IsotopeDistributionCache& iso,
const double product_mz,
std::vector<std::pair<double, double> >& isotopesSpec,
const double charge = 1.,
const int nr_isotopes = 4,
const double mannmass = 1.00048);

/// simulate spectrum from AASequence
OPENMS_DLLAPI void simulateSpectrumFromAASequence(const AASequence& aa,
OPENMS_DLLAPI void simulateSpectrumFromAASequence(IsotopeDistributionCache& iso,
const AASequence& aa,
Copy link
Owner

@cbielow cbielow May 28, 2020

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

it is probably a bit more intuitive to swap argument 1 and 2... because AASequence is the primary input

std::vector<double>& firstIsotopeMasses, //[out]
std::vector<std::pair<double, double> >& isotopeMasses, //[out]
TheoreticalSpectrumGenerator const * g,
Expand All @@ -137,7 +140,7 @@ namespace OpenMS
double charge = 1.);

/// given an experimental spectrum add isotope pattern.
OPENMS_DLLAPI void addIsotopes2Spec(const std::vector<std::pair<double, double> >& spec,
OPENMS_DLLAPI void addIsotopes2Spec(IsotopeDistributionCache& iso, const std::vector<std::pair<double, double> >& spec,
std::vector<std::pair<double, double> >& isotopeMasses, //[out]
double charge = 1.);

Expand Down
7 changes: 5 additions & 2 deletions src/openms/include/OpenMS/ANALYSIS/OPENSWATH/DIAPrescoring.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -44,6 +44,7 @@
#include <OpenMS/OPENSWATHALGO/DATAACCESS/TransitionExperiment.h>

#include <OpenMS/DATASTRUCTURES/DefaultParamHandler.h>
#include <OpenMS/ANALYSIS/OPENSWATH/DIAScoring.h>
Copy link
Owner

@cbielow cbielow May 28, 2020

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

no need for this header.
forward declare IsotopeDistributionCache as before


namespace OpenMS
{
Expand Down Expand Up @@ -84,7 +85,8 @@ namespace OpenMS
and compute manhattan distance and dotprod score between spectrum intensities
and simulated spectrum.
*/
void score(OpenSwath::SpectrumPtr spec,
void score(IsotopeDistributionCache& iso,
OpenSwath::SpectrumPtr spec,
const std::vector<OpenSwath::LightTransition>& lt,
double& dotprod,
double& manhattan);
Expand All @@ -93,7 +95,8 @@ namespace OpenMS
@brief Compute manhattan and dotprod score for all spectra which can be accessed by
the SpectrumAccessPtr for all transitions groups in the LightTargetedExperiment.
*/
void operator()(OpenSwath::SpectrumAccessPtr swath_ptr,
void operator()(IsotopeDistributionCache& iso,
OpenSwath::SpectrumAccessPtr swath_ptr,
OpenSwath::LightTargetedExperiment& transition_exp_used,
OpenSwath::IDataFrameWriter* ivw);
};
Expand Down
8 changes: 7 additions & 1 deletion src/openms/include/OpenMS/ANALYSIS/OPENSWATH/DIAScoring.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -42,6 +42,7 @@
#include <OpenMS/OPENSWATHALGO/DATAACCESS/DataStructures.h>
#include <OpenMS/OPENSWATHALGO/DATAACCESS/ITransition.h>
#include <OpenMS/OPENSWATHALGO/DATAACCESS/TransitionExperiment.h>
#include <OpenMS/FILTERING/DATAREDUCTION/IsotopeDistributionCache.h>

namespace OpenMS
{
Expand Down Expand Up @@ -148,7 +149,10 @@ namespace OpenMS
double& manhattan);
//@}

private:
//function to access private member iso_cache_
IsotopeDistributionCache& getCache();

private:

/// Copy constructor (algorithm class)
DIAScoring(const DIAScoring& rhs);
Expand Down Expand Up @@ -215,6 +219,8 @@ namespace OpenMS
bool dia_extraction_ppm_;
bool dia_centroided_;

IsotopeDistributionCache iso_cache_;

TheoreticalSpectrumGenerator * generator;
};
}
Expand Down
28 changes: 25 additions & 3 deletions src/openms/include/OpenMS/FILTERING/DATAREDUCTION/IsotopeDistributionCache.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -36,6 +36,7 @@

#include <OpenMS/KERNEL/StandardTypes.h>
#include <OpenMS/TRANSFORMATIONS/FEATUREFINDER/FeatureFinderAlgorithmPickedHelperStructs.h>
#include <OpenMS/CHEMISTRY/ISOTOPEDISTRIBUTION/IsotopeDistribution.h>

namespace OpenMS
{
Expand All @@ -47,16 +48,37 @@ namespace OpenMS
public:
typedef FeatureFinderAlgorithmPickedHelperStructs::TheoreticalIsotopePattern TheoreticalIsotopePattern;

IsotopeDistributionCache(double max_mass, double mass_window_width, double intensity_percentage = 0, double intensity_percentage_optional = 0);
/// @name Constructors and Destructors
//@{
/** Default constructor
*/
IsotopeDistributionCache();

/// Destructor
~IsotopeDistributionCache() = default;
//@}


void precalculateDistributionCache(Size num_begin, Size index);
Copy link
Owner

@cbielow cbielow May 28, 2020

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

can you document what the function does and what the parameters mean?


void renormalize( TheoreticalIsotopePattern& isotopes, IsotopeDistribution& isotope_dist);
Copy link
Owner

@cbielow cbielow May 28, 2020

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

docs here as well please


/// Returns the isotope distribution for a certain mass window
const TheoreticalIsotopePattern & getIsotopeDistribution(double mass) const;
const TheoreticalIsotopePattern& getIsotopeDistribution(double mass) ;

const IsotopeDistribution& getIntensity(double mass);
Copy link
Owner

@cbielow cbielow May 28, 2020

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

this is not really an intensity which is returned here, right?
Can you find a better name for the method?
Also document it


private:
/// Vector of pre-calculated isotope distributions for several mass windows
std::vector<TheoreticalIsotopePattern> isotope_distributions_;

std::vector<IsotopeDistribution> distribution_cache_;
Copy link
Owner

@cbielow cbielow May 28, 2020

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

can you document the four members?


double mass_window_width_;

double intensity_percentage_ ;
Copy link
Owner

@cbielow cbielow May 28, 2020

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

what is the percentage? (document)


double intensity_percentage_optional_;
};
}
}//namespace OpenMS

20 changes: 11 additions & 9 deletions src/openms/source/ANALYSIS/OPENSWATH/DIAHelper.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -42,6 +42,9 @@
#include <OpenMS/KERNEL/MSSpectrum.h>
#include <OpenMS/KERNEL/MSExperiment.h>

#include <OpenMS/FILTERING/DATAREDUCTION/IsotopeDistributionCache.h>
#include <include/OpenMS/ANALYSIS/OPENSWATH/DIAScoring.h>

#include <utility>
#include <boost/bind.hpp>

Expand Down Expand Up @@ -252,17 +255,15 @@ namespace OpenMS
}
} // end getBYSeries

void getAveragineIsotopeDistribution(const double product_mz,
void getAveragineIsotopeDistribution(IsotopeDistributionCache& iso,
const double product_mz,
std::vector<std::pair<double, double> >& isotopesSpec,
const double charge,
const int nr_isotopes,
const double mannmass)
{
typedef OpenMS::FeatureFinderAlgorithmPickedHelperStructs::TheoreticalIsotopePattern TheoreticalIsotopePattern;
// create the theoretical distribution
CoarseIsotopePatternGenerator solver(nr_isotopes);
TheoreticalIsotopePattern isotopes;
auto d = solver.estimateFromPeptideWeight(product_mz * charge);
auto d = iso.getIntensity(product_mz * charge);
Copy link
Owner

@cbielow cbielow May 28, 2020

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

can be const auto& d = .... to avoid the copy?


double mass = product_mz;
for (IsotopeDistribution::Iterator it = d.begin(); it != d.end(); ++it)
Expand All @@ -273,29 +274,30 @@ namespace OpenMS
} //end of dia_isotope_corr_sub

//simulate spectrum from AASequence
void simulateSpectrumFromAASequence(const AASequence& aa,
void simulateSpectrumFromAASequence(IsotopeDistributionCache& iso,
const AASequence& aa,
std::vector<double>& firstIsotopeMasses, //[out]
std::vector<std::pair<double, double> >& isotopeMasses, //[out]
TheoreticalSpectrumGenerator const * generator, double charge)
{
getTheorMasses(aa, firstIsotopeMasses, generator, charge);
for (std::size_t i = 0; i < firstIsotopeMasses.size(); ++i)
{
getAveragineIsotopeDistribution(firstIsotopeMasses[i], isotopeMasses,
getAveragineIsotopeDistribution(iso, firstIsotopeMasses[i], isotopeMasses,
charge);
}
}

//given an experimental spectrum add isotope pattern.
void addIsotopes2Spec(const std::vector<std::pair<double, double> >& spec,
void addIsotopes2Spec(IsotopeDistributionCache& iso, const std::vector<std::pair<double, double> >& spec,
std::vector<std::pair<double, double> >& isotopeMasses, //[out]
double charge)
{

for (std::size_t i = 0; i < spec.size(); ++i)
{
std::vector<std::pair<double, double> > isotopes;
Copy link
Owner

@cbielow cbielow May 28, 2020

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

just as a suggestion for speed (even though this was not part of the PR):
move std::vector<std::pair<double, double> > isotopes; in front of the loop and just clear it inside to safe a lot of allocations.

getAveragineIsotopeDistribution(spec[i].first, isotopes, charge);
getAveragineIsotopeDistribution(iso, spec[i].first, isotopes, charge);
for (Size j = 0; j < isotopes.size(); ++j)
{
isotopes[j].second *= spec[i].second; //multiple isotope intensity by spec intensity
Expand Down
10 changes: 6 additions & 4 deletions src/openms/source/ANALYSIS/OPENSWATH/DIAPrescoring.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -72,7 +72,8 @@ namespace OpenMS
}
}

void DiaPrescore::operator()(OpenSwath::SpectrumAccessPtr swath_ptr,
void DiaPrescore::operator()(IsotopeDistributionCache& iso,
OpenSwath::SpectrumAccessPtr swath_ptr,
OpenSwath::LightTargetedExperiment& transition_exp_used,
OpenSwath::IDataFrameWriter* ivw)
{
Expand Down Expand Up @@ -114,7 +115,7 @@ namespace OpenMS
double score1;
double score2;
//OpenSwath::LightPeptide pep;
score(spec, beg->second, score1, score2);
score(iso, spec, beg->second, score1, score2);

score1v.push_back(score1);
score2v.push_back(score2);
Expand All @@ -127,7 +128,8 @@ namespace OpenMS
} //end of forloop over spectra
}

void DiaPrescore::score(OpenSwath::SpectrumPtr spec,
void DiaPrescore::score(IsotopeDistributionCache& iso,
OpenSwath::SpectrumPtr spec,
const std::vector<OpenSwath::LightTransition>& lt,
double& dotprod,
double& manhattan)
Expand All @@ -137,7 +139,7 @@ namespace OpenMS
std::vector<double> firstIstotope, theomasses;
DIAHelpers::extractFirst(res, firstIstotope);
std::vector<std::pair<double, double> > spectrum, spectrum2;
DIAHelpers::addIsotopes2Spec(res, spectrum, nr_charges_);
DIAHelpers::addIsotopes2Spec(iso, res, spectrum, nr_charges_);
spectrum2.resize(spectrum.size());
std::copy(spectrum.begin(), spectrum.end(), spectrum2.begin());
//std::cout << spectrum.size() << std::endl;
Expand Down
36 changes: 8 additions & 28 deletions src/openms/source/ANALYSIS/OPENSWATH/DIAScoring.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -268,7 +268,12 @@ namespace OpenMS
double& dotprod, double& manhattan)
{
OpenMS::DiaPrescore dp(dia_extract_window_, dia_nr_isotopes_, dia_nr_charges_);
dp.score(spectrum, transitions, dotprod, manhattan);
dp.score(iso_cache_, spectrum, transitions, dotprod, manhattan);
}

IsotopeDistributionCache& DIAScoring::getCache()
{
return iso_cache_;
}

///////////////////////////////////////////////////////////////////////////
Expand Down Expand Up @@ -394,38 +399,13 @@ namespace OpenMS
// create the theoretical distribution from the sum formula
EmpiricalFormula empf(sum_formula);
isotope_dist = empf.getIsotopeDistribution(CoarseIsotopePatternGenerator(dia_nr_isotopes_));
iso_cache_.renormalize(isotopes, isotope_dist);
}
else
{
// create the theoretical distribution from the peptide weight
CoarseIsotopePatternGenerator solver(dia_nr_isotopes_ + 1);
isotope_dist = solver.estimateFromPeptideWeight(std::fabs(product_mz * putative_fragment_charge));
}


for (IsotopeDistribution::Iterator it = isotope_dist.begin(); it != isotope_dist.end(); ++it)
{
isotopes.intensity.push_back(it->getIntensity());
isotopes = iso_cache_.getIsotopeDistribution(std::fabs(product_mz * putative_fragment_charge));
}
isotopes.optional_begin = 0;
isotopes.optional_end = dia_nr_isotopes_;

// scale the distribution to a maximum of 1
double max = 0.0;
for (Size i = 0; i < isotopes.intensity.size(); ++i)
{
if (isotopes.intensity[i] > max)
{
max = isotopes.intensity[i];
}
}
isotopes.max = max;
for (Size i = 0; i < isotopes.intensity.size(); ++i)
{
isotopes.intensity[i] /= max;
}
isotopes.trimmed_left = 0;

// score the pattern against a theoretical one
double int_score = OpenSwath::cor_pearson(isotopes_int.begin(), isotopes_int.end(), isotopes.intensity.begin());
if (boost::math::isnan(int_score))
Expand Down
Loading