26 #include "../../exception/exceptionnotpossible.h"
36 : m_ms2MedianFilter(10), m_ms2MeanFilter(15), m_ms1MeanFilter(1)
45 [](
const double &
a,
const double &
b) { return (a < b); });
50 : m_ms2MedianFilter(other.m_ms2MedianFilter),
51 m_ms2MeanFilter(other.m_ms2MeanFilter),
52 m_ms1MeanFilter(other.m_ms1MeanFilter)
76 return *(mcsp_msrunId.get());
83 return m_ms2MedianFilter;
91 m_ms2MedianFilter = ms2MedianFilter;
98 return m_ms2MeanFilter;
106 m_ms2MeanFilter = ms2MeanFilter;
109 template <
typename T>
113 return m_ms1MeanFilter;
120 m_ms1MeanFilter = ms1MeanFilter;
124 const std::vector<MsRunRetentionTimeSeamarkPoint<T>> &
132 const std::vector<double> &
135 return m_alignedRetentionTimeVector;
142 return m_valuesCorrected;
145 const std::vector<double> &
148 return m_ms1RetentionTimeVector;
157 getCommonDeltaRt(common_points, other_seamarks);
158 return common_points;
164 std::size_t ms2_spectrum_index)
168 msp_msrunReader.get()->acquireDevice();
172 msp_msrunReader.get()->qualifiedMassSpectrum(ms2_spectrum_index,
false);
178 m_allMs2Points.push_back(ms2point);
186 double retentionTime,
187 double precursorIntensity)
193 m_allMs2Points.push_back(ms2point);
203 if(m_allMs2Points.size() == 0)
209 if(m_retentionTimeReferenceMethod ==
210 ComputeRetentionTimeReference::maximum_intensity)
214 std::sort(m_allMs2Points.begin(),
215 m_allMs2Points.end(),
217 if(a.entityHash == b.entityHash)
219 return (a.precursorIntensity > b.precursorIntensity);
221 return (
a.entityHash <
b.entityHash);
225 std::unique(m_allMs2Points.begin(),
226 m_allMs2Points.end(),
228 return (a.entityHash == b.entityHash);
231 auto it = m_allMs2Points.begin();
234 m_seamarks.push_back(
235 {it->entityHash, it->retentionTime, it->precursorIntensity});
239 m_allMs2Points.clear();
241 std::sort(m_seamarks.begin(),
245 return (a.entityHash < b.entityHash);
258 auto it = other_seamarks.begin();
262 while((it != other_seamarks.end()) &&
263 (it->entityHash < seamark.entityHash))
267 if(it == other_seamarks.end())
269 if(it->entityHash == seamark.entityHash)
272 seamark.retentionTime, seamark.retentionTime - it->retentionTime));
277 if((m_ms2MedianFilter.getHalfWindowSize() * 2 + 1) >= delta_rt.size())
280 QObject::tr(
"ERROR : MS2 alignment of MS run '%1' (%2)' not possible : "
281 "\ntoo few MS2 points (%3) in common")
282 .arg(msp_msrunReader.get()->getMsRunId().get()->getXmlId())
283 .arg(msp_msrunReader.get()->getMsRunId().get()->getFileName())
284 .arg(delta_rt.size()));
288 if((m_ms2MeanFilter.getHalfWindowSize() * 2 + 1) >= delta_rt.size())
291 QObject::tr(
"ERROR : MS2 alignment of MS run '%1' (%2)' not possible : "
292 "\ntoo few MS2 points (%3) in common")
293 .arg(msp_msrunReader.get()->getMsRunId().get()->getXmlId())
294 .arg(msp_msrunReader.get()->getMsRunId().get()->getFileName())
295 .arg(delta_rt.size()));
313 return m_alignedRetentionTimeVector.front();
315 return m_ms1RetentionTimeVector.front();
323 return m_alignedRetentionTimeVector.back();
325 return m_ms1RetentionTimeVector.back();
332 double original_retention_time)
const
334 if(m_alignedRetentionTimeVector.size() < 3)
337 QObject::tr(
"ERROR : too few aligned points to compute aligned "
338 "retention time (%1)")
339 .arg(m_ms1RetentionTimeVector.size()));
341 if(m_alignedRetentionTimeVector.size() != m_ms1RetentionTimeVector.size())
344 QObject::tr(
"ERROR : m_alignedRetentionTimeVector.size() %1 != %2 "
345 "m_ms1RetentionTimeVector.size()")
346 .arg(m_alignedRetentionTimeVector.size())
347 .arg(m_ms1RetentionTimeVector.size()));
350 std::find_if(m_ms1RetentionTimeVector.begin(),
351 m_ms1RetentionTimeVector.end(),
352 [original_retention_time](
const double &rt_point) {
353 return original_retention_time < rt_point;
355 double rt1_a, rt2_a, rt1_b, rt2_b;
356 if(it_plus == m_ms1RetentionTimeVector.end())
360 if(it_plus == m_ms1RetentionTimeVector.begin())
364 auto it_minus = it_plus - 1;
369 double ratio = (original_retention_time - rt1_a) / (rt2_a - rt1_a);
371 auto itref = m_alignedRetentionTimeVector.begin() +
372 std::distance(m_ms1RetentionTimeVector.begin(), it_minus);
378 return (((rt2_b - rt1_b) * ratio) + rt1_b);
384 double aligned_retention_time)
const
386 if(m_alignedRetentionTimeVector.size() < 3)
389 QObject::tr(
"ERROR : too few aligned points to compute aligned "
390 "retention time (%1)")
391 .arg(m_ms1RetentionTimeVector.size()));
393 if(m_alignedRetentionTimeVector.size() != m_ms1RetentionTimeVector.size())
396 QObject::tr(
"ERROR : m_alignedRetentionTimeVector.size() %1 != %2 "
397 "m_ms1RetentionTimeVector.size()")
398 .arg(m_alignedRetentionTimeVector.size())
399 .arg(m_ms1RetentionTimeVector.size()));
401 auto it_plus = std::find_if(m_alignedRetentionTimeVector.begin(),
402 m_alignedRetentionTimeVector.end(),
403 [aligned_retention_time](
const double &rt_point) {
404 return aligned_retention_time < rt_point;
406 double rt1_a, rt2_a, rt1_b, rt2_b;
407 if(it_plus == m_alignedRetentionTimeVector.end())
411 if(it_plus == m_alignedRetentionTimeVector.begin())
415 auto it_minus = it_plus - 1;
420 double ratio = (aligned_retention_time - rt1_a) / (rt2_a - rt1_a);
422 auto itref = m_ms1RetentionTimeVector.begin() +
423 std::distance(m_alignedRetentionTimeVector.begin(), it_minus);
429 return (((rt2_b - rt1_b) * ratio) + rt1_b);
433 const std::vector<MsRunRetentionTimeSeamarkPoint<T>>
436 std::vector<MsRunRetentionTimeSeamarkPoint<T>> other_seamarks = m_seamarks;
437 for(
auto &seamark : other_seamarks)
439 seamark.retentionTime =
440 translateOriginal2AlignedRetentionTime(seamark.retentionTime);
442 return other_seamarks;
449 return (m_alignedRetentionTimeVector.size() > 0);
458 std::vector<MsRunRetentionTimeSeamarkPoint<T>> other_seamarks;
459 if(msrun_retention_time_reference.
isAligned())
465 other_seamarks = msrun_retention_time_reference.
getSeamarks();
468 if((m_ms1MeanFilter.getHalfWindowSize() * 2 + 1) >=
469 m_ms1RetentionTimeVector.size())
472 QObject::tr(
"ERROR : MS1 alignment of MS run '%1' (%2)' not possible : "
473 "\ntoo few MS1 points (%3)")
474 .arg(msp_msrunReader.get()->getMsRunId().get()->getXmlId())
475 .arg(msp_msrunReader.get()->getMsRunId().get()->getFileName())
476 .arg(m_ms1RetentionTimeVector.size()));
479 qDebug() << m_seamarks[0].entityHash <<
" " << m_seamarks[0].retentionTime
480 <<
" " << other_seamarks[0].entityHash
481 << other_seamarks[0].retentionTime <<
" ";
484 getCommonDeltaRt(common_points, other_seamarks);
488 qDebug() << common_points.front().x <<
" " << common_points.front().y;
489 m_ms2MedianFilter.
filter(common_points);
491 m_ms2MeanFilter.
filter(common_points);
494 qDebug() << common_points.front().x <<
" " << common_points.front().y;
499 first_point.
x = m_ms1RetentionTimeVector.front() - (double)1;
500 if(first_point.
x < 0)
505 m_ms1RetentionTimeVector.front() -
508 common_points.push_back(first_point);
511 last_point.
x = m_ms1RetentionTimeVector.back() + 1;
512 last_point.
y = m_ms1RetentionTimeVector.back() -
514 common_points.push_back(last_point);
515 common_points.
sortX();
519 m_alignedRetentionTimeVector.clear();
521 qDebug() << common_points.front().x <<
" " << common_points.front().y;
523 Trace ms1_aligned_points;
525 linearRegressionMs2toMs1(ms1_aligned_points, common_points);
530 m_ms1MeanFilter.filter(ms1_aligned_points);
535 for(
DataPoint &data_point : ms1_aligned_points)
537 data_point.y = (data_point.x - data_point.y);
543 double correction_parameter =
544 (m_ms1RetentionTimeVector.back() - m_ms1RetentionTimeVector.front()) /
545 (ms1_aligned_points.size());
547 correction_parameter = correction_parameter / (double)4;
548 correctNewTimeValues(ms1_aligned_points, correction_parameter);
550 m_alignedRetentionTimeVector = ms1_aligned_points.
yValues();
553 return ms1_aligned_points;
560 const Trace &common_points)
564 std::vector<DataPoint>::const_iterator itms2 = common_points.begin();
565 std::vector<DataPoint>::const_iterator itms2next = itms2 + 1;
566 if(itms2next == common_points.end())
570 QObject::tr(
"ERROR : MS1 alignment of MS run '%1' (%2)' not possible : "
571 "\ntoo few common points (%3)")
572 .arg(msp_msrunReader.get()->getMsRunId().get()->getXmlId())
573 .arg(msp_msrunReader.get()->getMsRunId().get()->getFileName())
574 .arg(common_points.size()));
576 qDebug() <<
"() itms2->x=" << itms2->x <<
" itms2->y=" << itms2->y;
578 for(
double &original_rt_point : m_ms1RetentionTimeVector)
581 ms1_point.
x = original_rt_point;
583 while(ms1_point.
x > itms2next->x)
589 double ratio = (itms2next->x - itms2->x);
592 ratio = (ms1_point.
x - itms2->x) / ratio;
602 ms1_point.
y = itms2->y + ((itms2next->y - itms2->y) * ratio);
606 ms1_aligned_points.push_back(ms1_point);
613 double correction_parameter)
616 m_valuesCorrected = 0;
617 auto new_it(ms1_aligned_points.begin());
618 auto new_nextit(ms1_aligned_points.begin());
620 for(; new_nextit != ms1_aligned_points.end(); ++new_nextit, ++new_it)
622 if(new_nextit->y < new_it->y)
625 new_nextit->y = new_it->y + correction_parameter;
635 return msp_msrunReader;
641 const std::vector<double> &aligned_times)
644 if(aligned_times.size() == m_ms1RetentionTimeVector.size())
646 m_alignedRetentionTimeVector = aligned_times;
650 if(aligned_times.size() == m_ms1RetentionTimeVector.size() * 2)
652 m_alignedRetentionTimeVector = m_ms1RetentionTimeVector;
653 for(std::size_t i = 0; i < m_ms1RetentionTimeVector.size(); i++)
655 if(aligned_times[2 * i] != m_ms1RetentionTimeVector[i])
659 "ERROR : aligned_times (size=%1) vector does not have "
660 "required size (size=%2)")
661 .arg(aligned_times.size())
662 .arg(m_ms1RetentionTimeVector.size()));
664 m_alignedRetentionTimeVector[i] = aligned_times[(2 * i) + 1];
670 QObject::tr(
"ERROR : aligned_times (size=%1) vector does not have "
671 "required size (size=%2)")
672 .arg(aligned_times.size())
673 .arg(m_ms1RetentionTimeVector.size()));
685 std::vector<MsRunRetentionTimeSeamarkPoint<T>> other_seamarks;
686 if(msrun_retention_time_reference.
isAligned())
692 other_seamarks = msrun_retention_time_reference.
getSeamarks();
696 qDebug() << m_seamarks[0].entityHash <<
" " << m_seamarks[0].retentionTime
697 <<
" " << other_seamarks[0].entityHash
698 << other_seamarks[0].retentionTime <<
" ";
701 getCommonDeltaRt(common_points, other_seamarks);
702 return common_points;
mean filter apply mean of y values inside the window : this results in a kind of smoothing
MS run identity MsRunId identifies an MS run with a unique ID (XmlId) and contains eventually informa...
void setMs2MedianFilter(const FilterMorphoMedian &ms2MedianFilter)
std::vector< double > m_alignedRetentionTimeVector
std::vector< PeptideMs2Point > m_allMs2Points
const FilterMorphoMean & getMs1MeanFilter() const
const FilterMorphoMedian & getMs2MedianFilter() const
std::size_t m_valuesCorrected
void computeSeamarks()
convert PeptideMs2Point into Peptide seamarks this is required before computing alignment
void setMs1MeanFilter(const FilterMorphoMean &ms1MeanFilter)
Trace getCommonDeltaRt(const std::vector< MsRunRetentionTimeSeamarkPoint< T >> &other_seamarks) const
void setMs2MeanFilter(const FilterMorphoMean &ms2MeanFilter)
pappso::MsRunReaderSPtr msp_msrunReader
double getFrontRetentionTimeReference() const
MsRunRetentionTime(MsRunReaderSPtr msrun_reader_sp)
const std::vector< double > & getAlignedRetentionTimeVector() const
get aligned retention time vector
void setAlignedRetentionTimeVector(const std::vector< double > &aligned_times)
std::vector< MsRunRetentionTimeSeamarkPoint< T > > m_seamarks
const std::vector< double > & getMs1RetentionTimeVector() const
get orginal retention time vector (not aligned)
const std::vector< MsRunRetentionTimeSeamarkPoint< T > > & getSeamarks() const
ComputeRetentionTimeReference m_retentionTimeReferenceMethod
const std::vector< MsRunRetentionTimeSeamarkPoint< T > > getSeamarksReferences() const
std::vector< double > m_ms1RetentionTimeVector
const FilterMorphoMean & getMs2MeanFilter() const
pappso::MsRunReaderSPtr getMsRunReaderSPtr() const
const MsRunId & getMsRunId() const
double getBackRetentionTimeReference() const
std::size_t getNumberOfCorrectedValues() const
void addPeptideAsSeamark(const T &peptide_id, std::size_t ms2_spectrum_index)
collects all peptide evidences of a given MSrun seamarks has to be converted to peptide retention tim...
pappso::MsRunIdCstSPtr mcsp_msrunId
Class representing a fully specified mass spectrum.
pappso_double getPrecursorIntensity(bool *ok=nullptr) const
Get the intensity of the precursor ion.
pappso_double getRtInSeconds() const
Get the retention time in seconds.
A simple container of DataPoint instances.
virtual Trace & filter(const FilterInterface &filter) final
apply a filter on this trace
std::vector< pappso_double > yValues() const
tries to keep as much as possible monoisotopes, removing any possible C13 peaks and changes multichar...
std::shared_ptr< MsRunReader > MsRunReaderSPtr
double precursorIntensity