3 // Copyright (C) 2006-2007, Stowers Institute for Medical Research
5 // This program is free software; you can redistribute it and/or modify
6 // it under the terms of the GNU General Public License as published by
7 // the Free Software Foundation; either version 2 of the License, or
8 // (at your option) any later version.
10 // This program is distributed in the hope that it will be useful,
11 // but WITHOUT ANY WARRANTY; without even the implied warranty of
12 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
13 // GNU General Public License for more details.
15 // You should have received a copy of the GNU General Public License along
16 // with this program; if not, write to the Free Software Foundation, Inc.,
17 // 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
36 class mass_regime_parameters
{
42 // residue char (incl '['/']') -> mass (per regime) with any fixed mod
43 std::vector
<double> fixed_residue_mass
;
46 // FIX: want these class members to all be static, but had SWIG trouble. The
47 // member 'the' is a handle to the singleton instance, for now.
50 bool estimate_only
; // just estimate work required
51 bool show_progress
; // running progress on stderr
53 std::vector
<double> ln_factorial
; // factorial[n] == ln(n!)
55 // deuterium not (yet) implemented
58 std::vector
<mass_regime_parameters
> parent_mass_regime
;
59 std::vector
<mass_regime_parameters
> fragment_mass_regime
;
62 double parent_mass_tolerance_1
; // for +1
63 double parent_mass_tolerance_max
; // for +N (typically +3)
65 unsigned int minimum_peptide_length
;
67 double fragment_mass_tolerance
;
68 int intensity_class_count
;
70 static parameters the
;
76 int sequence_index
; // index of this sequence in the database FIX:nn
77 int sequence_offset
; // offset of this run's start in sequence
79 std::vector
<int> cleavage_points
;
80 std::string name
; // e.g., initial defline word
83 sequence_run(const int sequence_index
, const int sequence_offset
,
84 const std::string sequence
, std::vector
<int> cleavage_points
,
85 const std::string name
)
86 : sequence_index(sequence_index
), sequence_offset(sequence_offset
),
87 sequence(sequence
), cleavage_points(cleavage_points
),
93 // information about the search that's passed in to--but not changed by--the
95 class search_context
{
97 unsigned int mod_count
;
98 int mass_regime_index
;
100 std::string pca_residues
;
103 double N_delta
; // 0 if none
104 double C_delta
; // 0 if none
106 // maps residue char to list of indices for _count/_delta
107 std::vector
< std::vector
<int> > delta_bag_lookup
;
108 std::vector
<double> delta_bag_delta
;
109 std::vector
<int> delta_bag_count
;
111 double parent_fixed_mass
;
112 double fragment_N_fixed_mass
;
113 double fragment_C_fixed_mass
;
115 // general search parameters
116 int maximum_missed_cleavage_sites
;
117 bool nonspecific_cleavage
;
119 std::vector
<sequence_run
> sequence_runs
;
127 int intensity_class
; // 0..N; -1 == unclassified
129 explicit peak(double mz
=0, double intensity
=0, int intensity_class
=-1)
130 : mz(mz
), intensity(intensity
), intensity_class(intensity_class
) {
133 char *__repr__() const;
135 static bool less_mz(peak x
, peak y
) { return x
.mz
< y
.mz
; }
137 static bool greater_intensity(peak x
, peak y
) {
138 return x
.intensity
> y
.intensity
;
141 static double get_mz(double mass
, int charge
) {
142 const parameters
&CP
= parameters::the
;
143 return mass
/charge
+ CP
.proton_mass
;
148 // A spectrum searched on multiple charges will be turned into separate
149 // spectra having differing id's, but the same physical id.
153 double mass
; // [M + H+], aka "neutral mass" (protonated)
155 static const int MAX_SUPPORTED_CHARGE
= 15;
156 static const int MAX_THEORETICAL_FRAGMENTS
= 256;
157 std::vector
<peak
> peaks
; // always ordered by increasing m/z!
159 std::vector
<int> intensity_class_counts
; // number of peaks in each class
161 // This is the mz range of peaks, but with a buffer on each end of size
162 // fragment_mass_tolerance.
166 // These are caches of the data in peaks above. They speed up search.
167 double *peak_mz_cache
;
168 int *peak_intensity_class_cache
;
170 double total_ion_current
;
172 // This is the total number of bins, 2*fragment_tolerance mz wide in peak
173 // range. So, for example, if the fragment tolerance is 1 and the min and
174 // max mz are 200 and 1800, total peak bins would be 800.
176 int empty_peak_bins
; // number of bins that aren't occupied by a peak
179 int file_id
; // file # of spectra's source
180 int id
; // unique for all spectra
181 int physical_id
; // FIX: unneeded?
183 unsigned long comparisons
; // times scored against theoretical spectra
185 // Construct an empty spectrum.
186 explicit spectrum(double mass
=0,
187 int charge
=0) : mass(mass
), charge(charge
), min_peak_mz(0),
188 max_peak_mz(0), peak_mz_cache(0),
189 peak_intensity_class_cache(0),
190 total_ion_current(0),
191 total_peak_bins(0), empty_peak_bins(0),
192 file_id(-1), physical_id(-1),
195 if (charge
> MAX_SUPPORTED_CHARGE
)
196 throw std::invalid_argument("attempt to create a spectrum with greater"
197 " than supported charge");
201 char *__repr__() const;
204 void set_peaks_from_matrix(const std::vector
< std::vector
<double> > &m
) {
205 peaks
.resize(m
.size());
206 std::vector
<peak
>::iterator p_it
= peaks
.begin();
207 for (std::vector
< std::vector
<double> >::const_iterator it
= m
.begin();
208 it
!= m
.end(); it
++, p_it
++) {
210 throw std::invalid_argument("invalid matrix (must be size N x 2)");
212 p_it
->intensity
= (*it
)[1];
217 // Read spectra from file in ms2 format, tagging them with file_id. Before
218 // reading, seek to absolute position offset_begin. If offset_end != -1,
219 // any spectra that begin after position offset_end in the file are not
221 static std::vector
<spectrum
>
222 read_spectra_from_ms2(FILE *f
, const int file_id
, const long offset_begin
,
223 const long offset_end
);
226 // Filter peaks to limit their number according to TIC_cutoff_proportion,
227 // and to remove those too large to be fragment products. (Peaks are
228 // assumed to be ordered by increasing mz.)
229 void filter_peaks(double TIC_cutoff_proportion
,
230 double parent_mass_tolerance_max
);
232 // Classify peaks and update class_counts.
233 void classify(int intensity_class_count
, double intensity_class_ratio
,
234 double fragment_mass_tolerance
);
236 // Update the *_cache fields from the peaks field.
237 void update_peak_cache();
239 void clear_peak_cache() {
241 //delete[] peak_mz_cache;
244 if (peak_intensity_class_cache
) {
245 //delete[] peak_intensity_class_cache;
246 peak_intensity_class_cache
= 0;
250 // Strictly speaking, we should release the cache on destruction. This
251 // would mean implementing bug-prone copy constructor and assignment
252 // constructor functions. Instead, skip it. This leaks memory if
253 // update_peak_cache() is called multiple times, but we don't do that.
255 // ~spectrum() { clear_peak_cache(); }
257 // Store the list of spectra that search_peptide will search against, and
258 // also build spectrum_mass_index.
259 static void set_searchable_spectra(const std::vector
<spectrum
> &spectra
);
261 // Search sequence runs for matches according to the context against the
262 // spectra. Updates score_stats and the number of candidate spectra found.
263 static void search_runs(const search_context
&context
, score_stats
&stats
);
265 // conceptually these are 'protected:', but we're taking it easy here
267 void set_id() { id
= next_id
++; assert(next_id
> 0); }
270 static int next_physical_id
;
272 static std::vector
<spectrum
> searchable_spectra
;
273 // parent mass -> searchable_spectra index
274 static std::multimap
<double,
275 std::vector
<spectrum
>::size_type
> spectrum_mass_index
;
279 struct mass_trace_item
{
284 bool operator==(const mass_trace_item
&x
) const {
285 return this->position
== x
.position
and this->delta
== x
.delta
;
290 // This is everything we want to remember about a match, so that we can report
296 std::string peptide_sequence
;
297 double predicted_parent_mass
;
299 int mass_regime_index
;
302 std::vector
<mass_trace_item
> mass_trace
;
304 // these are vectors of results for storing peptides found in multiple
305 // sequence locations
306 std::vector
<int> peptide_begin
; // absolute position within locus
307 std::vector
<std::string
> sequence_name
;
309 match() : score(0), spectrum_index(-1), predicted_parent_mass(0),
310 mass_regime_index(-1), conjunct_index(-1), pca_delta(0) { }
314 // This holds information about the best peptide/spectrum matches found, and
315 // some statistics. Essentially, it holds the results of the search process.
318 score_stats(int spectrum_count
, int best_result_count
)
319 : candidate_spectrum_count(0), combinations_searched(0) {
320 best_matches
.resize(spectrum_count
);
321 assert(best_result_count
>= 1);
322 for (int i
=0; i
<spectrum_count
; i
++)
323 best_matches
[i
].resize(best_result_count
);
326 // spectrum index -> list of match (only top N matches per spectrum) The
327 // best match lists are of fixed length (typically 5), and ordered
328 // best-first. Trailing entries with score 0 are null matches (to simplify
330 std::vector
< std::vector
<match
> > best_matches
;
333 // How many times was a spectrum scored against a peptide?
334 unsigned long long candidate_spectrum_count
; // may be > 2^32
335 // How many different mod combinations were searched? (e.g., A*ABC and AA*BC
336 // are two distinct combinations)
337 int combinations_searched
; // FIX: could be > 2^31?