ProteoWizard
SpectrumList_3D_Test.cpp
Go to the documentation of this file.
1 //
2 // $Id: SpectrumList_3D_Test.cpp 6388 2014-06-13 17:13:45Z chambm $
3 //
4 //
5 // Original author: Matt Chambers <matt.chambers <a.t> vanderbilt.edu>
6 //
7 // Copyright 2008 Vanderbilt University - Nashville, TN 37232
8 //
9 // Licensed under the Apache License, Version 2.0 (the "License");
10 // you may not use this file except in compliance with the License.
11 // You may obtain a copy of the License at
12 //
13 // http://www.apache.org/licenses/LICENSE-2.0
14 //
15 // Unless required by applicable law or agreed to in writing, software
16 // distributed under the License is distributed on an "AS IS" BASIS,
17 // WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
18 // See the License for the specific language governing permissions and
19 // limitations under the License.
20 //
21 
27 #include "SpectrumList_3D.hpp"
28 #include "boost/foreach_field.hpp"
29 
30 using namespace pwiz::util;
31 using namespace pwiz::cv;
32 using namespace pwiz::msdata;
33 using namespace pwiz::analysis;
34 
35 ostream* os_ = 0;
36 
37 template<typename MapT>
38 typename MapT::const_iterator find_nearest(MapT const& m, typename MapT::key_type const& query, typename MapT::key_type const& tolerance)
39 {
40  typename MapT::const_iterator cur, min, max, best;
41 
42  min = m.lower_bound(query - tolerance);
43  max = m.lower_bound(query + tolerance);
44 
45  if (min == m.end() || fabs(query - min->first) > tolerance)
46  return m.end();
47  else if (min == max)
48  return min;
49  else
50  best = min;
51 
52  double minDiff = fabs(query - best->first);
53  for (cur = min; cur != max; ++cur)
54  {
55  double curDiff = fabs(query - cur->first);
56  if (curDiff < minDiff)
57  {
58  minDiff = curDiff;
59  best = cur;
60  }
61  }
62  return best;
63 }
64 
65 
66 // drift time
67 // m/z 100 200 300 400 500 600
68 // 123.4 0 1 2 1 0 0
69 // 234.5 1 2 3 2 1 0
70 // 345.6 2 3 4 3 2 1
71 // 456.7 1 2 3 2 1 0
72 // 567.8 0 1 2 1 0 0
73 // translates to:
74 // "100,123.4,0 200,123.4,1 300,123.4,2 400,123.4,1 500,123.4,0 600,123.4,0 100,234.5,1 200,123.4,2 300,123.4,3 400,123.4,2 500,123.4,1 600,123.4,0 etc."
75 
76 void test(const MSDataFile& msd, double scanStartTime, const char* driftTimeRanges, const char* expectedSpectrumTable)
77 {
78  boost::icl::interval_set<double> driftTimeRangesSet;
79  vector<string> tokens, tokens2;
80  bal::split(tokens, driftTimeRanges, bal::is_any_of(" "));
81  BOOST_FOREACH(const string& token, tokens)
82  {
83  bal::split(tokens2, token, bal::is_any_of("-"));
84  driftTimeRangesSet.add(boost::icl::continuous_interval<double>(lexical_cast<double>(tokens2[0]), lexical_cast<double>(tokens2[1])));
85  }
86 
87  Spectrum3D expectedSpectrum3d;
88  bal::split(tokens, expectedSpectrumTable, bal::is_any_of(" "));
89  BOOST_FOREACH(const string& token, tokens)
90  {
91  bal::split(tokens2, token, bal::is_any_of(","));
92  expectedSpectrum3d[lexical_cast<double>(tokens2[0])][lexical_cast<double>(tokens2[1])] = lexical_cast<double>(tokens2[2]);
93  }
94 
96  Spectrum3DPtr resultSpectrum3d = sl.spectrum3d(scanStartTime, driftTimeRangesSet);
97 
98  unit_assert(!resultSpectrum3d->empty());
99  Spectrum3D inResultButNotExpected, inExpectedButNotResult;
100  BOOST_FOREACH_FIELD((double actualDriftTime)(const Spectrum3D::value_type::second_type& resultSpectrum), *resultSpectrum3d)
101  BOOST_FOREACH_FIELD((double expectedDriftTime)(const Spectrum3D::value_type::second_type& expectedSpectrum), expectedSpectrum3d)
102  {
103  if (find_nearest(*resultSpectrum3d, expectedDriftTime, 1e-5) == resultSpectrum3d->end())
104  inExpectedButNotResult[expectedDriftTime] = Spectrum3D::value_type::second_type();
105  if (find_nearest(expectedSpectrum3d, actualDriftTime, 1e-5) == expectedSpectrum3d.end())
106  inResultButNotExpected[actualDriftTime] = Spectrum3D::value_type::second_type();
107 
108  if (fabs(actualDriftTime - expectedDriftTime) < 1e-5)
109  BOOST_FOREACH_FIELD((double actualMz)(double actualIntensity), resultSpectrum)
110  BOOST_FOREACH_FIELD((double expectedMz)(double expectedIntensity), expectedSpectrum)
111  {
112  if (find_nearest(resultSpectrum, expectedMz, 1e-4) == resultSpectrum.end())
113  inExpectedButNotResult[expectedDriftTime][expectedMz] = expectedIntensity;
114  if (find_nearest(expectedSpectrum, actualMz, 1e-4) == expectedSpectrum.end())
115  inResultButNotExpected[actualDriftTime][actualMz] = actualIntensity;
116  }
117  }
118 
119  if (os_ && !inResultButNotExpected.empty())
120  {
121  *os_ << "Extra points in the result were not expected:\n";
122  BOOST_FOREACH_FIELD((double driftTime)(const Spectrum3D::value_type::second_type& spectrum), inResultButNotExpected)
123  {
124  *os_ << driftTime << ":";
125  BOOST_FOREACH_FIELD((double mz)(double intensity), spectrum)
126  *os_ << " " << mz << "," << intensity;
127  *os_ << endl;
128  }
129  }
130 
131  if (os_ && !inExpectedButNotResult.empty())
132  {
133  *os_ << "Missing points in the result that were expected:\n";
134  BOOST_FOREACH_FIELD((double driftTime)(const Spectrum3D::value_type::second_type& spectrum), inExpectedButNotResult)
135  {
136  *os_ << driftTime << ":";
137  BOOST_FOREACH_FIELD((double mz)(double intensity), spectrum)
138  *os_ << " " << mz << "," << intensity;
139  *os_ << endl;
140  }
141  }
142 
143  unit_assert(inResultButNotExpected.empty());
144  unit_assert(inExpectedButNotResult.empty());
145 }
146 
147 
148 void parseArgs(const vector<string>& args, vector<string>& rawpaths)
149 {
150  for (size_t i = 1; i < args.size(); ++i)
151  {
152  if (args[i] == "-v") os_ = &cout;
153  else if (bal::starts_with(args[i], "--")) continue;
154  else rawpaths.push_back(args[i]);
155  }
156 }
157 
158 
159 int main(int argc, char* argv[])
160 {
161  TEST_PROLOG(argc, argv)
162 
163  try
164  {
165  vector<string> args(argv, argv+argc);
166  vector<string> rawpaths;
167  parseArgs(args, rawpaths);
168 
169  ExtendedReaderList readerList;
170 
171  BOOST_FOREACH(const string& filepath, rawpaths)
172  {
173  if (bal::ends_with(filepath, "ImsSynthCCS.d"))
174  {
175  MSDataFile msd(filepath, &readerList);
176  // scans 529, 554, 557, 615, 638
177  test(msd, 2, "29.6-29.7 33.8-33.9 34.4-34.41 34.57-34.58 44.23-44.24 48.13-48.14",
178  "33.89830,99.9954,0 34.57627,99.9954,0 29.66101,99.9954,0 34.40678,99.9954,0 44.23729,99.9954,0 48.13559,99.9954,0 " // I don't know why these isolated points show up
179  "33.89830,174.03522,0 33.89830,174.04278,1 33.89830,174.05034,1 33.89830,174.05790,1 33.89830,174.06546,1 33.89830,174.07302,0 33.89830,174.35288,0 33.89830,174.36044,1 "
180  "33.89830,174.36801,1 33.89830,174.37558,1 33.89830,174.38314,1 33.89830,174.39071,1 33.89830,174.39828,1 33.89830,174.40585,1 33.89830,174.41342,1 34.57627,177.05704,0 "
181  "34.57627,177.06466,1 34.57627,177.07229,1 34.57627,177.07991,1 34.57627,177.08754,2 34.40678,177.09517,0 34.57627,177.09517,3 34.40678,177.10279,1 34.57627,177.10279,3 "
182  "34.40678,177.11042,1 34.57627,177.11042,4 34.40678,177.11804,1 34.57627,177.11804,4 34.40678,177.12567,1 34.57627,177.12567,5 34.40678,177.13330,1 34.57627,177.13330,5 "
183  "34.40678,177.14092,1 34.57627,177.14092,5 34.40678,177.14855,1 34.57627,177.14855,4 34.40678,177.15618,1 34.57627,177.15618,4 34.40678,177.16381,0 34.57627,177.16381,3 "
184  "34.57627,177.17143,2 34.57627,177.17906,1 34.57627,177.18669,1 34.57627,177.19432,1 34.57627,177.20195,0 34.57627,177.59120,0 34.57627,177.59884,1 34.57627,177.60648,1 "
185  "34.57627,177.61411,1 34.57627,177.62175,1 34.57627,177.62939,1 34.57627,177.63703,1 34.57627,177.64466,1 34.57627,177.65230,1 34.57627,177.65994,1 34.57627,177.66758,1 "
186  "34.57627,177.67522,0 29.66101,177.71341,0 29.66101,177.72105,1 29.66101,177.72869,1 29.66101,177.73633,1 29.66101,177.74397,2 29.66101,177.75161,2 29.66101,177.75925,3 "
187  "29.66101,177.76689,4 29.66101,177.77453,4 29.66101,177.78217,4 29.66101,177.78981,4 29.66101,177.79745,4 29.66101,177.80509,4 29.66101,177.81274,3 29.66101,177.82038,3 "
188  "29.66101,177.82802,2 29.66101,177.83566,1 29.66101,177.84330,1 29.66101,177.85094,1 29.66101,177.85859,0 29.66101,178.07264,0 29.66101,178.08029,1 29.66101,178.08793,1 "
189  "29.66101,178.09558,1 29.66101,178.10323,1 29.66101,178.11088,1 34.57627,178.11088,0 29.66101,178.11852,2 34.57627,178.11852,1 29.66101,178.12617,1 34.57627,178.12617,1 "
190  "29.66101,178.13382,1 34.57627,178.13382,1 29.66101,178.14147,1 34.57627,178.14147,1 29.66101,178.14912,1 34.57627,178.14912,0 29.66101,178.15677,1 29.66101,178.16442,1 "
191  "29.66101,178.17206,0 29.66101,178.39396,0 29.66101,178.40161,1 29.66101,178.40927,1 29.66101,178.41692,1 29.66101,178.42458,1 29.66101,178.43223,2 29.66101,178.43989,2 "
192  "29.66101,178.44754,2 29.66101,178.45520,2 29.66101,178.46285,2 29.66101,178.47051,2 29.66101,178.47816,1 29.66101,178.48582,1 29.66101,178.49347,1 29.66101,178.50113,1 "
193  "29.66101,178.50879,0 29.66101,178.76154,0 29.66101,178.76920,1 29.66101,178.77686,1 29.66101,178.78452,1 29.66101,178.79219,1 29.66101,178.79985,1 29.66101,178.80751,1 "
194  "34.57627,398.99589,0 34.57627,399.00734,1 34.57627,399.01878,1 34.57627,399.03023,1 34.57627,399.04168,1 34.57627,399.05313,2 34.40678,399.06457,0 34.57627,399.06457,2 "
195  "34.40678,399.07602,1 34.57627,399.07602,3 34.40678,399.08747,1 34.57627,399.08747,4 34.40678,399.09892,1 34.57627,399.09892,5 34.40678,399.11036,1 34.57627,399.11036,5 "
196  "34.40678,399.12181,1 34.57627,399.12181,6 34.40678,399.13326,1 34.57627,399.13326,7 34.40678,399.14471,1 34.57627,399.14471,7 34.40678,399.15616,1 34.57627,399.15616,7 "
197  "34.40678,399.16761,1 34.57627,399.16761,7 34.40678,399.17906,1 34.57627,399.17906,7 34.40678,399.19051,1 34.57627,399.19051,7 34.40678,399.20196,1 34.57627,399.20196,7 "
198  "34.40678,399.21341,1 34.57627,399.21341,6 34.40678,399.22486,1 34.57627,399.22486,5 34.40678,399.23631,1 34.57627,399.23631,5 34.40678,399.24776,1 34.57627,399.24776,4 "
199  "34.40678,399.25921,1 34.57627,399.25921,3 34.57627,399.27066,2 34.57627,399.28211,2 34.57627,399.29356,1 34.57627,399.30501,1 34.57627,399.31646,1 34.57627,399.32791,1 "
200  "34.57627,399.33936,0 34.57627,400.04968,0 34.57627,400.06114,1 34.57627,400.07260,1 34.57627,400.08406,1 34.57627,400.09552,1 34.57627,400.10699,1 34.57627,400.11845,1 "
201  "34.57627,400.12991,2 34.57627,400.14138,2 34.57627,400.15284,2 34.57627,400.16430,2 34.57627,400.17577,2 34.57627,400.18723,2 34.57627,400.19869,2 34.57627,400.21016,2 "
202  "34.57627,400.22162,1 34.57627,400.23309,1 34.57627,400.24455,1 34.57627,400.25601,1 34.57627,400.26748,1 34.57627,400.27894,1 48.13559,608.09487,0 48.13559,608.10900,1 "
203  "44.23729,608.12313,0 48.13559,608.12313,1 44.23729,608.13726,1 48.13559,608.13726,1 44.23729,608.15139,1 48.13559,608.15139,1 44.23729,608.16553,1 48.13559,608.16553,1 "
204  "44.23729,608.17966,1 48.13559,608.17966,1 44.23729,608.19379,1 48.13559,608.19379,1 44.23729,608.20792,1 48.13559,608.20792,1 44.23729,608.22206,1 48.13559,608.22206,1 "
205  "44.23729,608.23619,1 48.13559,608.23619,2 44.23729,608.25032,1 48.13559,608.25032,2 44.23729,608.26445,1 48.13559,608.26445,2 44.23729,608.27859,1 48.13559,608.27859,2 "
206  "44.23729,608.29272,1 48.13559,608.29272,2 44.23729,608.30685,1 48.13559,608.30685,2 44.23729,608.32099,1 48.13559,608.32099,2 44.23729,608.33512,1 48.13559,608.33512,1 "
207  "44.23729,608.34926,1 48.13559,608.34926,1 44.23729,608.36339,1 48.13559,608.36339,1 44.23729,608.37753,1 48.13559,608.37753,1 44.23729,608.39166,1 48.13559,608.39166,1 "
208  "44.23729,608.40579,1 48.13559,608.40579,1 48.13559,608.41993,1 48.13559,608.43406,1 48.13559,608.44820,0 48.13559,609.19760,0 48.13559,609.21175,1 48.13559,609.22589,1 "
209  "48.13559,609.24004,1 48.13559,609.25418,1 48.13559,609.26832,1 48.13559,609.28247,1 48.13559,609.29661,1 48.13559,609.31076,1 48.13559,609.32491,1 48.13559,609.33905,1");
210  }
211  // TODO: add a Waters test
212  /*else if (bal::ends_with(filepath, "GFP_IMS_TEST.raw"))
213  {
214  MSDataFile msd(filepath, &readerList);
215  // scans 529, 554, 557, 615, 638
216  test(msd, 2.03500008, "13.29-13.30 13.22-13.23 13.15-13.16",
217  "13.15698,1,0 13.22623,1,0 13.29547,1,0");
218  }*/
219  }
220  }
221  catch (exception& e)
222  {
223  TEST_FAILED(e.what())
224  }
225  catch (...)
226  {
227  TEST_FAILED("Caught unknown exception.")
228  }
229 
231 }
boost::container::flat_map< double, boost::container::flat_map< double, float > > Spectrum3D
void parseArgs(const vector< string > &args, vector< string > &rawpaths)
#define TEST_EPILOG
Definition: unit.hpp:182
float lexical_cast(const std::string &str)
default ReaderList, extended to include vendor readers
Run run
a run in mzML should correspond to a single, consecutive and coherent set of scans on an instrument...
Definition: MSData.hpp:882
double mz(double neutralMass, int protonDelta, int electronDelta=0, int neutronDelta=0)
Definition: Ion.hpp:78
MSData object plus file I/O.
Definition: MSDataFile.hpp:40
boost::shared_ptr< Spectrum3D > Spectrum3DPtr
SpectrumListPtr spectrumListPtr
all mass spectra and the acquisitions underlying them are described and attached here. Subsidiary data arrays are also both described and attached here.
Definition: MSData.hpp:823
#define TEST_FAILED(x)
Definition: unit.hpp:176
ostream * os_
#define TEST_PROLOG(argc, argv)
Definition: unit.hpp:174
MapT::const_iterator find_nearest(MapT const &m, typename MapT::key_type const &query, typename MapT::key_type const &tolerance)
SpectrumList implementation that can create 3D spectra of ion mobility drift time and m/z...
int main(int argc, char *argv[])
#define unit_assert(x)
Definition: unit.hpp:85
Definition: cv.hpp:91
void test(const MSDataFile &msd, double scanStartTime, const char *driftTimeRanges, const char *expectedSpectrumTable)