ProteoWizard
SpectrumList_DemuxTest.cpp
Go to the documentation of this file.
1//
2// $Id$
3//
4//
5// Original author: Austin Keller <atkeller .@. uw.edu>
6//
7// Licensed under the Apache License, Version 2.0 (the "License");
8// you may not use this file except in compliance with the License.
9// You may obtain a copy of the License at
10//
11// http://www.apache.org/licenses/LICENSE-2.0
12//
13// Unless required by applicable law or agreed to in writing, software
14// distributed under the License is distributed on an "AS IS" BASIS,
15// WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
16// See the License for the specific language governing permissions and
17// limitations under the License.
18//
19
30#include "pwiz_tools/common/FullReaderList.hpp"
32#include <boost/make_shared.hpp>
33
34#define _VERIFY_EXACT_SPECTRUM
35
36using namespace pwiz::msdata;
37using namespace pwiz::util;
38using namespace pwiz::analysis;
39
40ostream* os_ = 0;
41
42const size_t TEST_SPECTRUM_OVERLAP = 134;
46
47const size_t TEST_SPECTRUM_MSX = 105;
51
52struct DemuxTest {
53 typedef boost::shared_ptr<MSData> MSDataPtr;
60
61 MSDPair GenerateSpectrumList(const string& inputFile,
62 bool demux = false,
64
65 void GetMask(const vector<double>& original, const vector<double>& derived, vector<size_t>& mask) const;
66};
67
69 bool demux,
70 const SpectrumList_Demux::Params& params) const
71{
72 FullReaderList readers;
73 MSDataPtr msdPtr = boost::make_shared<MSDataFile>(inputFile, &readers);
74 IntegerSet levelsToCentroid(1, 2);
75 SpectrumListPtr centroidedPtr(
76 new SpectrumList_PeakPicker(msdPtr->run.spectrumListPtr,
77 PeakDetectorPtr(boost::make_shared<LocalMaximumPeakDetector>(3)),
78 true,
79 levelsToCentroid));
80 msdPtr->filterApplied();
81
82 if (!demux)
83 return MSDPair(msdPtr, centroidedPtr);
84
85 SpectrumListPtr demuxList(new SpectrumList_Demux(centroidedPtr, params));
86 msdPtr->filterApplied();
87 msdPtr->run.spectrumListPtr = demuxList;
88 return MSDPair(msdPtr, demuxList);
89}
90
91void DemuxTest::GetMask(const vector<double>& original, const vector<double>& derived, vector<size_t>& mask) const
92{
93 unit_assert(std::is_sorted(original.begin(), original.end()));
94 unit_assert(std::is_sorted(derived.begin(), derived.end()));
95 mask.clear();
96 auto originalIt = original.begin();
97 for (auto derivedIt = derived.begin(); derivedIt != derived.end(); ++derivedIt)
98 {
99 for (; originalIt != original.end(); ++originalIt)
100 {
101 if (abs(*originalIt - *derivedIt) < 1.0e-5)
102 {
103 mask.push_back(originalIt - original.begin());
104 break;
105 }
106 unit_assert(*originalIt < *derivedIt);
107 }
108 }
109 unit_assert_operator_equal(derived.size(), mask.size());
110}
111
112void testOverlapOnly(const string& filepath)
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}
264
265void testMSXOnly(const string& filepath)
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}
425
426
427void parseArgs(const vector<string>& args, vector<string>& rawpaths)
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}
436
437
438int main(int argc, char* argv[])
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}
Helper functions for demultiplexing Helper functions include nice methods of accessing CV parameters ...
int main(int argc, char *argv[])
const size_t TEST_SPECTRUM_OVERLAP_DEMUX_INDEX
const size_t NUM_DECONV_IN_TEST_SPECTRUM_MSX
const size_t TEST_SPECTRUM_MSX
const size_t TEST_SPECTRUM_MSX_ORIGINAL
void testMSXOnly(const string &filepath)
const size_t NUM_DECONV_IN_TEST_SPECTRUM_OVERLAP
const size_t TEST_SPECTRUM_OVERLAP_ORIGINAL
const size_t TEST_SPECTRUM_MSX_DEMUX_INDEX
ostream * os_
const size_t TEST_SPECTRUM_OVERLAP
void testOverlapOnly(const string &filepath)
void parseArgs(const vector< string > &args, vector< string > &rawpaths)
SpectrumList decorator implementation that can demultiplex spectra of several precursor windows acqui...
SpectrumList implementation to replace peak profiles with picked peaks.
default ReaderList, extended to include vendor readers
a virtual container of integers, accessible via an iterator interface, stored as union of intervals
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
boost::shared_ptr< PeakDetector > PeakDetectorPtr
bool TryGetOriginalIndex(const msdata::SpectrumIdentity &spectrumIdentity, size_t &index)
Tries to read the original index of the spectrum before demultiplexing using the SpectrumIdentity of ...
boost::shared_ptr< SpectrumList > SpectrumListPtr
Definition MSData.hpp:711
MSDPair(MSDataPtr msdata, SpectrumListPtr spectrumList)
MSDPair GenerateSpectrumList(const string &inputFile, bool demux=false, const SpectrumList_Demux::Params &params=SpectrumList_Demux::Params()) const
void GetMask(const vector< double > &original, const vector< double > &derived, vector< size_t > &mask) const
boost::shared_ptr< MSData > MSDataPtr
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 TEST_EPILOG
Definition unit.hpp:183
#define TEST_FAILED(x)
Definition unit.hpp:177
#define unit_assert_equal(x, y, epsilon)
Definition unit.hpp:99
#define unit_assert_operator_equal(expected, actual)
Definition unit.hpp:92
#define TEST_PROLOG(argc, argv)
Definition unit.hpp:175