timsframe.cpp
Go to the documentation of this file.
/**
* \file pappsomspp/vendors/tims/timsframe.cpp
* \date 23/08/2019
* \author Olivier Langella
* \brief handle a single Bruker's TimsTof frame
*/
/*******************************************************************************
* Copyright (c) 2019 Olivier Langella <Olivier.Langella@u-psud.fr>.
*
* This file is part of the PAPPSOms++ library.
*
* PAPPSOms++ is free software: you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation, either version 3 of the License, or
* (at your option) any later version.
*
* PAPPSOms++ is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with PAPPSOms++. If not, see <http://www.gnu.org/licenses/>.
*
******************************************************************************/
#include "timsframe.h"
#include "../../../pappsomspp/pappsoexception.h"
#include "../../../pappsomspp/exception/exceptionoutofrange.h"
#include <QDebug>
#include <QtEndian>
#include <memory>
#include <solvers.h>
#include "xicextractor/timsdirectxicextractor.h"
namespace pappso
{
TimsFrame::XicComputeStructure::XicComputeStructure(
const TimsFrame *fram_p, const XicCoordTims &xic_struct)
{
xic_ptr = xic_struct.xicSptr.get();
mobilityIndexBegin = xic_struct.scanNumBegin;
mobilityIndexEnd = xic_struct.scanNumEnd;
mzIndexLowerBound =
fram_p->getMzCalibrationInterfaceSPtr().get()->getTofIndexFromMz(
xic_struct.mzRange.lower()); // convert mz to raw digitizer value
mzIndexUpperBound =
fram_p->getMzCalibrationInterfaceSPtr().get()->getTofIndexFromMz(
xic_struct.mzRange.upper());
tmpIntensity = 0;
}
TimsFrame::TimsFrame(std::size_t timsId, quint32 scanNum)
: TimsFrameBase(timsId, scanNum)
{
// m_timsDataFrame.resize(10);
}
TimsFrame::TimsFrame(std::size_t timsId,
quint32 scanNum,
char *p_bytes,
std::size_t len)
: TimsFrameBase(timsId, scanNum)
{
// langella@themis:~/developpement/git/bruker/cbuild$
// ./src/sample/timsdataSamplePappso
// /gorgone/pappso/fichiers_fabricants/Bruker/Demo_TimsTOF_juin2019/Samples/1922001/1922001-1_S-415_Pep_Pur-1ul_Slot1-10_1_2088.d/
qDebug() << timsId;
m_timsDataFrame.resize(len);
if(p_bytes != nullptr)
{
unshufflePacket(p_bytes);
}
else
{
if(m_scanNumber == 0)
{
throw pappso::PappsoException(
QObject::tr("TimsFrame::TimsFrame(%1,%2,nullptr,%3) FAILED")
.arg(m_timsId)
.arg(m_scanNumber)
.arg(len));
}
}
}
TimsFrame::TimsFrame(const TimsFrame &other) : TimsFrameBase(other)
{
}
TimsFrame::~TimsFrame()
{
}
void
TimsFrame::unshufflePacket(const char *src)
{
qDebug();
quint64 len = m_timsDataFrame.size();
if(len % 4 != 0)
{
throw pappso::PappsoException(
QObject::tr("TimsFrame::unshufflePacket error: len%4 != 0"));
}
quint64 nb_uint4 = len / 4;
char *dest = m_timsDataFrame.data();
quint64 src_offset = 0;
for(quint64 j = 0; j < 4; j++)
{
for(quint64 i = 0; i < nb_uint4; i++)
{
dest[(i * 4) + j] = src[src_offset];
src_offset++;
}
}
qDebug();
}
std::size_t
TimsFrame::getNbrPeaks(std::size_t scanNum) const
{
if(m_timsDataFrame.size() == 0)
return 0;
/*
if(scanNum == 0)
{
quint32 res = (*(quint32 *)(m_timsDataFrame.constData() + 4)) -
(*(quint32 *)(m_timsDataFrame.constData()-4));
return res / 2;
}*/
if(scanNum == (m_scanNumber - 1))
{
auto nb_uint4 = m_timsDataFrame.size() / 4;
std::size_t cumul = 0;
for(quint32 i = 0; i < m_scanNumber; i++)
{
cumul += (*(quint32 *)(m_timsDataFrame.constData() + (i * 4)));
}
return (nb_uint4 - cumul) / 2;
}
checkScanNum(scanNum);
// quint32 *res = (quint32 *)(m_timsDataFrame.constData() + (scanNum * 4));
// qDebug() << " res=" << *res;
return (*(quint32 *)(m_timsDataFrame.constData() + ((scanNum + 1) * 4))) / 2;
}
std::size_t
TimsFrame::getScanOffset(std::size_t scanNum) const
{
std::size_t offset = 0;
for(std::size_t i = 0; i < (scanNum + 1); i++)
{
offset += (*(quint32 *)(m_timsDataFrame.constData() + (i * 4)));
}
return offset;
}
std::vector<quint32>
TimsFrame::getScanIndexList(std::size_t scanNum) const
{
qDebug();
checkScanNum(scanNum);
std::vector<quint32> scan_tof;
if(m_timsDataFrame.size() == 0)
return scan_tof;
scan_tof.resize(getNbrPeaks(scanNum));
std::size_t offset = getScanOffset(scanNum);
qint32 previous = -1;
for(std::size_t i = 0; i < scan_tof.size(); i++)
{
scan_tof[i] =
(*(quint32 *)(m_timsDataFrame.constData() + (offset * 4) + (i * 8))) +
previous;
previous = scan_tof[i];
}
qDebug();
return scan_tof;
}
std::vector<quint32>
TimsFrame::getScanIntensities(std::size_t scanNum) const
{
qDebug();
checkScanNum(scanNum);
std::vector<quint32> scan_intensities;
if(m_timsDataFrame.size() == 0)
return scan_intensities;
scan_intensities.resize(getNbrPeaks(scanNum));
std::size_t offset = getScanOffset(scanNum);
for(std::size_t i = 0; i < scan_intensities.size(); i++)
{
scan_intensities[i] = (*(quint32 *)(m_timsDataFrame.constData() +
(offset * 4) + (i * 8) + 4));
}
qDebug();
return scan_intensities;
}
void
TimsFrame::cumulateScan(std::size_t scanNum,
std::map<quint32, quint32> &accumulate_into) const
{
qDebug();
if(m_timsDataFrame.size() == 0)
return;
// checkScanNum(scanNum);
std::size_t size = getNbrPeaks(scanNum);
std::size_t offset = getScanOffset(scanNum);
qint32 previous = -1;
for(std::size_t i = 0; i < size; i++)
{
quint32 x =
(*(quint32 *)((m_timsDataFrame.constData() + (offset * 4) + (i * 8))) +
previous);
quint32 y = (*(quint32 *)(m_timsDataFrame.constData() + (offset * 4) +
(i * 8) + 4));
previous = x;
auto ret = accumulate_into.insert(std::pair<quint32, quint32>(x, y));
if(ret.second == false)
{
// already existed : cumulate
ret.first->second += y;
}
}
qDebug();
}
Trace
TimsFrame::cumulateScanToTrace(std::size_t scanNumBegin,
std::size_t scanNumEnd) const
{
qDebug();
Trace new_trace;
try
{
if(m_timsDataFrame.size() == 0)
return new_trace;
std::map<quint32, quint32> raw_spectrum;
// double local_accumulationTime = 0;
std::size_t imax = scanNumEnd + 1;
qDebug();
for(std::size_t i = scanNumBegin; i < imax; i++)
{
qDebug() << i;
cumulateScan(i, raw_spectrum);
qDebug() << i;
// local_accumulationTime += m_accumulationTime;
}
qDebug();
pappso::DataPoint data_point_cumul;
MzCalibrationInterface *mz_calibration_p =
getMzCalibrationInterfaceSPtr().get();
for(std::pair<quint32, quint32> pair_tof_intensity : raw_spectrum)
{
data_point_cumul.x =
mz_calibration_p->getMzFromTofIndex(pair_tof_intensity.first);
// normalization
data_point_cumul.y =
pair_tof_intensity.second * ((double)100.0 / m_accumulationTime);
new_trace.push_back(data_point_cumul);
}
new_trace.sortX();
qDebug();
}
catch(std::exception &error)
{
qDebug() << QString(
"Failure in TimsFrame::cumulateScanToTrace %1 to %2 :\n %3")
.arg(scanNumBegin, scanNumEnd)
.arg(error.what());
}
return new_trace;
}
void
TimsFrame::cumulateScansInRawMap(std::map<quint32, quint32> &rawSpectrum,
std::size_t scanNumBegin,
std::size_t scanNumEnd) const
{
qDebug() << "begin scanNumBegin=" << scanNumBegin
<< " scanNumEnd=" << scanNumEnd;
if(m_timsDataFrame.size() == 0)
return;
try
{
std::size_t imax = scanNumEnd + 1;
qDebug();
for(std::size_t i = scanNumBegin; i < imax; i++)
{
qDebug() << i;
cumulateScan(i, rawSpectrum);
qDebug() << i;
// local_accumulationTime += m_accumulationTime;
}
}
catch(std::exception &error)
{
qDebug() << QString("Failure in %1 %2 to %3 :\n %4")
.arg(__FUNCTION__)
.arg(scanNumBegin)
.arg(scanNumEnd)
.arg(error.what());
}
qDebug() << "end";
}
pappso::MassSpectrumCstSPtr
TimsFrame::getMassSpectrumCstSPtr(std::size_t scanNum) const
{
qDebug();
return getMassSpectrumSPtr(scanNum);
}
pappso::MassSpectrumSPtr
TimsFrame::getMassSpectrumSPtr(std::size_t scanNum) const
{
qDebug() << " scanNum=" << scanNum;
checkScanNum(scanNum);
qDebug();
pappso::MassSpectrumSPtr mass_spectrum_sptr =
std::make_shared<pappso::MassSpectrum>();
// std::vector<DataPoint>
if(m_timsDataFrame.size() == 0)
return mass_spectrum_sptr;
qDebug();
std::size_t size = getNbrPeaks(scanNum);
std::size_t offset = getScanOffset(scanNum);
MzCalibrationInterface *mz_calibration_p =
getMzCalibrationInterfaceSPtr().get();
qint32 previous = -1;
qint32 tof_index;
// std::vector<quint32> index_list;
DataPoint data_point;
for(std::size_t i = 0; i < size; i++)
{
tof_index =
(*(quint32 *)((m_timsDataFrame.constData() + (offset * 4) + (i * 8))) +
previous);
data_point.y = (*(quint32 *)(m_timsDataFrame.constData() + (offset * 4) +
(i * 8) + 4));
// intensity normalization
data_point.y *= 100.0 / m_accumulationTime;
previous = tof_index;
// mz calibration
data_point.x = mz_calibration_p->getMzFromTofIndex(tof_index);
mass_spectrum_sptr.get()->push_back(data_point);
}
qDebug();
return mass_spectrum_sptr;
}
void
TimsFrame::extractTimsXicListInRtRange(
std::vector<XicCoordTims *>::iterator &itXicListbegin,
std::vector<XicCoordTims *>::iterator &itXicListend,
XicExtractMethod method) const
{
qDebug() << std::distance(itXicListbegin, itXicListend);
std::vector<TimsFrame::XicComputeStructure> tmp_xic_list;
for(auto it = itXicListbegin; it != itXicListend; it++)
{
tmp_xic_list.push_back(TimsFrame::XicComputeStructure(this, **it));
qDebug() << " tmp_xic_struct.mobilityIndexBegin="
<< tmp_xic_list.back().mobilityIndexBegin
<< " tmp_xic_struct.mobilityIndexEnd="
<< tmp_xic_list.back().mobilityIndexEnd;
qDebug() << " tmp_xic_struct.mzIndexLowerBound="
<< tmp_xic_list.back().mzIndexLowerBound
<< " tmp_xic_struct.mzIndexUpperBound="
<< tmp_xic_list.back().mzIndexUpperBound;
}
if(tmp_xic_list.size() == 0)
return;
/*
std::sort(tmp_xic_list.begin(), tmp_xic_list.end(), [](const TimsXicStructure
&a, const TimsXicStructure &b) { return a.mobilityIndexBegin <
b.mobilityIndexBegin;
});
*/
std::vector<std::size_t> unique_scan_num_list;
for(auto &&struct_xic : tmp_xic_list)
{
for(std::size_t scan = struct_xic.mobilityIndexBegin;
(scan <= struct_xic.mobilityIndexEnd) && (scan < m_scanNumber);
scan++)
{
unique_scan_num_list.push_back(scan);
}
}
std::sort(unique_scan_num_list.begin(), unique_scan_num_list.end());
auto it_scan_num_end =
std::unique(unique_scan_num_list.begin(), unique_scan_num_list.end());
auto it_scan_num = unique_scan_num_list.begin();
while(it_scan_num != it_scan_num_end)
{
TraceSPtr ms_spectrum = getRawTraceSPtr(*it_scan_num);
// qDebug() << ms_spectrum.get()->toString();
for(auto &&tmp_xic_struct : tmp_xic_list)
{
if(((*it_scan_num) >= tmp_xic_struct.mobilityIndexBegin) &&
((*it_scan_num) <= tmp_xic_struct.mobilityIndexEnd))
{
if(method == XicExtractMethod::max)
{
tmp_xic_struct.tmpIntensity +=
ms_spectrum.get()->maxY(tmp_xic_struct.mzIndexLowerBound,
tmp_xic_struct.mzIndexUpperBound);
qDebug() << "tmp_xic_struct.tmpIntensity="
<< tmp_xic_struct.tmpIntensity;
}
else
{
// sum
tmp_xic_struct.tmpIntensity +=
ms_spectrum.get()->sumY(tmp_xic_struct.mzIndexLowerBound,
tmp_xic_struct.mzIndexUpperBound);
qDebug() << "tmp_xic_struct.tmpIntensity="
<< tmp_xic_struct.tmpIntensity;
}
}
}
it_scan_num++;
}
for(auto &&tmp_xic_struct : tmp_xic_list)
{
if(tmp_xic_struct.tmpIntensity != 0)
{
qDebug() << tmp_xic_struct.xic_ptr;
tmp_xic_struct.xic_ptr->push_back(
{m_time, tmp_xic_struct.tmpIntensity});
}
}
qDebug();
}
pappso::TraceSPtr
TimsFrame::getRawTraceSPtr(std::size_t scanNum) const
{
// qDebug();
pappso::TraceSPtr trace_sptr = std::make_shared<pappso::Trace>();
// std::vector<DataPoint>
if(m_timsDataFrame.size() == 0)
return trace_sptr;
// qDebug();
std::size_t size = getNbrPeaks(scanNum);
std::size_t offset = getScanOffset(scanNum);
qint32 previous = -1;
std::vector<quint32> index_list;
for(std::size_t i = 0; i < size; i++)
{
DataPoint data_point(
(*(quint32 *)((m_timsDataFrame.constData() + (offset * 4) + (i * 8))) +
previous),
(*(quint32 *)(m_timsDataFrame.constData() + (offset * 4) + (i * 8) +
4)));
// intensity normalization
data_point.y *= 100.0 / m_accumulationTime;
previous = data_point.x;
trace_sptr.get()->push_back(data_point);
}
// qDebug();
return trace_sptr;
}
} // namespace pappso
pappso::MzCalibrationInterface
Definition: mzcalibrationinterface.h:43
pappso::MzCalibrationInterface::getMzFromTofIndex
virtual double getMzFromTofIndex(quint32 tof_index)=0
get m/z from time of flight raw index
pappso::MzRange::lower
pappso_double lower() const
Definition: mzrange.h:71
pappso::MzRange::upper
pappso_double upper() const
Definition: mzrange.h:77
pappso::PappsoException
Definition: pappsoexception.h:42
pappso::TimsFrameBase
Definition: timsframebase.h:48
pappso::TimsFrameBase::m_accumulationTime
double m_accumulationTime
accumulation time in milliseconds
Definition: timsframebase.h:207
pappso::TimsFrameBase::m_time
double m_time
retention time
Definition: timsframebase.h:213
pappso::TimsFrameBase::m_scanNumber
quint32 m_scanNumber
total number of scans contained in this frame
Definition: timsframebase.h:197
pappso::TimsFrameBase::m_timsId
std::size_t m_timsId
Tims frame database id (the SQL identifier of this frame)
Definition: timsframebase.h:203
pappso::TimsFrameBase::getMzCalibrationInterfaceSPtr
virtual const MzCalibrationInterfaceSPtr & getMzCalibrationInterfaceSPtr() const final
get the MzCalibration model to compute mz and TOF for this frame
Definition: timsframebase.cpp:421
pappso::TimsFrameBase::checkScanNum
bool checkScanNum(std::size_t scanNum) const
check that this scan number exists
Definition: timsframebase.cpp:101
pappso::TimsFrame
Definition: timsframe.h:53
pappso::TimsFrame::m_timsDataFrame
QByteArray m_timsDataFrame
Definition: timsframe.h:181
pappso::TimsFrame::getNbrPeaks
virtual std::size_t getNbrPeaks(std::size_t scanNum) const override
get the number of peaks in this spectrum need the binary file
Definition: timsframe.cpp:131
pappso::TimsFrame::~TimsFrame
~TimsFrame()
Definition: timsframe.cpp:98
pappso::TimsFrame::cumulateScan
virtual void cumulateScan(std::size_t scanNum, std::map< quint32, quint32 > &accumulate_into) const
cumulate a scan into a map
Definition: timsframe.cpp:222
pappso::TimsFrame::TimsFrame
TimsFrame(std::size_t timsId, quint32 scanNum, char *p_bytes, std::size_t len)
Definition: timsframe.cpp:63
pappso::TimsFrame::getMassSpectrumSPtr
virtual pappso::MassSpectrumSPtr getMassSpectrumSPtr(std::size_t scanNum) const override
get Mass spectrum with peaks for this scan number need the binary file
Definition: timsframe.cpp:359
pappso::TimsFrame::unshufflePacket
void unshufflePacket(const char *src)
unshuffle data packet of tims compression type 2
Definition: timsframe.cpp:104
pappso::TimsFrame::getScanOffset
std::size_t getScanOffset(std::size_t scanNum) const
get offset for this spectrum in the binary file
Definition: timsframe.cpp:161
pappso::TimsFrame::cumulateScansInRawMap
virtual void cumulateScansInRawMap(std::map< quint32, quint32 > &rawSpectrum, std::size_t scanNumBegin, std::size_t scanNumEnd) const override
cumulate scan list into a trace into a raw spectrum map
Definition: timsframe.cpp:316
pappso::TimsFrame::cumulateScanToTrace
virtual Trace cumulateScanToTrace(std::size_t scanNumBegin, std::size_t scanNumEnd) const override
cumulate scan list into a trace
Definition: timsframe.cpp:260
pappso::TimsFrame::getScanIndexList
virtual std::vector< quint32 > getScanIndexList(std::size_t scanNum) const
get raw index list for one given scan index are not TOF nor m/z, just index on digitizer
Definition: timsframe.cpp:173
pappso::TimsFrame::getScanIntensities
virtual std::vector< quint32 > getScanIntensities(std::size_t scanNum) const
get raw intensities without transformation from one scan it needs intensity normalization
Definition: timsframe.cpp:198
pappso::TimsFrame::extractTimsXicListInRtRange
void extractTimsXicListInRtRange(std::vector< XicCoordTims * >::iterator &itXicListbegin, std::vector< XicCoordTims * >::iterator &itXicListend, XicExtractMethod method) const
Definition: timsframe.cpp:412
pappso::TimsFrame::getMassSpectrumCstSPtr
virtual pappso::MassSpectrumCstSPtr getMassSpectrumCstSPtr(std::size_t scanNum) const
get the mass spectrum corresponding to a scan number
Definition: timsframe.cpp:352
pappso::TimsFrame::getRawTraceSPtr
virtual pappso::TraceSPtr getRawTraceSPtr(std::size_t scanNum) const
get the raw index tof_index and intensities (normalized)
Definition: timsframe.cpp:504
pappso::Trace
A simple container of DataPoint instances.
Definition: trace.h:148
pappso::Trace::sortX
void sortX()
Definition: trace.cpp:905
pappso
tries to keep as much as possible monoisotopes, removing any possible C13 peaks and changes multichar...
Definition: aa.cpp:39
pappso::PeptideIonCter::y
@ y
pappso::PeptideIonCter::x
@ x
pappso::TraceSPtr
std::shared_ptr< Trace > TraceSPtr
Definition: trace.h:135
pappso::MassSpectrumCstSPtr
std::shared_ptr< const MassSpectrum > MassSpectrumCstSPtr
Definition: massspectrum.h:55
pappso::MassSpectrumSPtr
std::shared_ptr< MassSpectrum > MassSpectrumSPtr
Definition: massspectrum.h:54
pappso::XicExtractMethod
XicExtractMethod
Definition: types.h:200
pappso::XicExtractMethod::max
@ max
maximum of intensities
pappso::DataPoint
Definition: datapoint.h:21
pappso::DataPoint::x
pappso_double x
Definition: datapoint.h:22
pappso::DataPoint::y
pappso_double y
Definition: datapoint.h:23
pappso::TimsFrame::XicComputeStructure
Definition: timsframe.h:167
pappso::TimsFrame::XicComputeStructure::tmpIntensity
double tmpIntensity
Definition: timsframe.h:177
pappso::TimsFrame::XicComputeStructure::XicComputeStructure
XicComputeStructure(const TimsFrame *fram_p, const XicCoordTims &xic_struct)
Definition: timsframe.cpp:41
pappso::TimsFrame::XicComputeStructure::xic_ptr
Xic * xic_ptr
Definition: timsframe.h:172
pappso::TimsFrame::XicComputeStructure::mobilityIndexBegin
std::size_t mobilityIndexBegin
Definition: timsframe.h:173
pappso::TimsFrame::XicComputeStructure::mzIndexUpperBound
std::size_t mzIndexUpperBound
Definition: timsframe.h:176
pappso::TimsFrame::XicComputeStructure::mzIndexLowerBound
std::size_t mzIndexLowerBound
Definition: timsframe.h:175
pappso::TimsFrame::XicComputeStructure::mobilityIndexEnd
std::size_t mobilityIndexEnd
Definition: timsframe.h:174
pappso::XicCoordTims
coordinates of the XIC to extract and the resulting XIC after extraction
Definition: xiccoordtims.h:51
pappso::XicCoordTims::scanNumEnd
std::size_t scanNumEnd
mobility index end
Definition: xiccoordtims.h:91
pappso::XicCoordTims::scanNumBegin
std::size_t scanNumBegin
mobility index begin
Definition: xiccoordtims.h:87
pappso::XicCoord::xicSptr
XicSPtr xicSptr
extracted xic
Definition: xiccoord.h:113
pappso::XicCoord::mzRange
MzRange mzRange
the mass to extract
Definition: xiccoord.h:103
timsdirectxicextractor.h
minimum functions to extract XICs from Tims Data
timsframe.h
handle a single Bruker's TimsTof frame