source: trunk/mzcms/parsers.py @ 35

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

Renamed annotations to psms in parsers

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