26 #include "../../exception/exceptionnotpossible.h"
38 : m_ms2MedianFilter(10), m_ms2MeanFilter(15), m_ms1MeanFilter(1)
47 [](
const double &
a,
const double &
b) { return (a < b); });
52 : m_ms2MedianFilter(other.m_ms2MedianFilter),
53 m_ms2MeanFilter(other.m_ms2MeanFilter),
54 m_ms1MeanFilter(other.m_ms1MeanFilter)
78 return *(mcsp_msrunId.get());
85 return m_ms2MedianFilter;
93 m_ms2MedianFilter = ms2MedianFilter;
100 return m_ms2MeanFilter;
108 m_ms2MeanFilter = ms2MeanFilter;
111 template <
typename T>
115 return m_ms1MeanFilter;
122 m_ms1MeanFilter = ms1MeanFilter;
126 const std::vector<MsRunRetentionTimeSeamarkPoint<T>> &
134 const std::vector<double> &
137 return m_alignedRetentionTimeVector;
144 return m_valuesCorrected;
147 const std::vector<double> &
150 return m_ms1RetentionTimeVector;
159 getCommonDeltaRt(common_points, other_seamarks);
160 return common_points;
166 std::size_t ms2_spectrum_index)
170 msp_msrunReader.get()->acquireDevice();
174 msp_msrunReader.get()->qualifiedMassSpectrum(ms2_spectrum_index,
false);
180 m_allMs2Points.push_back(ms2point);
188 double retentionTime,
189 double precursorIntensity)
195 m_allMs2Points.push_back(ms2point);
205 if(m_allMs2Points.size() == 0)
211 if(m_retentionTimeReferenceMethod ==
212 ComputeRetentionTimeReference::maximum_intensity)
216 std::sort(m_allMs2Points.begin(),
217 m_allMs2Points.end(),
219 if(a.entityHash == b.entityHash)
221 return (a.precursorIntensity > b.precursorIntensity);
223 return (
a.entityHash <
b.entityHash);
227 std::unique(m_allMs2Points.begin(),
228 m_allMs2Points.end(),
230 return (a.entityHash == b.entityHash);
233 auto it = m_allMs2Points.begin();
236 m_seamarks.push_back(
237 {it->entityHash, it->retentionTime, it->precursorIntensity});
241 m_allMs2Points.clear();
243 std::sort(m_seamarks.begin(),
247 return (a.entityHash < b.entityHash);
260 auto it = other_seamarks.begin();
264 while((it != other_seamarks.end()) &&
265 (it->entityHash < seamark.entityHash))
269 if(it == other_seamarks.end())
271 if(it->entityHash == seamark.entityHash)
274 seamark.retentionTime, seamark.retentionTime - it->retentionTime));
279 if((m_ms2MedianFilter.getHalfWindowSize() * 2 + 1) >= delta_rt.size())
282 QObject::tr(
"ERROR : MS2 alignment of MS run '%1' (%2)' not possible : "
283 "\ntoo few MS2 points (%3) in common")
284 .arg(msp_msrunReader.get()->getMsRunId().get()->getXmlId())
285 .arg(msp_msrunReader.get()->getMsRunId().get()->getFileName())
286 .arg(delta_rt.size()));
290 if((m_ms2MeanFilter.getHalfWindowSize() * 2 + 1) >= delta_rt.size())
293 QObject::tr(
"ERROR : MS2 alignment of MS run '%1' (%2)' not possible : "
294 "\ntoo few MS2 points (%3) in common")
295 .arg(msp_msrunReader.get()->getMsRunId().get()->getXmlId())
296 .arg(msp_msrunReader.get()->getMsRunId().get()->getFileName())
297 .arg(delta_rt.size()));
315 return m_alignedRetentionTimeVector.front();
317 return m_ms1RetentionTimeVector.front();
325 return m_alignedRetentionTimeVector.back();
327 return m_ms1RetentionTimeVector.back();
334 double original_retention_time)
const
336 if(m_alignedRetentionTimeVector.size() < 3)
339 QObject::tr(
"ERROR : too few aligned points to compute aligned "
340 "retention time (%1)")
341 .arg(m_ms1RetentionTimeVector.size()));
343 if(m_alignedRetentionTimeVector.size() != m_ms1RetentionTimeVector.size())
346 QObject::tr(
"ERROR : m_alignedRetentionTimeVector.size() %1 != %2 "
347 "m_ms1RetentionTimeVector.size()")
348 .arg(m_alignedRetentionTimeVector.size())
349 .arg(m_ms1RetentionTimeVector.size()));
352 std::find_if(m_ms1RetentionTimeVector.begin(),
353 m_ms1RetentionTimeVector.end(),
354 [original_retention_time](
const double &rt_point) {
355 return original_retention_time < rt_point;
357 double rt1_a, rt2_a, rt1_b, rt2_b;
358 if(it_plus == m_ms1RetentionTimeVector.end())
362 if(it_plus == m_ms1RetentionTimeVector.begin())
366 auto it_minus = it_plus - 1;
371 double ratio = (original_retention_time - rt1_a) / (rt2_a - rt1_a);
373 auto itref = m_alignedRetentionTimeVector.begin() +
374 std::distance(m_ms1RetentionTimeVector.begin(), it_minus);
380 return (((rt2_b - rt1_b) * ratio) + rt1_b);
386 double aligned_retention_time)
const
388 if(m_alignedRetentionTimeVector.size() < 3)
391 QObject::tr(
"ERROR : too few aligned points to compute aligned "
392 "retention time (%1)")
393 .arg(m_ms1RetentionTimeVector.size()));
395 if(m_alignedRetentionTimeVector.size() != m_ms1RetentionTimeVector.size())
398 QObject::tr(
"ERROR : m_alignedRetentionTimeVector.size() %1 != %2 "
399 "m_ms1RetentionTimeVector.size()")
400 .arg(m_alignedRetentionTimeVector.size())
401 .arg(m_ms1RetentionTimeVector.size()));
403 auto it_plus = std::find_if(m_alignedRetentionTimeVector.begin(),
404 m_alignedRetentionTimeVector.end(),
405 [aligned_retention_time](
const double &rt_point) {
406 return aligned_retention_time < rt_point;
408 double rt1_a, rt2_a, rt1_b, rt2_b;
409 if(it_plus == m_alignedRetentionTimeVector.end())
413 if(it_plus == m_alignedRetentionTimeVector.begin())
417 auto it_minus = it_plus - 1;
422 double ratio = (aligned_retention_time - rt1_a) / (rt2_a - rt1_a);
424 auto itref = m_ms1RetentionTimeVector.begin() +
425 std::distance(m_alignedRetentionTimeVector.begin(), it_minus);
431 return (((rt2_b - rt1_b) * ratio) + rt1_b);
435 const std::vector<MsRunRetentionTimeSeamarkPoint<T>>
438 std::vector<MsRunRetentionTimeSeamarkPoint<T>> other_seamarks = m_seamarks;
439 for(
auto &seamark : other_seamarks)
441 seamark.retentionTime =
442 translateOriginal2AlignedRetentionTime(seamark.retentionTime);
444 return other_seamarks;
451 return (m_alignedRetentionTimeVector.size() > 0);
460 std::vector<MsRunRetentionTimeSeamarkPoint<T>> other_seamarks;
461 if(msrun_retention_time_reference.
isAligned())
467 other_seamarks = msrun_retention_time_reference.
getSeamarks();
470 if((m_ms1MeanFilter.getHalfWindowSize() * 2 + 1) >=
471 m_ms1RetentionTimeVector.size())
474 QObject::tr(
"ERROR : MS1 alignment of MS run '%1' (%2)' not possible : "
475 "\ntoo few MS1 points (%3)")
476 .arg(msp_msrunReader.get()->getMsRunId().get()->getXmlId())
477 .arg(msp_msrunReader.get()->getMsRunId().get()->getFileName())
478 .arg(m_ms1RetentionTimeVector.size()));
481 qDebug() << m_seamarks[0].entityHash <<
" " << m_seamarks[0].retentionTime
482 <<
" " << other_seamarks[0].entityHash
483 << other_seamarks[0].retentionTime <<
" ";
486 getCommonDeltaRt(common_points, other_seamarks);
490 qDebug() << common_points.front().x <<
" " << common_points.front().y;
491 m_ms2MedianFilter.
filter(common_points);
493 m_ms2MeanFilter.
filter(common_points);
496 qDebug() << common_points.front().x <<
" " << common_points.front().y;
501 first_point.
x = m_ms1RetentionTimeVector.front() - (double)1;
502 if(first_point.
x < 0)
507 m_ms1RetentionTimeVector.front() -
510 common_points.push_back(first_point);
513 last_point.
x = m_ms1RetentionTimeVector.back() + 1;
514 last_point.
y = m_ms1RetentionTimeVector.back() -
516 common_points.push_back(last_point);
517 common_points.
sortX();
521 m_alignedRetentionTimeVector.clear();
523 qDebug() << common_points.front().x <<
" " << common_points.front().y;
525 Trace ms1_aligned_points;
527 linearRegressionMs2toMs1(ms1_aligned_points, common_points);
532 m_ms1MeanFilter.filter(ms1_aligned_points);
537 for(
DataPoint &data_point : ms1_aligned_points)
539 data_point.y = (data_point.x - data_point.y);
545 double correction_parameter =
546 (m_ms1RetentionTimeVector.back() - m_ms1RetentionTimeVector.front()) /
547 (ms1_aligned_points.size());
549 correction_parameter = correction_parameter / (double)4;
550 correctNewTimeValues(ms1_aligned_points, correction_parameter);
552 m_alignedRetentionTimeVector = ms1_aligned_points.
yValues();
555 return ms1_aligned_points;
562 const Trace &common_points)
566 std::vector<DataPoint>::const_iterator itms2 = common_points.begin();
567 std::vector<DataPoint>::const_iterator itms2next = itms2 + 1;
568 if(itms2next == common_points.end())
572 QObject::tr(
"ERROR : MS1 alignment of MS run '%1' (%2)' not possible : "
573 "\ntoo few common points (%3)")
574 .arg(msp_msrunReader.get()->getMsRunId().get()->getXmlId())
575 .arg(msp_msrunReader.get()->getMsRunId().get()->getFileName())
576 .arg(common_points.size()));
578 qDebug() <<
"() itms2->x=" << itms2->x <<
" itms2->y=" << itms2->y;
580 for(
double &original_rt_point : m_ms1RetentionTimeVector)
583 ms1_point.
x = original_rt_point;
585 while(ms1_point.
x > itms2next->x)
591 double ratio = (itms2next->x - itms2->x);
594 ratio = (ms1_point.
x - itms2->x) / ratio;
604 ms1_point.
y = itms2->y + ((itms2next->y - itms2->y) * ratio);
608 ms1_aligned_points.push_back(ms1_point);
615 double correction_parameter)
618 m_valuesCorrected = 0;
619 auto new_it(ms1_aligned_points.begin());
620 auto new_nextit(ms1_aligned_points.begin());
622 for(; new_nextit != ms1_aligned_points.end(); ++new_nextit, ++new_it)
624 if(new_nextit->y < new_it->y)
627 new_nextit->y = new_it->y + correction_parameter;
637 return msp_msrunReader;
643 const std::vector<double> &aligned_times)
646 if(aligned_times.size() == m_ms1RetentionTimeVector.size())
648 m_alignedRetentionTimeVector = aligned_times;
652 if(aligned_times.size() == m_ms1RetentionTimeVector.size() * 2)
654 m_alignedRetentionTimeVector = m_ms1RetentionTimeVector;
655 for(std::size_t i = 0; i < m_ms1RetentionTimeVector.size(); i++)
657 if(aligned_times[2 * i] != m_ms1RetentionTimeVector[i])
661 "ERROR : aligned_times (size=%1) vector does not have "
662 "required size (size=%2)")
663 .arg(aligned_times.size())
664 .arg(m_ms1RetentionTimeVector.size()));
666 m_alignedRetentionTimeVector[i] = aligned_times[(2 * i) + 1];
672 QObject::tr(
"ERROR : aligned_times (size=%1) vector does not have "
673 "required size (size=%2)")
674 .arg(aligned_times.size())
675 .arg(m_ms1RetentionTimeVector.size()));
687 std::vector<MsRunRetentionTimeSeamarkPoint<T>> other_seamarks;
688 if(msrun_retention_time_reference.
isAligned())
694 other_seamarks = msrun_retention_time_reference.
getSeamarks();
698 qDebug() << m_seamarks[0].entityHash <<
" " << m_seamarks[0].retentionTime
699 <<
" " << other_seamarks[0].entityHash
700 << other_seamarks[0].retentionTime <<
" ";
703 getCommonDeltaRt(common_points, other_seamarks);
704 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