Skip to content

Commit

Permalink
[FIX] remove references to removed type 'MutableEigenMatrixXdPtr'
Browse files Browse the repository at this point in the history
  • Loading branch information
@hendrikweisser
hendrikweisser committed Mar 13, 2024
1 parent 36aedc7 commit 2b6b90d
Show file tree
Hide file tree
Showing 5 changed files with 52 additions and 67 deletions.
20 changes: 9 additions & 11 deletions src/openms/include/OpenMS/ANALYSIS/QUANTITATION/IsobaricIsotopeCorrector.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -54,18 +54,18 @@ namespace OpenMS
private:
/**
* @brief Fills the input vector for the Eigen/NNLS step given the ConsensusFeature.
*
*
* Warning, assumes that the consensusMap and its ConsensusFeatures have exactly the same cardinality as the
* number of channels as in the quantitation method and are in the same order as the channels.
*
*
* I.e. for a TMT16plex, although the whole ConsensusMap has 160 potential map_index values because we had 10 files,
* every ConsensusFeature is only allowed to have exactly 16 map_index values (one for each channel) and they are
* in the same order as the channels in the quantitation method.
*
* @param[out] b
* @param[out] m_b
* @param[in] cf
* @param[in] cm
*
* @param[out] b
* @param[out] m_b
* @param[in] cf
* @param[in] cm
*/
static void fillInputVector_(Eigen::VectorXd& b,
Matrix<double>& m_b,
Expand All @@ -79,9 +79,8 @@ namespace OpenMS
*/
static void solveNNLS_(const Matrix<double>& correction_matrix,
const Matrix<double>& m_b, Matrix<double>& m_x);
// @jpfeuffer why shared_ptr needed here?
static void solveNNLS_(std::shared_ptr<Eigen::MatrixXd> & correction_matrix, std::vector<double> & b, std::vector<double> & x);
static void solveNNLS_(std::shared_ptr<const Eigen::MatrixXd> & correction_matrix, std::vector<double> & b, std::vector<double> & x);
// @jpfeuffer why shared_ptr needed here?
static void solveNNLS_(Matrix<double>::EigenMatrixType& correction_matrix, std::vector<double>& b, std::vector<double>& x);

/**
@brief
Expand Down Expand Up @@ -109,4 +108,3 @@ namespace OpenMS
const Matrix<double>& m_x);
};
} // namespace

3 changes: 1 addition & 2 deletions src/openms/include/OpenMS/MATH/MISC/NonNegativeLeastSquaresSolver.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -49,8 +49,7 @@ namespace OpenMS
@throws Exception::InvalidParameters if Matrix dimensions do not fit
*/
static Int solve(MutableEigenMatrixXdPtr & A, std::vector<double> & b, std::vector<double> & x);
static Int solve(Matrix<double>::EigenMatrixType& A, std::vector<double>& b, std::vector<double>& x);
};

} // namespace OpenMS

28 changes: 8 additions & 20 deletions src/openms/source/ANALYSIS/QUANTITATION/IsobaricIsotopeCorrector.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -26,12 +26,12 @@ namespace OpenMS
{

void
IsobaricIsotopeCorrector::correctIsotopicImpurities(std::vector<double> & intensities,
const IsobaricQuantitationMethod* quant_method)
IsobaricIsotopeCorrector::correctIsotopicImpurities(std::vector<double>& intensities,
const IsobaricQuantitationMethod* quant_method)
{
std::vector<double> res(quant_method->getNumberOfChannels());
// we need to copy anyway because NNLS will modify the input Matrix
MutableEigenMatrixXdPtr m(convertOpenMSMatrix2MutableEigenMatrixXd(quant_method->getIsotopeCorrectionMatrix()));
Matrix<double>::EigenMatrixType m = quant_method->getIsotopeCorrectionMatrix().getEigenMatrix();
solveNNLS_(m, intensities, res);
intensities = res;
}
Expand All @@ -53,7 +53,7 @@ namespace OpenMS
}
Matrix<double> m_x(c, 1.);
solveNNLS_(m, m_b, m_x);
for (Size i = 0; i < c; ++i)
{
Expand Down Expand Up @@ -84,13 +84,13 @@ namespace OpenMS

// workaround: TMT11plex has a special case where the correction matrix is the identity matrix
if (quant_method->getMethodName() != "tmt11plex")
{
{
throw Exception::InvalidParameter(__FILE__, __LINE__, OPENMS_PRETTY_FUNCTION,
"IsobaricIsotopeCorrector: The given isotope correction matrix is an identity matrix leading to no correction. "
"Please provide a valid isotope_correction matrix as it was provided with the sample kit!");
}
}

Eigen::FullPivLU<Eigen::MatrixXd> ludecomp(correction_matrix.getEigenMatrix());
Eigen::VectorXd b;
b.resize(quant_method->getNumberOfChannels());
Expand Down Expand Up @@ -134,7 +134,7 @@ namespace OpenMS
computeStats_(m_x, e_mx, cf_intensity, quant_method, stats);

/* Try this when time permits
std::vector<double> b = getIntensities_(quant_method, consensus_map_in[i], consensus_map_in);
//solve
Expand Down Expand Up @@ -207,7 +207,7 @@ namespace OpenMS
}

void
IsobaricIsotopeCorrector::solveNNLS_(std::shared_ptr<Eigen::MatrixXd> & correction_matrix, std::vector<double> & b, std::vector<double> & x)
IsobaricIsotopeCorrector::solveNNLS_(Matrix<double>::EigenMatrixType& correction_matrix, std::vector<double>& b, std::vector<double>& x)
{
Int status = NonNegativeLeastSquaresSolver::solve(correction_matrix, b, x);
if (status != NonNegativeLeastSquaresSolver::SOLVED)
Expand All @@ -216,18 +216,6 @@ namespace OpenMS
}
}

void
IsobaricIsotopeCorrector::solveNNLS_(std::shared_ptr<const Eigen::MatrixXd> & correction_matrix, std::vector<double> & b, std::vector<double> & x)
{
Eigen::MatrixXd copy = *correction_matrix;
auto copy_ptr = std::make_shared<Eigen::MatrixXd>(copy);
Int status = NonNegativeLeastSquaresSolver::solve(copy_ptr, b, x);
if (status != NonNegativeLeastSquaresSolver::SOLVED)
{
throw Exception::FailedAPICall(__FILE__, __LINE__, OPENMS_PRETTY_FUNCTION, "IsobaricIsotopeCorrector: Failed to find least-squares fit!");
}
}

void
IsobaricIsotopeCorrector::computeStats_(const std::vector<double>& m_x,
const Eigen::MatrixXd& x, const float cf_intensity,
Expand Down
18 changes: 9 additions & 9 deletions src/openms/source/MATH/MISC/NonNegativeLeastSquaresSolver.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -12,11 +12,11 @@
namespace OpenMS
{

Int NonNegativeLeastSquaresSolver::solve(MutableEigenMatrixXdPtr & A, std::vector<double> & b, std::vector<double> & x)
Int NonNegativeLeastSquaresSolver::solve(Matrix<double>::EigenMatrixType& A, std::vector<double>& b, std::vector<double>& x)
{
// this needs to be int (not Int, Size or anything else), because the external nnls constructor expects it this way!
int a_rows = (int)A->rows();
int a_cols = (int)A->cols();
int a_rows = (int)A.rows();
int a_cols = (int)A.cols();

if (a_rows != b.size())
{
Expand All @@ -26,17 +26,17 @@ namespace OpenMS
x.resize(a_cols); // description says it does not have to be initialized

// translate A to array a (column major order)
double * a_vec = A->data();
double* a_vec = A.data();

// translate b
double * b_vec = b.data();
double* b_vec = b.data();

// prepare solution array (directly copied from example)
double * x_vec = x.data();
double* x_vec = x.data();
double rnorm;
double * w = new double[a_cols];
double * zz = new double[a_rows];
int * indx = new int[a_cols];
double* w = new double[a_cols];
double* zz = new double[a_rows];
int* indx = new int[a_cols];
int mode;

NNLS::nnls_(a_vec, &a_rows, &a_rows, &a_cols, b_vec, x_vec, &rnorm, w, zz, indx, &mode);
Expand Down
50 changes: 25 additions & 25 deletions src/topp/IsobaricWorkflow.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -75,32 +75,32 @@ using namespace std;
For intensity, the closest non-zero m/z signal to the theoretical position is taken as reporter ion abundance.
The position (RT, m/z) of the consensus centroid is the precursor position in MS1 (from the MS2 spectrum);
the consensus sub-elements correspond to the theoretical channel m/z (with m/z values of 113-121 Th for iTRAQ and 126-131 Th for TMT, respectively).
For all labeling techniques, the search radius (@p reporter_mass_shift) should be set as small as possible, to avoid picking up false-positive ions as reporters.
Usually, Orbitraps deliver precision of about 0.0001 Th at this low mass range. Low intensity reporters might have a slightly higher deviation.
By default, the mass range is set to ~0.002 Th, which should be sufficient for all instruments (~15 ppm).
The tool will throw an Exception if you set it below 0.0001 Th (~0.7ppm).
The tool will also throw an Exception if you set @p reporter_mass_shift > 0.003 Th for TMT-10plex and TMT-11plex, since this could
lead to ambiguities with neighbouring channels (which are ~0.006 Th apart in most cases).
For quality control purposes, the tool reports the median distance between the theoretical vs. observed reporter ion peaks in each channel.
The search radius is fixed to 0.5 Th (regardless of the user defined search radius). This allows to track calibration issues.
For TMT-10plex, these results are automatically omitted if they could be confused with a neighbouring channel, i.e.
exceed the tolerance to a neighbouring channel with the same nominal mass (C/N channels).
If the distance is too large, you might have a m/z calibration problem (see @ref TOPP_InternalCalibration).
@note If none of the reporter ions can be detected in an MSn scan, a consensus feature will still be generated,
@note If none of the reporter ions can be detected in an MSn scan, a consensus feature will still be generated,
but the intensities of the overall feature and of all its sub-elements will be zero.
(If desired, such features can be removed by applying an intensity filter in @ref TOPP_FileFilter.)
However, if the spectrum is completely empty (no ions whatsoever), no consensus feature will be generated.
Isotope correction is done using non-negative least squares (NNLS), i.e.:@n
Minimize ||Ax - b||, subject to x >= 0, where b is the vector of observed reporter intensities (with "contaminating" isotope species),
Minimize ||Ax - b||, subject to x >= 0, where b is the vector of observed reporter intensities (with "contaminating" isotope species),
A is a correction matrix (as supplied by the manufacturer of the labeling kit) and x is the desired vector of corrected (real) reporter intensities.
Other software tools solve this problem by using an inverse matrix multiplication, but this can yield entries in x which are negative.
Other software tools solve this problem by using an inverse matrix multiplication, but this can yield entries in x which are negative.
In a real sample, this solution cannot possibly be true, so usually negative values (= negative reporter intensities) are set to zero.
However, a negative result usually means that noise was not properly accounted for in the calculation.
We thus use NNLS to get a non-negative solution, without the need to truncate negative values.
We thus use NNLS to get a non-negative solution, without the need to truncate negative values.
In the (usual) case that inverse matrix multiplication yields only positive values, our NNLS will give the exact same optimal solution.
The correction matrices can be found (and changed) in the INI file (parameter @p correction_matrix of the corresponding labeling method).
Expand All @@ -123,8 +123,8 @@ using namespace std;
After the quantitation, you may want to annotate the consensus features with corresponding peptide identifications,
obtained from an identification pipeline. Use @ref TOPP_IDMapper to perform the annotation, but make sure to set
suitably small RT and m/z tolerances for the mapping. Since the positions of the consensus features reported here
are taken from the precursor of the MS2 (also if quant was done in MS3), it should be possible to achieve a
suitably small RT and m/z tolerances for the mapping. Since the positions of the consensus features reported here
are taken from the precursor of the MS2 (also if quant was done in MS3), it should be possible to achieve a
perfect one-to-one matching of every identification (from MS2) to a single consensus feature.
Note that quantification will be solely on peptide level after this stage. In order to obtain protein quantities,
Expand Down Expand Up @@ -280,7 +280,7 @@ class TOPPIsobaricWorkflow :
} while (next_ms1_spec < exp.size());
ms1_purity = PrecursorPurity::computeInterpolatedPrecursorPurity(id_spec_idx, ms1_spec_idx, next_ms1_spec, exp, max_precursor_isotope_deviation)[0];
}

if (has_ms3)
{
id_purity = ms1_purity;
Expand All @@ -296,7 +296,7 @@ class TOPPIsobaricWorkflow :

/**
* @brief Fills an existing ConsensusFeature with all kinds of information of an identified and isobarically quantified peptide.
*
*
* @param[out] cf the ConsensusFeature to fill
* @param pep information about the PSM (object will be moved)
* @param exp the MSExperiment to extract information about spectra
Expand Down Expand Up @@ -403,7 +403,7 @@ class TOPPIsobaricWorkflow :
quant_method->setParameters(getParam_().copy(quant_method->getMethodName() + ":", true));

bool calc_id_purity = getParam_().getValue("calculate_id_purity").toBool();


Param extract_param(getParam_().copy("extraction:", true));
IsobaricChannelExtractor channel_extractor(quant_method.get());
Expand Down Expand Up @@ -444,14 +444,14 @@ class TOPPIsobaricWorkflow :
#pragma omp parallel for num_threads(max_parallel_files)
#endif
*/

for (Size i = 0; i < in_mz.size(); ++i)
{
//ConsensusMap& cur_cmap = all_cmaps[i];
ConsensusMap cur_cmap;
const String& mz_file = in_mz[i];
const String& id_file = in_id[i];

// load mzML
PeakMap exp;
mzml_file.load(mz_file, exp);
Expand Down Expand Up @@ -496,7 +496,7 @@ class TOPPIsobaricWorkflow :
{
column.second.filename = mz_file;
}

#ifdef _OPENMP
#pragma omp parallel for /*num_threads(inner_threads)*/
#endif
Expand All @@ -518,11 +518,11 @@ class TOPPIsobaricWorkflow :
{
throw Exception::InvalidValue(__FILE__, __LINE__, OPENMS_PRETTY_FUNCTION, "MS3 spectrum expected but not found.", String(exp[quant_spec_idx].getMSLevel()));
}

std::vector<double> itys = channel_extractor.extractSingleSpec(quant_spec_idx, exp, channel_qc);

// TODO if itys are all zero we can actually skip correction and quantification
auto m = convertOpenMSMatrix2MutableEigenMatrixXd(correction_matrix);
auto m = correction_matrix.getEigenMatrix();
std::vector<double> corrected(itys.size(), 0.);
NonNegativeLeastSquaresSolver::solve(m, itys, corrected);
fillConsensusFeature_(cur_cmap[pep_idx], pep, exp, id_spec_idx, quant_spec_idx, corrected, quant_method, quant_purity, id_purity, min_reporter_intensity, i);
Expand Down Expand Up @@ -558,7 +558,7 @@ class TOPPIsobaricWorkflow :
//IsobaricNormalizer::normalize(cur_cmap);

// TODO cleanup, reset?

if (cmap.empty())
{
cmap = std::move(cur_cmap);
Expand All @@ -576,7 +576,7 @@ class TOPPIsobaricWorkflow :
{
cmap.insert(cmap.end(), std::make_move_iterator(cur_cmap.begin()), std::make_move_iterator(cur_cmap.end()));
}*/

//-------------------------------------------------------------
// writing output
//-------------------------------------------------------------
Expand Down Expand Up @@ -610,17 +610,17 @@ class TOPPIsobaricWorkflow :
if (inferred_proteins.getIndistinguishableProteins().empty())
{
throw Exception::MissingInformation(
__FILE__,
__LINE__,
OPENMS_PRETTY_FUNCTION,
__FILE__,
__LINE__,
OPENMS_PRETTY_FUNCTION,
"No information on indistinguishable protein groups found.");
}

prot_quantifier.quantifyProteins(inferred_proteins);

auto const & protein_quants = prot_quantifier.getProteinResults();
if (protein_quants.empty())
{
{
OPENMS_LOG_WARN << "Warning: No proteins were quantified." << endl;
}

Expand Down

0 comments on commit 2b6b90d

Please sign in to comment.