source: trunk/mzcms/parsers.py @ 25

Last change on this file since 25 was 25, checked in by j@…, 10 years ago

Pseudocode for parse_dats

File size: 10.0 KB
Line 
1"""Mascot dat format to InsPecT output converter.
2
3"""
4import sys
5import re
6import csv
7from optparse import OptionParser
8from string import Template
9from collections import defaultdict
10
11SEPARATOR = '--gc0p4Jq0M2Yt08jU534c0p'
12SECTION_TEMPLATE = Template(
13        r'Content-Type: application/x-Mascot; name="$section"')
14
15PEP_REGEX = re.compile(r'^q(\d+)'    # nativeid
16                        '_p\d+'      # rank
17                        '=\d+?'      # missed cleavages
18                        ',[\.\d]+?'  # peptide mass
19                        ',[-\.\d]+?' # delta mass
20                        ',\d+?'      # n ions matches
21                        ',(\w+?)'    # peptide string
22                        ',\d+?'      # peaks used for Ions1
23                        ',(\d+?)'    # modstring
24                        ',[\.\d]+?'  # score
25                        ',\d+?'      # ion series found
26                        ',\d+?'      # peaks used for Ions2
27                        ',\d+?'      # peaks used for Ions3
28                        ';(.+)$'     # protein accessions string
29                        )
30#fraction_regex = re.compile(r'fraction%3a%20(\d+)')
31RAWFN_REGEX = re.compile(r'raw%3a%20(.+raw)')
32
33class DatParser(object):
34    """Mascot Dat Parser
35    """
36    def __init__(self, decoy_str, scans_in_title=False):
37        self.decoy_regex = re.compile(decoy_str)
38        self.scans_in_title = scans_in_title
39        if scans_in_title:
40            self.scan_regex = re.compile(r'scan%3a%20(\d+)')
41
42    def _parse_spectra_fn(self, datfile):
43        """Parses the spectra file name used for Mascot search
44        """
45        for line in datfile:
46            if line.startswith("FILE="):
47                #XXX: There must be something in the standard library
48                # to normalize this properly
49                full_path = line.strip().split('=')[1]
50                norm_path = full_path.replace("\\", "/")
51                return norm_path.split("/")[-1]
52
53    def _parse_frag_mode(self, datfile):
54        """Get the fragmentation type by looking at dat parameters
55        """
56        # TODO: Handle CID and error for other instruments
57        for line in datfile:
58            if line.startswith("INSTRUMENT="):
59                if line.strip().split('=')[1] == "ETD-TRAP":
60                    return 'ETD'
61                else:
62                    return 'CID'
63
64    def _parse_annotations(self, datfile):
65        """Parses the peptide section of mascot
66        """
67        nativeid_annotations = defaultdict(list)
68        # Seek until Mascot peptide section
69        peptide_section = SECTION_TEMPLATE.substitute(section="peptides")
70        for line in datfile:
71            if line.strip() == peptide_section:
72                break
73        for line in datfile:
74            match = re.match(PEP_REGEX, line)
75            if match:
76                nativeid = int(match.group(1))
77                peptide_str = match.group(2)
78                mods_str = match.group(3)
79                accs = match.group(4)
80                peptide = apply_mods(peptide_str, mods_str)
81                is_decoy = check_decoy(accs, self.decoy_regex)
82                # Assuming that after a peptide line, always the next line is
83                # the line with flanking Aa
84                # q4_p2_terms=E,K:E,K:E,K:E,K:E,K
85                line = datfile.next().strip()
86                # (('E', 'K'), ('E', 'K'), ('E', 'K'))
87                flanking_pairs = set(tuple(x.split(','))
88                        for x in line.split("=")[1].split(":"))
89                for flanking_pair in flanking_pairs:
90                    annotation = '.'.join((
91                            flanking_pair[0], peptide, flanking_pair[1]))
92                    annotation = annotation.replace('-', '*')
93                    nativeid_annotations[nativeid].append(dict(
94                            Annotation=annotation, IsDecoy=is_decoy))
95            elif line.strip() == SEPARATOR:
96                return nativeid_annotations
97
98    def _parse_spectra(self, datfile):
99        """Gets the data coming from the input file from Mascot dat
100           format
101        """
102        nativeid_spectra = defaultdict(dict)
103        nativeid = 1
104        for line in datfile:
105            nativeid_str = ''.join(("query", str(nativeid)))
106            query_line = SECTION_TEMPLATE.substitute(section=nativeid_str)
107            if line.strip() == query_line:
108                current_id_mgf = nativeid_spectra[nativeid]
109                nativeid += 1
110            elif line.startswith("title="):
111                line = line.strip()
112                #TODO: handle fractions
113                #fraction = int(re.search(fraction_regex, line).group(1))
114                #current_id_mgf["Fraction"] = fraction
115                try:
116                    rawfn = re.search(RAWFN_REGEX, line).group(1)
117                except AttributeError:
118                    sys.exit("It seems there is no raw: field in TITLE")
119                rawfn = rawfn.replace('%2e', '.')
120                current_id_mgf["RawFile"] = rawfn
121                if self.scans_in_title:
122                    try:
123                        scan = int(re.search(self.scan_regex, line).group(1))
124                    except AttributeError:
125                        sys.exit("Default scan regex for TITLE field"
126                                 " doesn't match anything")
127                    current_id_mgf['Scan'] = scan
128            elif line.startswith("charge="):
129                charge = line.strip().split('=')[1]
130                current_id_mgf['Charge'] = charge
131            elif not self.scans_in_title and line.startswith("scans="):
132                scan = int(line.strip().split('=')[2])
133                current_id_mgf['Scan'] = scan
134        return nativeid_spectra
135
136    def parse(self, datfile):
137        """Takes a dat file and returns a dictionary of psms.
138        """
139        # TODO: handle frag_mode from mgf data
140        pkl_fn = self._parse_spectra_fn(datfile)
141        frag_mode = self._parse_frag_mode(datfile)
142        nativeid_annotations = self._parse_annotations(datfile)
143        nativeid_spectra = self._parse_spectra(datfile)
144        for native_id, annotations in nativeid_annotations.items():
145            for annotation in annotations:
146                annotation['#SpectraFile'] = pkl_fn
147                annotation['FragMode'] = frag_mode
148                annotation.update(nativeid_spectra[native_id])
149                yield annotation
150
151def apply_mods(peptide_str, mod_str):
152    """Takes a peptide string and applies mascot modification string. It
153       only supports C+57.
154    """
155    mod_str = mod_str[1:-1]
156    return ''.join((match != '0' and (p + '+57') or p for
157          p, match in zip(peptide_str, mod_str)))
158
159def check_decoy(accs, decoy_regex):
160    """Returns True if there is at least one target protein in Mascot
161       protein accession string.
162    """
163    prots = (y[1] for y in (x.split('"') for x in accs.split(',')))
164    for prot in prots:
165        match = re.search(decoy_regex, prot)
166        if not match:
167            return False
168    return True
169
170def write_csv(psms, out_file):
171    """Takes a dictionary of PSMs and writes a modified InsPecT output.
172    """
173    fieldnames = ("#SpectraFile",
174                  "RawFile",
175                  "Fraction",
176                  "Scan",
177                  "Annotation",
178                  "Charge",
179                  "FragMode",
180                  "IsDecoy")
181    header = dict()
182    for fieldname in fieldnames:
183        header[fieldname] = fieldname
184    writer = csv.DictWriter(out_file, fieldnames, dialect='excel-tab')
185    for psm in psms:
186        writer.writerow(psm)
187
188def parse_dats(dats_dir, proteins_factory=dict, peptides_factory=dict,
189               spectra_factory=dict, psms_factory=dict):
190    """Parses all the dat files in dats_dir
191    """
192    dats = get_dat_files(dats_dir)
193    proteins = proteins_factory()
194    peptides = peptides_factory()
195    spectra = spectra_factory()
196    psms = psms_factory()
197    dat_parser = DatParser(decoy_str, scans_in_title)
198    for dat in dats:
199        dat_proteins, dat_peptides, dat_spectra, dat_psms = dat_parser.parse(file)
200        proteins.update(dat_proteins)
201        peptides.update(dat_peptides)
202        spectra.update(dat_spectra)
203        psms.update(dat_psms)
204    return proteins, peptides, spectra, psms
205
206def main():
207    usage = "usage: %prog [options] mascotresult1.dat [mascotresult2.dat ...]"
208    optparser = OptionParser(usage)
209    optparser.add_option("-o",
210                      "--output",
211                      dest="out_fn",
212                      help="Output file [default: stdout]",
213                      metavar="FILE")
214    optparser.add_option("-s",
215                      "--scans-in-title",
216                      dest="scans_in_title",
217                      type="int",
218                      default="1",
219                      help="1 to parse the scan number from the TITLE field, "
220                           "0 to parse the scan number from the SCAN field "
221                           "[default: %default]")
222    optparser.add_option("-d",
223                      "--decoy-str",
224                      dest="decoy_str",
225                      default="^IPI:REV_:IPI",
226                      help="1 to parse the scan number from the TITLE field, "
227                           "0 to parse the scan number from the SCAN field "
228                           "[default: %default]")
229    (options, args) = optparser.parse_args()
230
231    if len(args) == 0:
232        optparser.print_usage()
233
234    else:
235        for arg in args:
236            if not arg.endswith(".dat"):
237                optparser.error("Provide only mascot dat files")
238            else:
239                if options.out_fn:
240                    out_fn = options.out_file
241                    out_file = open(out_fn, 'wb')
242                else:
243                    out_file = sys.stdout
244                for arg in args:
245                    with open(arg) as in_file:
246                        dat_parser = DatParser(options.decoy_str,
247                                scans_in_title = options.scans_in_title)
248                        psms = dat_parser.parse(in_file)
249                        write_csv(psms, out_file)
250
251if __name__ == '__main__':
252    main()
Note: See TracBrowser for help on using the repository browser.