ProteoWizard
MSnReaderTest.cpp
Go to the documentation of this file.
1 //
2 // $Id: MSnReaderTest.cpp 4129 2012-11-20 00:05:37Z chambm $
3 //
4 //
5 // Original author: Barbara Frewen <frewen@u.washington.edu>
6 //
7 // Copyright 2008 Spielberg Family Center for Applied Proteomics
8 // Cedars Sinai Medical Center, Los Angeles, California 90048
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 
26 #include "boost/filesystem/path.hpp"
27 #include "boost/filesystem/fstream.hpp"
29 #include <cstring>
30 
31 
32 using namespace pwiz::data;
33 using namespace pwiz::msdata;
34 using namespace pwiz::util;
35 namespace bfs = boost::filesystem;
36 
37 ostream* os_ = 0;
38 
40 {
41  size_t index;
43  double precursor_mz;
44  size_t num_peaks;
45  double first_peak_mz;
46  double intensity;
47  string charge_state;
48  string possible_charges; // Space delimited list
49 };
50 
52 {
53  {4, 122, 576.82, 139, 176.381, 0.9, "", "2 3"},
54  {6, 125, 785.72, 76, 333.224, 0.8, "", "2 3"},
55  {2, 120, 508.95, 82, 261.342, 0.7, "", "2 3"}
56 };
57 
58 void checkSpectrumInfo(SpectrumPtr s, size_t idx, int msLevel)
59 {
60  // Get the current spectrum info
61  SpectrumInfo* spec_info = new SpectrumInfo();
62  unit_assert(spec_info);
63  spec_info->SpectrumInfo::update(*s, true);
64 
65  // Validate spectrum info matches what we expect
66  unit_assert(spec_info->index == testSpectrum[idx].index);
67  unit_assert(spec_info->scanNumber == testSpectrum[idx].scanNumber);
68  unit_assert(spec_info->msLevel == msLevel);
69  if (msLevel > 1) {
70  unit_assert_equal(spec_info->precursors[0].mz, testSpectrum[idx].precursor_mz, 5e-2);
71  }
72  unit_assert(spec_info->data.size() == testSpectrum[idx].num_peaks);
73  unit_assert_equal(spec_info->data.at(0).mz, testSpectrum[idx].first_peak_mz, 5e-1);
74  unit_assert_equal(spec_info->data.at(0).intensity, testSpectrum[idx].intensity, 5e-1);
75 
76  // Print spectrum info
77  if (os_)
78  {
79  *os_ << "spectrum index: " << spec_info->index << "\t"
80  << "scan number: " << spec_info->scanNumber << "\t"
81  << "level: " << spec_info->msLevel << "\t";
82  if (msLevel > 1) {
83  *os_ << "precursor mz: " << spec_info->precursors[0].mz << "\t";
84  }
85  *os_ << "num peaks: " << spec_info->data.size() << "\t"
86  << "first peak mz: " << spec_info->data.at(0).mz << "\t"
87  << "intenisity: " << spec_info->data.at(0).intensity << "\t"
88  << "possible charges: ";
89  }
90 
91  if (msLevel > 1)
92  {
93  Precursor& precur = s->precursors[0];
94  SelectedIon& si = precur.selectedIons[0];
95 
96  // Since we are expecting "possible charge states", there shouldn't be a
97  // MS_charge_state!
99 
100  // Check the possible charge states (expecting 2, values = 2,3)
101  size_t numChargeStates = 0;
102  BOOST_FOREACH(CVParam& param, si.cvParams)
103  {
104  if (param.cvid == MS_possible_charge_state)
105  {
106  // Assume charge is single digit
107  unit_assert(string::npos != testSpectrum[idx].possible_charges.find(param.value));
108  numChargeStates++;
109 
110  if (os_)
111  {
112  *os_ << param.value << " ";
113  }
114  }
115  }
116  unit_assert(numChargeStates == 2);
117  }
118 
119  if (os_)
120  {
121  *os_ << "\n";
122  }
123 }
124 
125 void test(const bfs::path& datadir, int msLevel)
126 {
127  if (os_) *os_ << "test()\n";
128 
129  vector<string> filenames;
130  if (msLevel == 1) {
131  filenames.push_back("10-spec.ms1");
132  filenames.push_back("10-spec.bms1");
133  filenames.push_back("10-spec.cms1");
134  } else if (msLevel == 2) {
135  filenames.push_back("10-spec.ms2");
136  filenames.push_back("10-spec.bms2");
137  filenames.push_back("10-spec.bms2.gz");
138  filenames.push_back("10-spec.cms2");
139  } else {
140  cerr << "Invalid MS level." << endl;
141  return;
142  }
143 
144  // look up these spectrum indexes
145  size_t indexes[] = {4, 6, 2};
146  size_t num_spec = sizeof(indexes)/sizeof(size_t);
147 
148  // open each file and look up the same spec
149  vector<string>::const_iterator file_it = filenames.begin();
150  for (; file_it != filenames.end(); ++file_it)
151  {
152  string filename = *file_it;
153  MSDataFile data_file((datadir / filename).string());
154  SpectrumListPtr all_spectra = data_file.run.spectrumListPtr;
155 
156  // initialize a SpectrumInfo object with a dummy spec
157  unit_assert(data_file.run.spectrumListPtr.get());
158 
159  if (os_)
160  {
161  *os_ << "Found " << data_file.run.spectrumListPtr->size() << " spectra"
162  << " in " << filename << endl;
163  }
164 
165  // for each index, get info and print
166  for (size_t i=0; i<num_spec; i++)
167  {
168  SpectrumPtr cur_spec = all_spectra->spectrum(indexes[i], true);
169  checkSpectrumInfo(cur_spec, i, msLevel);
170  }
171 
172  // do the same thing for a list of scan numbers
173  if (os_)
174  {
175  *os_ << "Use scan numbers to get the same spec." << endl;
176  }
177 
178  size_t scan_nums[] = {122, 125, 120};
179  num_spec = sizeof(scan_nums)/sizeof(size_t);
180  for (size_t i=0; i<num_spec; i++)
181  {
182  string id_str = "scan=" + boost::lexical_cast<string>(scan_nums[i]);
183  if (os_)
184  {
185  *os_ << "Looking for the " << i << "th scan num, " << scan_nums[i]
186  << ", id " << id_str << endl;
187  }
188 
189  size_t found_index = all_spectra->find(id_str);
190  if (os_)
191  {
192  *os_ << "found_index = " << found_index << endl;
193  }
194 
195  if (found_index == all_spectra->size())
196  {
197  if (os_)
198  {
199  *os_ << "Not found." << endl;
200  }
201  }
202  else
203  {
204  SpectrumPtr cur_spec = all_spectra->spectrum(found_index, true);
205  checkSpectrumInfo(cur_spec, i, msLevel);
206  }
207  }
208  }
209 }
210 
211 int main(int argc, char* argv[])
212 {
213  TEST_PROLOG(argc, argv)
214 
215  try
216  {
217  if (argc>1 && !strcmp(argv[1],"-v")) os_ = &cout;
218 
219  std::string buildparent(argv[0]);
220  size_t pos = buildparent.find("build");
221  if (string::npos==pos)
222  {
223  buildparent = __FILE__; // nonstandard build, maybe? try using source file name
224  // something like \ProteoWizard\pwiz\pwiz\data\msdata\RAMPAdapterTest.cpp
225  pos = buildparent.rfind("pwiz");
226  }
227  buildparent.resize(pos);
228  bfs::path example_data_dir = buildparent + "example_data/";
229  test(example_data_dir, 1);
230  test(example_data_dir, 2);
231 
232  }
233  catch (exception& e)
234  {
235  TEST_FAILED(e.what())
236  }
237  catch (...)
238  {
239  TEST_FAILED("Caught unknown exception.")
240  }
241 
243 }
244 
245