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.
35 class mass_regime_parameters
{
41 // residue char (incl '['/']') -> mass (per regime) with any fixed mod
42 std::vector
<double> fixed_residue_mass
;
45 // FIX: want these class members to all be static, but had SWIG trouble. The
46 // member 'the' is a handle to the singleton instance, for now.
49 bool estimate_only
; // just estimate work required
50 bool show_progress
; // running progress on stderr
52 std::vector
<double> ln_factorial
; // factorial[n] == ln(n!)
54 // deuterium not (yet) implemented
57 std::vector
<mass_regime_parameters
> parent_mass_regime
;
58 std::vector
<mass_regime_parameters
> fragment_mass_regime
;
61 double parent_mass_tolerance_1
; // for +1
62 double parent_mass_tolerance_max
; // for +N (typically +3)
64 int minimum_peptide_length
;
66 double fragment_mass_tolerance
;
67 int intensity_class_count
;
69 static parameters the
;
75 int sequence_index
; // index of this sequence in the database FIX:nn
76 int sequence_offset
; // offset of this run's start in sequence
78 std::vector
<int> cleavage_points
;
79 std::string name
; // e.g., initial defline word
82 sequence_run(const int sequence_index
, const int sequence_offset
,
83 const std::string sequence
, std::vector
<int> cleavage_points
,
84 const std::string name
)
85 : sequence_index(sequence_index
), sequence_offset(sequence_offset
),
86 sequence(sequence
), cleavage_points(cleavage_points
),
89 // generate non-specific cleavage points
90 assert(sequence
.size() < INT_MAX
);
91 if (this->cleavage_points
.empty())
92 for (int i
=0; i
<=static_cast<int>(sequence
.size()); i
++)
93 this->cleavage_points
.push_back(i
);
98 // information about the search that's passed in to--but not changed by--the
100 class search_context
{
103 int mass_regime_index
;
104 std::string pca_residues
;
107 double N_delta
; // 0 if none
108 double C_delta
; // 0 if none
110 // maps residue char to list of indices for _count/_delta
111 std::vector
< std::vector
<int> > delta_bag_lookup
;
112 std::vector
<double> delta_bag_delta
;
113 std::vector
<int> delta_bag_count
;
115 double parent_fixed_mass
;
116 double fragment_N_fixed_mass
;
117 double fragment_C_fixed_mass
;
119 // general search parameters
120 int maximum_missed_cleavage_sites
;
122 std::vector
<sequence_run
> sequence_runs
;
126 enum ion
{ ION_MIN
, ION_Y
=ION_MIN
, ION_B
, ION_MAX
}; // NB: [ION_MIN, ION_MAX)
128 inline void operator++(ion
&i
) { i
= ion(i
+ 1); }
129 inline void operator++(ion
&i
, int) { i
= ion(i
+ 1); }
136 int intensity_class
; // 0..N; -1 == unclassified
138 explicit peak(double mz
=0, double intensity
=0, int intensity_class
=-1)
139 : mz(mz
), intensity(intensity
), intensity_class(intensity_class
) {
142 char *__repr__() const;
144 static bool less_mz(peak x
, peak y
) { return x
.mz
< y
.mz
; }
145 // FIX: delete one of these
146 static bool less_intensity(peak x
, peak y
) {
147 return x
.intensity
< y
.intensity
;
149 static bool greater_intensity(peak x
, peak y
) {
150 return x
.intensity
> y
.intensity
;
153 static double get_mz(double mass
, int charge
) {
154 const parameters
&CP
= parameters::the
;
155 return mass
/charge
+ CP
.proton_mass
;
160 // A spectrum searched on multiple charges will be turned into separate
161 // spectra having differing id's, but the same physical id.
165 double mass
; // [M + H+], aka "neutral mass" (protonated)
167 static const int MAX_SUPPORTED_CHARGE
= 10;
168 std::vector
<peak
> peaks
; // always ordered by increasing m/z!
170 std::vector
<int> intensity_class_counts
; // number of peaks in each class
172 // This is the mz range of peaks, but with a buffer on each end of size
173 // fragment_mass_tolerance.
177 double total_ion_current
;
179 // This is the total number of bins, 2*fragment_tolerance mz wide in peak
180 // range. So, for example, if the fragment tolerance is 1 and the min and
181 // max mz are 200 and 1800, total peak bins would be 800.
183 int empty_peak_bins
; // number of bins that aren't occupied by a peak
186 int file_id
; // file # of spectra's source
187 int id
; // unique for all spectra
188 int physical_id
; // FIX: unneeded?
190 unsigned long comparisons
; // times scored against theoretical spectra
192 // Construct an empty spectrum.
193 explicit spectrum(double mass
=0,
194 int charge
=0) : mass(mass
), charge(charge
), min_peak_mz(0),
195 max_peak_mz(0), total_ion_current(0),
196 total_peak_bins(0), empty_peak_bins(0),
197 file_id(-1), physical_id(-1),
200 if (charge
> MAX_SUPPORTED_CHARGE
)
201 throw std::invalid_argument("attempt to create a spectrum with greater"
202 " than supported charge");
206 char *__repr__() const;
209 void set_peaks_from_matrix(const std::vector
< std::vector
<double> > &m
) {
210 peaks
.resize(m
.size());
211 std::vector
<peak
>::iterator p_it
= peaks
.begin();
212 for (std::vector
< std::vector
<double> >::const_iterator it
= m
.begin();
213 it
!= m
.end(); it
++, p_it
++) {
215 throw std::invalid_argument("invalid matrix (must be size N x 2)");
217 p_it
->intensity
= (*it
)[1];
222 // Read spectra from file in ms2 format, tagging them with file_id. Before
223 // reading, seek to absolute position offset_begin. If offset_end != -1,
224 // any spectra that begin after position offset_end in the file are not
226 static std::vector
<spectrum
>
227 read_spectra_from_ms2(FILE *f
, const int file_id
, const long offset_begin
,
228 const long offset_end
);
231 // Filter peaks to limit their number according to TIC_cutoff_proportion,
232 // and to remove those too large to be fragment products. (Peaks are
233 // assumed to be ordered by increasing mz.)
234 void filter_peaks(double TIC_cutoff_proportion
,
235 double parent_mass_tolerance_max
);
237 // Classify peaks and update class_counts.
238 void classify(int intensity_class_count
, double intensity_class_ratio
,
239 double fragment_mass_tolerance
);
241 // Store the list of spectra that search_peptide will search against, and
242 // also build spectrum_mass_index.
243 static void set_searchable_spectra(const std::vector
<spectrum
> &spectra
);
245 // Search sequence runs for matches according to the context against the
246 // spectra. Updates score_stats and the number of candidate spectra found.
247 static void search_runs(const search_context
&context
, score_stats
&stats
);
249 // conceptually these are 'protected:', but we're taking it easy here
251 void set_id() { id
= next_id
++; assert(next_id
> 0); }
254 static int next_physical_id
;
256 static std::vector
<spectrum
> searchable_spectra
;
257 // parent mass -> searchable_spectra index
258 static std::multimap
<double,
259 std::vector
<spectrum
>::size_type
> spectrum_mass_index
;
263 struct mass_trace_item
{
270 // This is everything we want to remember about a match, so that we can report
272 // FIX: Are any of these fields unneeded?
276 int missed_cleavage_count
; // FIX: unneeded
278 int peptide_begin
; // not run-relative
279 std::string peptide_sequence
;
280 double predicted_parent_mass
;
281 std::string sequence_name
;
282 std::vector
<mass_trace_item
> mass_trace
;
284 match() : score(0), missed_cleavage_count(-1), spectrum_index(-1),
285 peptide_begin(-1), predicted_parent_mass(0) {
290 // This holds information about the best peptide/spectrum matches found, and
291 // some statistics. Essentially, it holds the results of the search process.
294 score_stats(int spectrum_count
, int best_result_count
)
295 : candidate_spectrum_count(0), spectra_with_candidates(0),
296 combinations_searched(0) {
297 best_matches
.resize(spectrum_count
);
298 assert(best_result_count
>= 1);
299 for (int i
=0; i
<spectrum_count
; i
++)
300 best_matches
[i
].resize(best_result_count
);
303 // spectrum index -> list of match (only top N matches per spectrum) The
304 // best match lists are of fixed length (typically 5), and ordered
305 // best-first. Trailing entries with score 0 are null matches (to simplify
307 std::vector
< std::vector
<match
> > best_matches
;
310 // How many times was a spectrum scored against a peptide?
311 unsigned long long candidate_spectrum_count
; // may be > 2^32
312 // How many spectra have at a non-zero best score?
313 int spectra_with_candidates
;
314 // How many different mod combinations were searched? (e.g., A*ABC and AA*BC
315 // are two distinct combinations)
316 int combinations_searched
; // FIX: could be > 2^31?