ProteoWizard
SpectrumList_MZRefinerTest.cpp
Go to the documentation of this file.
1//
2// $Id$
3//
4//
5// Original author: Bryson Gibbons <bryson.gibbons@pnnl.gov>
6//
7// Copyright 2014 Pacific Northwest National Laboratory
8// Richland, WA 99352
9//
10// Licensed under the Apache License, Version 2.0 (the "License");
11// you may not use this file except in compliance with the License.
12// You may obtain a copy of the License at
13//
14// http://www.apache.org/licenses/LICENSE-2.0
15//
16// Unless required by applicable law or agreed to in writing, software
17// distributed under the License is distributed on an "AS IS" BASIS,
18// WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
19// See the License for the specific language governing permissions and
20// limitations under the License.
21//
22
23
28#include "boost/filesystem/path.hpp"
30#include <cstring>
31
32
33using namespace pwiz::cv;
34using namespace pwiz::msdata;
35using namespace pwiz::util;
36using namespace pwiz::analysis;
37namespace bfs = boost::filesystem;
38
39
40ostream* os_ = 0;
41
42// Check scan metadata and a small sample of the m/z data for high-res scans.
43void verifyScanInfo(const Spectrum& spectrum, const double& epsilon, double basePeakMZ, double lowestObservedMZ, double highestObservedMZ, int mzArrayIndex1, double mzArrayValue1, int mzArrayIndex2, double mzArrayValue2)
44{
45 unit_assert(spectrum.hasBinaryData());
46 const BinaryDataArrayPtr binaryData = spectrum.getMZArray();
47
48 if (os_)
49 {
50 *os_ << "[verifyScanInfo] " << spectrum.index << " " << spectrum.id << " "
51 << basePeakMZ << " " << lowestObservedMZ << " " << highestObservedMZ << " "
52 << mzArrayValue1 << " " << mzArrayValue2 << ": "
53 << spectrum.cvParam(MS_base_peak_m_z).value << " " << spectrum.cvParam(MS_lowest_observed_m_z).value << " "
54 << spectrum.cvParam(MS_highest_observed_m_z).value << " " << binaryData->data[mzArrayIndex1] << " "
55 << binaryData->data[mzArrayIndex2] << endl;
56 }
57
58 unit_assert_equal(spectrum.cvParam(MS_base_peak_m_z).valueAs<double>(), basePeakMZ, epsilon);
59 unit_assert_equal(spectrum.cvParam(MS_lowest_observed_m_z).valueAs<double>(), lowestObservedMZ, epsilon);
60 unit_assert_equal(spectrum.cvParam(MS_highest_observed_m_z).valueAs<double>(), highestObservedMZ, epsilon);
61 unit_assert_equal(binaryData->data[mzArrayIndex1], mzArrayValue1, epsilon);
62 unit_assert_equal(binaryData->data[mzArrayIndex2], mzArrayValue2, epsilon);
63}
64
65// Check scan precursor metadata for MS/MS scans
66void verifyPrecursorInfo(const Spectrum& spectrum, const double& epsilon, double precursorMZ, double isolationWindowTarget)
67{
68 unit_assert(!spectrum.precursors.empty());
69 const Precursor& precursor = spectrum.precursors[0];
70 unit_assert(!precursor.selectedIons.empty());
71 const SelectedIon& selectedIon = precursor.selectedIons[0];
72 unit_assert(!precursor.isolationWindow.empty());
73 const IsolationWindow& isoWindow = precursor.isolationWindow;
74
75 if (os_)
76 {
77 *os_ << "[verifyPrecursorInfo] " << spectrum.index << " " << spectrum.id << " "
78 << precursorMZ << " " << isolationWindowTarget << ": "
79 << selectedIon.cvParam(MS_selected_ion_m_z).value << " " << isoWindow.cvParam(MS_isolation_window_target_m_z).value << endl;
80 }
81
82 unit_assert_equal(selectedIon.cvParam(MS_selected_ion_m_z).valueAs<double>(), precursorMZ, epsilon);
83 unit_assert_equal(isoWindow.cvParam(MS_isolation_window_target_m_z).valueAs<double>(), isolationWindowTarget, epsilon);
84}
85
86void testShift(const bfs::path& datadir)
87{
88 MSDataFile msd((datadir / "JD_06232014_sample4_C.mzML").string());
89
90 unit_assert(msd.run.spectrumListPtr.get() && msd.run.spectrumListPtr->size()==610);
91 if (os_) *os_ << "original spectra:\n";
92 // Provided mzML file is high-res/high-res
93 double epsilon = 1e-4;
94 // MS1 scans 0, 224, 398 (0 and 224
95 verifyScanInfo(*msd.run.spectrumListPtr->spectrum(0, true), epsilon, 371.09958, 300.14306, 1568.55126, 30, 303.64633, 1200, 416.24838);
96 verifyScanInfo(*msd.run.spectrumListPtr->spectrum(224, true), epsilon, 558.30688, 301.05908, 1522.72473, 200, 407.26425, 1500, 724.32824);
97 verifyScanInfo(*msd.run.spectrumListPtr->spectrum(10, true), epsilon, 530.32782, 74.06039, 887.42852, 41, 188.11117, 93, 442.22839);
98 verifyPrecursorInfo(*msd.run.spectrumListPtr->spectrum(10), epsilon, 530.26684, 530.27);
99 verifyScanInfo(*msd.run.spectrumListPtr->spectrum(173, true), epsilon, 141.10162, 87.05542, 1187.53137, 63, 248.15817, 116, 887.44793);
100 verifyPrecursorInfo(*msd.run.spectrumListPtr->spectrum(173), epsilon, 629.30160, 629.3);
101 verifyScanInfo(*msd.run.spectrumListPtr->spectrum(346, true), epsilon, 848.45895, 116.00368, 1454.73327, 16, 185.16418, 95, 862.43109);
102 verifyPrecursorInfo(*msd.run.spectrumListPtr->spectrum(346), epsilon, 840.45480, 840.45);
103 verifyScanInfo(*msd.run.spectrumListPtr->spectrum(470, true), epsilon, 249.15857, 119.04895, 1402.77331, 23, 217.08113, 102, 1154.59863);
104 verifyPrecursorInfo(*msd.run.spectrumListPtr->spectrum(470), epsilon, 838.96706, 838.97);
105 verifyScanInfo(*msd.run.spectrumListPtr->spectrum(551, true), epsilon, 1048.55047, 155.08105, 1321.67761, 50, 368.19134, 104, 941.96954);
106 verifyPrecursorInfo(*msd.run.spectrumListPtr->spectrum(551), epsilon, 739.69935, 740.03);
107
108 shared_ptr<SpectrumList_MZRefiner> spectrumListMZRefined(
109 new SpectrumList_MZRefiner(msd, (datadir / "JD_06232014_sample4_C.mzid").string(), "specEValue", "-1e-10", IntegerSet(1, 2)));
110
111 unit_assert(spectrumListMZRefined->size() == 610);
112 if (os_) *os_ << "refined spectra:\n";
113 epsilon = 1e-2; // Increase the tolerance a little bit for the refined results (add some forgiveness with algorithm updates...)
114 verifyScanInfo(*msd.run.spectrumListPtr->spectrum(0, true), epsilon, 371.10060, 300.14388, 1568.55631, 30, 303.64715, 1200, 416.24951);
115 verifyScanInfo(*msd.run.spectrumListPtr->spectrum(224, true), epsilon, 558.30841, 301.05990, 1522.72962, 200, 407.26538, 1500, 724.33007);
116 verifyScanInfo(*msd.run.spectrumListPtr->spectrum(10, true), epsilon, 530.32928, 74.06059, 887.43126, 41, 188.11169, 93, 442.22961);
117 verifyPrecursorInfo(*msd.run.spectrumListPtr->spectrum(10), epsilon, 530.26830, 530.27145);
118 verifyScanInfo(*msd.run.spectrumListPtr->spectrum(173, true), epsilon, 141.10200, 87.05566, 1187.53519, 63, 248.15885, 116, 887.45068);
119 verifyPrecursorInfo(*msd.run.spectrumListPtr->spectrum(173), epsilon, 629.30333, 629.30172);
120 verifyScanInfo(*msd.run.spectrumListPtr->spectrum(346, true), epsilon, 848.46155, 116.00400, 1454.73795, 16, 185.16468, 95, 862.43371);
121 verifyPrecursorInfo(*msd.run.spectrumListPtr->spectrum(346), epsilon, 840.45738, 840.45257);
122 verifyScanInfo(*msd.run.spectrumListPtr->spectrum(470, true), epsilon, 249.15926, 119.04927, 1402.77782, 23, 217.08172, 102, 1154.60229);
123 verifyPrecursorInfo(*msd.run.spectrumListPtr->spectrum(470), epsilon, 838.96963, 838.97257);
124 verifyScanInfo(*msd.run.spectrumListPtr->spectrum(551, true), epsilon, 1048.55384, 155.08147, 1321.68186, 50, 368.19235, 104, 941.97253);
125 verifyPrecursorInfo(*msd.run.spectrumListPtr->spectrum(551), epsilon, 739.70141, 740.03206);
126
127 bfs::remove(datadir / "JD_06232014_sample4_C.mzRefinement.tsv");
128}
129
130
131void test(const bfs::path& datadir)
132{
133 testShift(datadir);
134}
135
136
137int main(int argc, char* argv[])
138{
139 TEST_PROLOG(argc, argv)
140
141 try
142 {
143 bfs::path datadir = ".";
144
145 // grab the parent directory for the test files.
146 for (int i=1; i<argc; i++)
147 {
148 if (!strcmp(argv[i],"-v"))
149 os_ = &cout;
150 else
151 // hack to allow running unit test from a different directory:
152 // Jamfile passes full path to specified input file.
153 // we want the path, so we can ignore filename
154 datadir = bfs::path(argv[i]).branch_path();
155 }
156
157 if (argc>1 && !strcmp(argv[1],"-v")) os_ = &cout;
158 test(datadir);
159 }
160 catch (exception& e)
161 {
162 TEST_FAILED(e.what())
163 }
164 catch (...)
165 {
166 TEST_FAILED("Caught unknown exception.")
167 }
168
170}
171
172
int main(int argc, char *argv[])
void testShift(const bfs::path &datadir)
void verifyScanInfo(const Spectrum &spectrum, const double &epsilon, double basePeakMZ, double lowestObservedMZ, double highestObservedMZ, int mzArrayIndex1, double mzArrayValue1, int mzArrayIndex2, double mzArrayValue2)
void verifyPrecursorInfo(const Spectrum &spectrum, const double &epsilon, double precursorMZ, double isolationWindowTarget)
ostream * os_
a virtual container of integers, accessible via an iterator interface, stored as union of intervals
MS_highest_observed_m_z
highest observed m/z: Highest m/z value observed in the m/z array.
Definition cv.hpp:2187
MS_lowest_observed_m_z
lowest observed m/z: Lowest m/z value observed in the m/z array.
Definition cv.hpp:2190
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_selected_ion_m_z
selected ion m/z: Mass-to-charge ratio of an selected ion.
Definition cv.hpp:2901
MS_base_peak_m_z
base peak m/z: M/z value of the signal of highest intensity in the mass spectrum.
Definition cv.hpp:2118
const double epsilon
Definition DiffTest.cpp:41
boost::shared_ptr< BinaryDataArray > BinaryDataArrayPtr
Definition MSData.hpp:417
A container that wraps DemuxWindow to preserve the full precision window boundaries.
value_type valueAs() const
templated value access with type conversion
bool empty() const
returns true iff the element contains no params or param groups
CVParam cvParam(CVID cvid) const
finds cvid in the container:
MSData object plus file I/O.
Run run
a run in mzML should correspond to a single, consecutive and coherent set of scans on an instrument.
Definition MSData.hpp:886
The method of precursor ion selection and activation.
Definition MSData.hpp:312
std::vector< SelectedIon > selectedIons
this list of precursor ions that were selected.
Definition MSData.hpp:329
IsolationWindow isolationWindow
this element captures the isolation (or 'selection') window configured to isolate one or more precurs...
Definition MSData.hpp:326
SpectrumListPtr spectrumListPtr
all mass spectra and the acquisitions underlying them are described and attached here....
Definition MSData.hpp:827
The structure that captures the generation of a peak list (including the underlying acquisitions)
Definition MSData.hpp:506
bool hasBinaryData() const
returns true iff has nonnull and nonempty BinaryDataArrayPtr
Definition MSData.hpp:535
BinaryDataArrayPtr getMZArray() const
get m/z array (may be null)
std::vector< Precursor > precursors
list and descriptions of precursors to the spectrum currently being described.
Definition MSData.hpp:520
std::string id
a unique identifier for this spectrum. It should be expected that external files may use this identif...
Definition MSData.hpp:476
size_t index
the zero-based, consecutive index of the spectrum in the SpectrumList.
Definition MSData.hpp:473
#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 TEST_PROLOG(argc, argv)
Definition unit.hpp:175