1 #!/usr/bin/env greylag-python
3 '''Generate a set of sqt files from the specified greylag search result
4 file. The sqt files will be created in the current directory, with names
5 corresponding to the searched ms2 files.
8 from __future__
import with_statement
11 greylag, Copyright (C) 2006-2007, Stowers Institute for Medical Research
13 This program is free software; you can redistribute it and/or modify
14 it under the terms of the GNU General Public License as published by
15 the Free Software Foundation; either version 2 of the License, or
16 (at your option) any later version.
18 This program is distributed in the hope that it will be useful,
19 but WITHOUT ANY WARRANTY; without even the implied warranty of
20 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
21 GNU General Public License for more details.
23 You should have received a copy of the GNU General Public License along
24 with this program; if not, write to the Free Software Foundation, Inc.,
25 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
37 from pprint
import pprint
42 print >> sys
.stderr
, 'warning:', s
44 sys
.exit('error: ' + s
)
47 def print_header(f
, r
):
48 print >> f
, "H\tSQTGenerator\tgreylag"
49 print >> f
, "H\tSQTGeneratorVersion\t%s" % __version__
50 for db
in r
['databases']:
51 print >> f
, "H\tDatabase\t%s" % db
52 for pk
, pv
in sorted(r
['parameters'].items()):
53 print >> f
, "H\tParameter\t%s\t%s" % (pk
, pv
)
54 # How should these be specified?
55 # What does DTASelect do with these?
56 # H StaticMod C=160.1388
60 def print_spectrum(f
, sp_matches
):
61 """Print the lines associated with the given spectrum."""
62 spectrum
, sp_best_matches
= sp_matches
64 scan_low
, scan_high
, _rest
= spectrum
['name'].split('.', 2)
65 print >> f
, '\t'.join(str(v
) for v
in
66 ["S", scan_low
, scan_high
, spectrum
['charge'], 0,
67 'honk', spectrum
['mass'],
68 round(math
.log(spectrum
['total_ion_current']), 4),
69 0, spectrum
['comparisons']])
71 best_scores
= sorted(list(set(m
['score'] for m
in sp_best_matches
)))
73 rank_map
= dict(zip(best_scores
, range(1,len(best_scores
)+1)))
75 best_score
= best_scores
[0]
76 for match
in sp_best_matches
:
77 score
= match
['score']
80 rank
= rank_map
[score
]
81 score_delta
= (best_score
- score
) / best_score
83 # Note: scores, being log probability, are non-positive, but we
84 # flip the sign in the SQT output
85 print >> f
, '\t'.join(str(v
) for v
in
87 round(match
['predicted_parent_mass'], 5),
88 round(score_delta
, 4), round(-score
, 4), 0,
89 # 1 of 2 ions found--keep DTASelect happy
91 "-.%s.-" % match
['peptide_sequence'], 'U'])
93 assert len(match
['sequence_name']) == len(match
['peptide_begin'])
94 for sn
, pb
in zip(match
['sequence_name'], match
['peptide_begin']):
95 print >> f
, 'L\t%s\t%s' % (sn
, pb
)
98 def main(args
=sys
.argv
[1:]):
99 parser
= optparse
.OptionParser(usage
=
100 "usage: %prog [options] <result-file>",
101 description
=__doc__
, version
=__version__
)
102 pa
= parser
.add_option
103 pa("-d", "--output-directory", dest
="output_directory",
104 help="directory where output files are written [default is '.']",
106 pa("-v", "--verbose", action
="store_true", dest
="verbose",
108 pa("--dump", action
="store_true", dest
="dump",
109 help="dump the result file (for debugging)")
110 pa("--copyright", action
="store_true", dest
="copyright",
111 help="print copyright and exit")
112 (options
, args
) = parser
.parse_args(args
=args
)
114 if options
.copyright
:
122 with contextlib
.closing(gzip
.open(args
[0])) as r_file
:
123 r
= cPickle
.load(r_file
)
129 spectrum_fns
= r
['spectrum files']
130 assert len(spectrum_fns
) == len(set(spectrum_fns
)) # check uniqueness
132 matches
= r
['matches'].items()
134 # order matches by (spectrum filename, spectrum name)
135 def filename_specname_order(m
):
138 return sp
['file_id'], sp
['name']
140 matches
.sort(key
=filename_specname_order
)
142 for spectrum_n
, spectrum_fn
in enumerate(spectrum_fns
):
143 assert os
.path
.dirname(spectrum_fn
) == ''
144 sqt_fn
= os
.path
.splitext(spectrum_fn
)[0] + '.sqt'
145 with
open(sqt_fn
, 'w') as sqtf
:
146 print_header(sqtf
, r
)
147 for match
in matches
:
148 if match
[1][0]['file_id'] == spectrum_n
:
149 print_spectrum(sqtf
, match
[1])
152 if __name__
== '__main__':