source: trunk/mzcms/parsers.py @ 27

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

specified default mgf regex to parse

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