ProteoWizard
Classes | Macros | Functions | Variables
SpectrumList_DemuxTest.cpp File Reference
#include "pwiz/utility/misc/Std.hpp"
#include "pwiz/utility/misc/Filesystem.hpp"
#include "pwiz/utility/misc/unit.hpp"
#include "pwiz/analysis/spectrum_processing/SpectrumList_Demux.hpp"
#include "pwiz/analysis/spectrum_processing/SpectrumList_PeakPicker.hpp"
#include "pwiz/analysis/demux/DemuxHelpers.hpp"
#include "pwiz/data/msdata/MSData.hpp"
#include "pwiz/data/msdata/MSDataFile.hpp"
#include "pwiz/data/msdata/Serializer_mzML.hpp"
#include "pwiz/data/msdata/Diff.hpp"
#include "pwiz_tools/common/FullReaderList.hpp"
#include <pwiz/utility/misc/IntegerSet.hpp>
#include <boost/make_shared.hpp>

Go to the source code of this file.

Classes

struct  DemuxTest
 
struct  DemuxTest::MSDPair
 

Macros

#define _VERIFY_EXACT_SPECTRUM
 

Functions

void testOverlapOnly (const string &filepath)
 
void testMSXOnly (const string &filepath)
 
void parseArgs (const vector< string > &args, vector< string > &rawpaths)
 
int main (int argc, char *argv[])
 

Variables

ostream * os_ = 0
 
const size_t TEST_SPECTRUM_OVERLAP = 134
 
const size_t TEST_SPECTRUM_OVERLAP_ORIGINAL = 67
 
const size_t NUM_DECONV_IN_TEST_SPECTRUM_OVERLAP = 2
 
const size_t TEST_SPECTRUM_OVERLAP_DEMUX_INDEX = 134
 
const size_t TEST_SPECTRUM_MSX = 105
 
const size_t TEST_SPECTRUM_MSX_ORIGINAL = 21
 
const size_t NUM_DECONV_IN_TEST_SPECTRUM_MSX = 5
 
const size_t TEST_SPECTRUM_MSX_DEMUX_INDEX = 105
 

Macro Definition Documentation

◆ _VERIFY_EXACT_SPECTRUM

#define _VERIFY_EXACT_SPECTRUM

Definition at line 34 of file SpectrumList_DemuxTest.cpp.

Function Documentation

◆ testOverlapOnly()

void testOverlapOnly ( const string &  filepath)

Definition at line 112 of file SpectrumList_DemuxTest.cpp.

113{
114 // Select the appropriate overlap demux file
115 bfs::path overlapTestFile = filepath;
116
117 // Create output file in the same directory
118 bfs::path testOutputFile = "OverlapTestOutput.mzML";
119
121
122 // Create reader for spectrum without demux
123 auto originalSpectrumList = test.GenerateSpectrumList(overlapTestFile.string());
124 SpectrumList_Demux::Params demuxParams;
125 demuxParams.optimization = DemuxOptimization::OVERLAP_ONLY;
126 auto demuxList = test.GenerateSpectrumList(overlapTestFile.string(), true, demuxParams);
127
128 // Find the original spectrum for this demux spectrum
129 auto demuxID = demuxList.spectrumList->spectrumIdentity(TEST_SPECTRUM_OVERLAP);
130 size_t originalIndex;
131 unit_assert(TryGetOriginalIndex(demuxID, originalIndex));
132
133 {
134 // Verify that the original spectrum was matched with the demux spectrum ids
135 auto originalSpectrumId = originalSpectrumList.spectrumList->spectrumIdentity(TEST_SPECTRUM_OVERLAP_ORIGINAL);
136 size_t originalIndexFromDemux;
137 unit_assert(TryGetOriginalIndex(originalSpectrumId, originalIndexFromDemux));
138 unit_assert_operator_equal(originalIndex, originalIndexFromDemux);
139 }
140
141 // Get original spectrum
142 auto originalSpectrum = originalSpectrumList.spectrumList->spectrum(TEST_SPECTRUM_OVERLAP_ORIGINAL, true);
143 auto originalMzs = originalSpectrum->getMZArray()->data;
144 auto originalIntensities = originalSpectrum->getIntensityArray()->data;
145
146 {
147 // Calculate summed intensites of the demux spectra
148 vector<double> peakSums(originalIntensities.size(), 0.0);
149 for (size_t i = 0, demuxIndex = TEST_SPECTRUM_OVERLAP; i < NUM_DECONV_IN_TEST_SPECTRUM_OVERLAP; ++i, ++demuxIndex)
150 {
151 auto demuxSpectrum = demuxList.spectrumList->spectrum(demuxIndex);
152 auto demuxIntensities = demuxSpectrum->getIntensityArray()->data;
153 auto demuxMzs = demuxSpectrum->getMZArray()->data;
154
155 vector<size_t> indexMask;
156 test.GetMask(originalMzs, demuxMzs, indexMask);
157
158 size_t j = 0;
159 for (auto index : indexMask)
160 {
161 peakSums[index] += demuxIntensities.at(j++);
162 }
163 }
164
165 // Verify that the demux spectra sum to the original spectrum
166 for (size_t i = 0; i < peakSums.size(); ++i)
167 {
168 unit_assert_equal(peakSums.at(i), originalIntensities.at(i), 1e-7);
169 }
170 }
171
172 // Verify that the spectrum window boundaries are set correctly
173 {
174 auto originalPrecursor = originalSpectrum->precursors[0];
175 double originalTarget = originalPrecursor.isolationWindow.cvParam(MS_isolation_window_target_m_z).valueAs<double>();
176 double originalLowerOffset = originalPrecursor.isolationWindow.cvParam(MS_isolation_window_lower_offset).valueAs<double>();
177 double originalUpperOffset = originalPrecursor.isolationWindow.cvParam(MS_isolation_window_upper_offset).valueAs<double>();
178 auto expectedOffset = (originalLowerOffset + originalUpperOffset) / (2.0 * static_cast<double>(NUM_DECONV_IN_TEST_SPECTRUM_OVERLAP));
179 auto windowStart = originalTarget - originalLowerOffset;
180 for (size_t i = 0, demuxIndex = TEST_SPECTRUM_OVERLAP; i < NUM_DECONV_IN_TEST_SPECTRUM_OVERLAP; ++i, ++demuxIndex)
181 {
182 double expectedTarget = windowStart + expectedOffset + 2.0 * expectedOffset * i;
183
184 auto demuxSpectrum = demuxList.spectrumList->spectrum(demuxIndex);
185 auto demuxPrecursor = demuxSpectrum->precursors[0];
186 double actualTarget = demuxPrecursor.isolationWindow.cvParam(MS_isolation_window_target_m_z).valueAs<double>();
187 double actualLowerOffset = demuxPrecursor.isolationWindow.cvParam(MS_isolation_window_lower_offset).valueAs<double>();
188 double actualUpperOffset = demuxPrecursor.isolationWindow.cvParam(MS_isolation_window_upper_offset).valueAs<double>();
189
190 // We expect the boundaries to vary based on the minimum window size. Adjacent boundaries
191 // are merged and averaged when within this window size threshold. So we only check for agreement
192 // to within this precision
193 const double minimumWindowSize = 0.01;
194
195 unit_assert_equal(expectedTarget, actualTarget, minimumWindowSize / 2.0);
196 unit_assert_equal(expectedOffset, actualLowerOffset, minimumWindowSize);
197 unit_assert_equal(expectedOffset, actualUpperOffset, minimumWindowSize);
198 }
199 }
200
201#ifdef _VERIFY_EXACT_SPECTRUM
202 // Verify that the intensity values are as expected for a demux spectrum
203
204 // TODO These are the Skyline intensities for this spectrum. It would be good to verify that they are close
205 // TODO to the actually used test values (uncommented) before removing them from the code.
206 /*vector<size_t> intensityIndices = { 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 290, 291, 292, 293, 294, 295, 296 };
207 vector<double> intensityValues =
208 {
209 0.0, 0.0, 0.0, 0.0, 4545.85, 15660.49,
210 35050.01, 56321.66, 62715.75, 43598.31, 23179.42,
211 2745.94, 3870.54, 4060.16, 3148.17, 1656.38,
212 0.0, 0.0
213 };*/
214
215 vector<double> intensityValues =
216 {
217 62715.75,
218 10856.38,
219 26514.10,
220 15964.11,
221 35976.23,
222 24815.48,
223 10131.85,
224 21044.27,
225 34393.21,
226 9127.96,
227 50067.90,
228 10287.26,
229 11103.65,
230 19305.24,
231 9583.66,
232 11572.70,
233 9995.09,
234 29599.00,
235 46296.34,
236 32724.88,
237 9292.13,
238 8167.25,
239 1111.66,
240 25497.61,
241 23860.40,
242 44635.87,
243 28415.64,
244 9848.89,
245 18376.83,
246 24337.12,
247 43483.74,
248 26286.20,
249 40075.65
250 };
251
252 auto demuxSpectrumAbsoluteCheck = demuxList.spectrumList->spectrum(TEST_SPECTRUM_OVERLAP_DEMUX_INDEX);
253 auto demuxIntensities = demuxSpectrumAbsoluteCheck->getIntensityArray()->data;
254 auto demuxMzs = demuxSpectrumAbsoluteCheck->getMZArray()->data;
255
256 // Verify intensities are equal
257 unit_assert_operator_equal(demuxIntensities.size(), intensityValues.size());
258 for (size_t i = 0; i < intensityValues.size(); ++i)
259 {
260 unit_assert_equal(demuxIntensities.at(i), intensityValues.at(i), 0.1);
261 }
262#endif
263}
const size_t TEST_SPECTRUM_OVERLAP_DEMUX_INDEX
const size_t NUM_DECONV_IN_TEST_SPECTRUM_OVERLAP
const size_t TEST_SPECTRUM_OVERLAP_ORIGINAL
const size_t TEST_SPECTRUM_OVERLAP
MS_isolation_window_lower_offset
isolation window lower offset: The extent of the isolation window in m/z below the isolation window t...
Definition cv.hpp:3183
MS_isolation_window_target_m_z
isolation window target m/z: The primary or reference m/z about which the isolation window is defined...
Definition cv.hpp:3180
MS_isolation_window_upper_offset
isolation window upper offset: The extent of the isolation window in m/z above the isolation window t...
Definition cv.hpp:3186
bool TryGetOriginalIndex(const msdata::SpectrumIdentity &spectrumIdentity, size_t &index)
Tries to read the original index of the spectrum before demultiplexing using the SpectrumIdentity of ...
MSDPair GenerateSpectrumList(const string &inputFile, bool demux=false, const SpectrumList_Demux::Params &params=SpectrumList_Demux::Params()) const
User-defined options for demultiplexing.
Optimization optimization
Optimizations can be chosen when experimental design is known.
#define unit_assert(x)
Definition unit.hpp:85
#define unit_assert_equal(x, y, epsilon)
Definition unit.hpp:99
#define unit_assert_operator_equal(expected, actual)
Definition unit.hpp:92

References DemuxTest::GenerateSpectrumList(), MS_isolation_window_lower_offset, MS_isolation_window_target_m_z, MS_isolation_window_upper_offset, NUM_DECONV_IN_TEST_SPECTRUM_OVERLAP, pwiz::analysis::SpectrumList_Demux::Params::optimization, test(), TEST_SPECTRUM_OVERLAP, TEST_SPECTRUM_OVERLAP_DEMUX_INDEX, TEST_SPECTRUM_OVERLAP_ORIGINAL, pwiz::analysis::TryGetOriginalIndex(), unit_assert, unit_assert_equal, and unit_assert_operator_equal.

Referenced by main().

◆ testMSXOnly()

void testMSXOnly ( const string &  filepath)

Definition at line 265 of file SpectrumList_DemuxTest.cpp.

266{
267 // Select the appropriate msx demux file
268 bfs::path msxTestFile = filepath;
269
270 // Create output file in the same directory
271 bfs::path testOutputFile = "MsxTestOutput.mzML";
272
274
275 // Create reader for spectrum without demux
276 auto originalSpectrumList = test.GenerateSpectrumList(msxTestFile.string());
277 SpectrumList_Demux::Params demuxParams;
278 auto demuxList = test.GenerateSpectrumList(msxTestFile.string(), true, demuxParams);
279
280 // Find the original spectrum for this demux spectrum
281 auto demuxID = demuxList.spectrumList->spectrumIdentity(TEST_SPECTRUM_MSX);
282 size_t originalIndex;
283 unit_assert(TryGetOriginalIndex(demuxID, originalIndex));
284
285 {
286 // Verify that the original spectrum was matched with the demux spectrum ids
287 auto originalSpectrumId = originalSpectrumList.spectrumList->spectrumIdentity(TEST_SPECTRUM_MSX_ORIGINAL);
288 size_t originalIndexFromDemux;
289 unit_assert(TryGetOriginalIndex(originalSpectrumId, originalIndexFromDemux));
290 unit_assert_operator_equal(originalIndex, originalIndexFromDemux);
291 }
292
293 // Get original spectrum
294 auto originalSpectrum = originalSpectrumList.spectrumList->spectrum(TEST_SPECTRUM_MSX_ORIGINAL, true);
295 auto originalMzs = originalSpectrum->getMZArray()->data;
296 auto originalIntensities = originalSpectrum->getIntensityArray()->data;
297
298 {
299 // Calculate summed intensites of the demux spectra
300 vector<double> peakSums(originalIntensities.size(), 0.0);
301 for (size_t i = 0, demuxIndex = TEST_SPECTRUM_MSX; i < NUM_DECONV_IN_TEST_SPECTRUM_MSX; ++i, ++demuxIndex)
302 {
303 auto demuxSpectrum = demuxList.spectrumList->spectrum(demuxIndex);
304 auto demuxIntensities = demuxSpectrum->getIntensityArray()->data;
305 auto demuxMzs = demuxSpectrum->getMZArray()->data;
306
307 vector<size_t> indexMask;
308 test.GetMask(originalMzs, demuxMzs, indexMask);
309
310 size_t j = 0;
311 for (auto index : indexMask)
312 {
313 peakSums[index] += demuxIntensities.at(j++);
314 }
315 }
316
317 // Verify that the demux spectra sum to the original spectrum
318 for (size_t i = 0; i < peakSums.size(); ++i)
319 {
320 unit_assert_equal(peakSums.at(i), originalIntensities.at(i), 1e-7);
321 }
322 }
323
324 // Verify that the spectrum window boundaries are set correctly
325 {
326 struct SimplePrecursor
327 {
328 double target;
329 double lowerOffset;
330 double upperOffset;
331
332 bool operator<(const SimplePrecursor& rhs) const { return this->target < rhs.target; }
333 };
334
335 vector<SimplePrecursor> originalPrecursors;
336 for (auto& precursor : originalSpectrum->precursors)
337 {
338 SimplePrecursor p;
339 p.target = precursor.isolationWindow.cvParam(MS_isolation_window_target_m_z).valueAs<double>();
340 p.lowerOffset = precursor.isolationWindow.cvParam(MS_isolation_window_lower_offset).valueAs<double>();
341 p.upperOffset = precursor.isolationWindow.cvParam(MS_isolation_window_upper_offset).valueAs<double>();
342 originalPrecursors.push_back(p);
343 }
344 sort(originalPrecursors.begin(), originalPrecursors.end());
345
346 for (size_t i = 0, demuxIndex = TEST_SPECTRUM_MSX; i < NUM_DECONV_IN_TEST_SPECTRUM_MSX; ++i, ++demuxIndex)
347 {
348 const auto& originalPrecursor = originalPrecursors.at(i);
349
350 auto demuxSpectrum = demuxList.spectrumList->spectrum(demuxIndex);
351 auto demuxPrecursor = demuxSpectrum->precursors[0];
352 double actualTarget = demuxPrecursor.isolationWindow.cvParam(MS_isolation_window_target_m_z).valueAs<double>();
353 double actualLowerOffset = demuxPrecursor.isolationWindow.cvParam(MS_isolation_window_lower_offset).valueAs<double>();
354 double actualUpperOffset = demuxPrecursor.isolationWindow.cvParam(MS_isolation_window_upper_offset).valueAs<double>();
355
356 // We expect the boundaries to vary based on the minimum window size. Adjacent boundaries
357 // are merged and averaged when within this window size threshold. So we only check for agreement
358 // to within this precision
359 const double minimumWindowSize = 0.01;
360
361 unit_assert_equal(originalPrecursor.target, actualTarget, minimumWindowSize / 2.0);
362 unit_assert_equal(originalPrecursor.lowerOffset, actualLowerOffset, minimumWindowSize);
363 unit_assert_equal(originalPrecursor.upperOffset, actualUpperOffset, minimumWindowSize);
364 }
365 }
366
367#ifdef _VERIFY_EXACT_SPECTRUM
368 // Verify that the intensity values are as expected for a demux spectrum
369
370 // TODO These are the Skyline intensities for this spectrum. It would be good to verify that they are close
371 // TODO to the actually used test values (uncommented) before removing them from the code.
372 /*vector<double> intensityValues =
373 {
374 0.0, 0.0, 0.0, 0.0, 142.95, 349.75,
375 542.87, 511.77, 248.4, 0.0, 49.28,
376 1033.65, 278.56, 0.0, 0.0, 0.0,
377 0.0, 0.0
378 };*/
379
380 vector<double> intensityValues =
381 {
382 931.31,
383 550.11,
384 650.53,
385 1870.50,
386 62.58,
387 2767.20,
388 4917.47,
389 1525.37,
390 923.80,
391 726.35,
392 1421.49,
393 1699.59,
394 3126.18,
395 25833.26,
396 23554.24,
397 10017.21,
398 900.55,
399 26146.96,
400 9478.34,
401 2643.12,
402 5988.79,
403 1562.70,
404 1952.92,
405 1392.36,
406 1354.70,
407 5745.34,
408 1891.37,
409 2545.78,
410 4131.52
411 };
412
413 auto demuxSpectrumAbsoluteCheck = demuxList.spectrumList->spectrum(TEST_SPECTRUM_MSX_DEMUX_INDEX);
414 auto demuxIntensities = demuxSpectrumAbsoluteCheck->getIntensityArray()->data;
415 auto demuxMzs = demuxSpectrumAbsoluteCheck->getMZArray()->data;
416
417 // Verify intensities are equal
418 unit_assert_operator_equal(demuxIntensities.size(), intensityValues.size());
419 for (size_t i = 0; i < intensityValues.size(); ++i)
420 {
421 unit_assert_equal(demuxIntensities.at(i), intensityValues.at(i), 0.1);
422 }
423#endif
424}
const size_t NUM_DECONV_IN_TEST_SPECTRUM_MSX
const size_t TEST_SPECTRUM_MSX
const size_t TEST_SPECTRUM_MSX_ORIGINAL
const size_t TEST_SPECTRUM_MSX_DEMUX_INDEX

References DemuxTest::GenerateSpectrumList(), MS_isolation_window_lower_offset, MS_isolation_window_target_m_z, MS_isolation_window_upper_offset, NUM_DECONV_IN_TEST_SPECTRUM_MSX, test(), TEST_SPECTRUM_MSX, TEST_SPECTRUM_MSX_DEMUX_INDEX, TEST_SPECTRUM_MSX_ORIGINAL, pwiz::analysis::TryGetOriginalIndex(), unit_assert, unit_assert_equal, and unit_assert_operator_equal.

Referenced by main().

◆ parseArgs()

void parseArgs ( const vector< string > &  args,
vector< string > &  rawpaths 
)

Definition at line 427 of file SpectrumList_DemuxTest.cpp.

428{
429 for (size_t i = 1; i < args.size(); ++i)
430 {
431 if (args[i] == "-v") os_ = &cout;
432 else if (bal::starts_with(args[i], "--")) continue;
433 else rawpaths.push_back(args[i]);
434 }
435}
ostream * os_

References os_.

Referenced by main().

◆ main()

int main ( int  argc,
char *  argv[] 
)

Definition at line 438 of file SpectrumList_DemuxTest.cpp.

439{
440 TEST_PROLOG(argc, argv)
441
442 try
443 {
444 vector<string> args(argv, argv + argc);
445 vector<string> rawpaths;
446 parseArgs(args, rawpaths);
447
448 ExtendedReaderList readerList;
449
450 BOOST_FOREACH(const string& filepath, rawpaths)
451 {
452 if (bal::ends_with(filepath, "MsxTest.mzML"))
453 {
454 testMSXOnly(filepath);
455 }
456 else if (bal::ends_with(filepath, "OverlapTest.mzML"))
457 {
458 testOverlapOnly(filepath);
459 }
460 }
461 }
462 catch (exception& e)
463 {
464 TEST_FAILED(e.what())
465 }
466 catch (...)
467 {
468 TEST_FAILED("Caught unknown exception.")
469 }
470
472}
void testMSXOnly(const string &filepath)
void testOverlapOnly(const string &filepath)
void parseArgs(const vector< string > &args, vector< string > &rawpaths)
default ReaderList, extended to include vendor readers
#define TEST_EPILOG
Definition unit.hpp:183
#define TEST_FAILED(x)
Definition unit.hpp:177
#define TEST_PROLOG(argc, argv)
Definition unit.hpp:175

References parseArgs(), TEST_EPILOG, TEST_FAILED, TEST_PROLOG, testMSXOnly(), and testOverlapOnly().

Variable Documentation

◆ os_

ostream* os_ = 0

Definition at line 40 of file SpectrumList_DemuxTest.cpp.

Referenced by parseArgs().

◆ TEST_SPECTRUM_OVERLAP

const size_t TEST_SPECTRUM_OVERLAP = 134

Definition at line 42 of file SpectrumList_DemuxTest.cpp.

Referenced by testOverlapOnly().

◆ TEST_SPECTRUM_OVERLAP_ORIGINAL

const size_t TEST_SPECTRUM_OVERLAP_ORIGINAL = 67

Definition at line 43 of file SpectrumList_DemuxTest.cpp.

Referenced by testOverlapOnly().

◆ NUM_DECONV_IN_TEST_SPECTRUM_OVERLAP

const size_t NUM_DECONV_IN_TEST_SPECTRUM_OVERLAP = 2

Definition at line 44 of file SpectrumList_DemuxTest.cpp.

Referenced by testOverlapOnly().

◆ TEST_SPECTRUM_OVERLAP_DEMUX_INDEX

const size_t TEST_SPECTRUM_OVERLAP_DEMUX_INDEX = 134

Definition at line 45 of file SpectrumList_DemuxTest.cpp.

Referenced by testOverlapOnly().

◆ TEST_SPECTRUM_MSX

const size_t TEST_SPECTRUM_MSX = 105

Definition at line 47 of file SpectrumList_DemuxTest.cpp.

Referenced by testMSXOnly().

◆ TEST_SPECTRUM_MSX_ORIGINAL

const size_t TEST_SPECTRUM_MSX_ORIGINAL = 21

Definition at line 48 of file SpectrumList_DemuxTest.cpp.

Referenced by testMSXOnly().

◆ NUM_DECONV_IN_TEST_SPECTRUM_MSX

const size_t NUM_DECONV_IN_TEST_SPECTRUM_MSX = 5

Definition at line 49 of file SpectrumList_DemuxTest.cpp.

Referenced by testMSXOnly().

◆ TEST_SPECTRUM_MSX_DEMUX_INDEX

const size_t TEST_SPECTRUM_MSX_DEMUX_INDEX = 105

Definition at line 50 of file SpectrumList_DemuxTest.cpp.

Referenced by testMSXOnly().