source: trunk/mzcms/parsers.py @ 56

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

Added update method for proteins folder

File size: 13.3 KB
Line 
1"""MS/MS search engine results parsers.
2
3"""
4import sys
5import os
6import re
7import csv
8import urllib
9import fnmatch
10from optparse import OptionParser
11from string import Template
12from collections import defaultdict
13from collections import namedtuple
14from itertools import izip
15
16SEPARATOR = '--gc0p4Jq0M2Yt08jU534c0p'
17SECTION_TEMPLATE = Template(
18        r'Content-Type: application/x-Mascot; name="$section"')
19
20PEP_REGEX = re.compile(r'^q(\d+)'      # nativeid
21                        '_p(\d+)'      # rank
22                        '=\d+?'        # missed cleavages
23                        ',([\.\d]+?)'  # peptide mass
24                        ',([-\.\d]+?)' # delta mass
25                        ',\d+?'        # n ions matches
26                        ',(\w+?)'      # peptide string
27                        ',\d+?'        # peaks used for Ions1
28                        ',(\d+?)'      # modstring
29                        ',([\.\d]+?)'  # score
30                        ',\d+?'        # ion series found
31                        ',\d+?'        # peaks used for Ions2
32                        ',\d+?'        # peaks used for Ions3
33                        ';(.+)$'       # protein accessions string
34                        )
35#fraction_regex = re.compile(r'fraction%3a%20(\d+)')
36
37MascotPsm = namedtuple('MascotPsm',
38                       'rank '
39                       'pep_mass '
40                       'delta_mass '
41                       'pep_seq '
42                       'score '
43                       'prot_ids'
44                       )
45
46class DatParser(object):
47    """Mascot Dat Parser
48    """
49    def __init__(self,
50                 protein_factory,
51                 peptide_factory,
52                 spectrum_factory,
53                 psm_factory,
54                 scan_str=r'FinneganScanNumber%3a%20(\d+)%20',
55                 rawfn_str=r'RawFile%3a%20(.+raw)',
56                 decoy_str=r'^REV_IPI',
57                 contaminant_str=r'^CON_IPI',
58                 # XXX: Fix default factories
59                 ):
60        self.scan_regex = re.compile(scan_str)
61        decoy_regex = re.compile(decoy_str)
62        contaminant_regex = re.compile(contaminant_str)
63        # Contaminants are not exactly decoy psms
64        self.non_target_regexes = (contaminant_regex, decoy_regex)
65        self.rawfn_regex = re.compile(rawfn_str)
66        self.protein_factory = protein_factory
67        self.peptide_factory = peptide_factory
68        self.spectrum_factory=spectrum_factory
69        self.psm_factory=psm_factory
70
71    def _parse_spectra_fn(self, dat_file):
72        """Parses the spectra file name used for Mascot search
73        """
74        for line in dat_file:
75            if line.startswith("FILE="):
76                full_path = line.strip().split('=')[1]
77                norm_path = full_path.replace('\\', '/')
78                return norm_path.split("/")[-1]
79
80    def _parse_frag_mode(self, dat_file):
81        """Get the fragmentation type by looking at dat parameters
82        """
83        # TODO: Handle CID and error for other instruments
84        for line in dat_file:
85            if line.startswith("INSTRUMENT="):
86                frag_line = line.strip().split('=')[1]
87                if frag_line == "ETD-TRAP":
88                    return 'ETD'
89                elif frag_line == "ESI-TRAP":
90                    return 'CID'
91                else:
92                    return frag_line
93
94    def _parse_summary(self, dat_file):
95        """Parse summary section
96        """
97        summspec = list()
98        summary_section = SECTION_TEMPLATE.substitute(section="summary")
99        for line in dat_file:
100            if line.strip() == summary_section:
101                break
102        for line in dat_file:
103            if line.startswith('qexp'):
104                r_mz, r_charge = line.strip().split('=')[1].split(',')
105                mz = float(r_mz)
106                charge = int(r_charge.rstrip('+'))
107                summspec.append((mz, charge))
108            elif line.strip() == SEPARATOR:
109                return summspec
110
111    def _parse_psms(self, dat_file):
112        """Parses the peptide section of mascot
113        """
114        nativeid_psms = defaultdict(list)
115        # Seek until Mascot peptide section
116        peptide_section = SECTION_TEMPLATE.substitute(section="peptides")
117        for line in dat_file:
118            if line.strip() == peptide_section:
119                break
120        for line in dat_file:
121            match = re.match(PEP_REGEX, line)
122            if match:
123                nativeid = int(match.group(1))
124                rank = int(match.group(2))
125                pep_mass = float(match.group(3))
126                delta_mass = float(match.group(4))
127                pep_seq = match.group(5)
128                mods_str = match.group(6)
129                score = float(match.group(7))
130                accs = match.group(8)
131                prot_ids = get_prot_ids(accs)
132                is_target = check_target(prot_ids,
133                        self.non_target_regexes)
134                if is_target and (rank == 1 or rank == 2):
135                    mascot_psm = MascotPsm(
136                            rank=rank,
137                            pep_mass=pep_mass,
138                            delta_mass=delta_mass,
139                            pep_seq=pep_seq,
140                            score=score,
141                            prot_ids=prot_ids
142                            )
143                    nativeid_psms[nativeid].append(mascot_psm)
144            elif line.strip() == SEPARATOR:
145                return nativeid_psms
146
147    def _parse_extraspec(self, dat_file):
148        """Gets the data coming from the input file from Mascot dat
149           format
150        """
151        extraspec = list()
152        for line in dat_file:
153            if line.startswith("title="):
154                line = line.strip()
155                try:
156                    quoted_rawfn = re.search(self.rawfn_regex, line).group(1)
157                    rawfn = urllib.unquote(quoted_rawfn)
158                except AttributeError:
159                    sys.exit("It seems the raw regex is not valid")
160                try:
161                    scan = int(re.search(self.scan_regex, line).group(1))
162                except AttributeError:
163                    sys.exit("Default scan regex for TITLE field"
164                             " doesn't match anything")
165                extraspec.append((rawfn, scan))
166        return extraspec
167
168    def parse(self, dat_file):
169        """Takes a dat file and returns a dictionary of psms.
170        """
171        pkl_fn = self._parse_spectra_fn(dat_file)
172        frag_mode = self._parse_frag_mode(dat_file)
173        summspec = self._parse_summary(dat_file)
174        nativeid_psms = self._parse_psms(dat_file)
175        extraspec = self._parse_extraspec(dat_file)
176        proteins = dict()
177        spectra = dict()
178        peptides = dict()
179        psms = dict()
180        for summ, extra in izip(summspec, extraspec):
181            spectrum = self.spectrum_factory(
182                    prec_mz=summ[0],
183                    prec_charge=summ[1],
184                    run=extra[0],
185                    scan=extra[1],
186                    peaks="TBI"
187                    )
188            spectra[(spectrum.scan, spectrum.run)] = spectrum
189        for nativeid, m_psms in nativeid_psms.iteritems():
190            for m_psm in m_psms:
191                peptide = self.peptide_factory(
192                        sequence=m_psm.pep_seq,
193                        mass=m_psm.pep_mass
194                        )
195                if peptide.sequence not in peptides:
196                    peptides[peptide.sequence] = peptide
197                run, spec = extraspec[nativeid - 1]
198                spec_id = spec, run
199                psm = self.psm_factory(
200                        score=m_psm.score,
201                        rank=m_psm.rank,
202                        delta_mass=m_psm.delta_mass,
203                        pep_ref=peptide,
204                        spec_ref=spectra[spec_id]
205                        )
206                psms[(peptide.sequence, spec_id)] = psm
207                for prot_id in m_psm.prot_ids:
208                    if prot_id not in proteins:
209                        protein = self.protein_factory(
210                                native_id=prot_id,
211                                sequence="TBI",
212                                pep_refs=[peptide]
213                                )
214                        proteins[prot_id] = protein
215                    else:
216                        if peptide not in proteins[prot_id].pep_refs:
217                            proteins[prot_id].pep_refs.append(peptide)
218        return proteins, peptides, spectra, psms
219
220def get_prot_ids(accs):
221    return tuple(x[1].strip('"')
222                 for x in (y.split(':')
223                 for y in accs.split(',"')
224                           )
225                 )
226
227def check_target(prot_ids, regexes):
228    """Returns True if there is at least one target protein in Mascot
229       protein accession string.
230    """
231    for prot_id in prot_ids:
232        for regex in regexes:
233            if re.match(regex, prot_id):
234                return False
235        return True
236
237def apply_mods(peptide_str, mod_str):
238    """Takes a peptide string and applies mascot modification string. It
239       only supports C+57.
240    """
241    mod_str = mod_str[1:-1]
242    return ''.join((match != '0' and (p + '+57') or p for
243          p, match in zip(peptide_str, mod_str)))
244
245def write_csv(psms, out_file):
246    """Takes a dictionary of PSMs and writes a modified InsPecT output.
247    """
248    fieldnames = ("#SpectraFile",
249                  "RawFile",
250                  "Fraction",
251                  "Scan",
252                  "Annotation",
253                  "Charge",
254                  "FragMode",
255                  "IsDecoy")
256    header = dict()
257    for fieldname in fieldnames:
258        header[fieldname] = fieldname
259    writer = csv.DictWriter(out_file, fieldnames, dialect='excel-tab')
260    for psm in psms:
261        writer.writerow(psm)
262
263# XXX: Use better defaults for containers and factories
264def parse_dats(dats_dir, proteins_container=dict, peptides_container=dict,
265               spectra_container=dict, psms_container=dict,
266               protein_factory=dict, peptide_factory=dict,
267               spectrum_factory=dict, psm_factory=dict):
268    """Parses all the dat files in dats_dir
269    """
270    proteins = proteins_container
271    peptides = peptides_container
272    spectra = spectra_container
273    psms = psms_container
274    dat_parser = DatParser(protein_factory=protein_factory,
275                           peptide_factory=peptide_factory,
276                           spectrum_factory=spectrum_factory,
277                           psm_factory=psm_factory)
278    import transaction
279    for root, dirs, files in os.walk(dats_dir):
280        for dat_fn in files:
281            if fnmatch.fnmatch(dat_fn, '*.dat'):
282                with open(os.path.join(root, dat_fn)) as dat_file:
283                    ser_proteins, ser_peptides, \
284                    ser_spectra, ser_psms = dat_parser.parse(dat_file)
285                    proteins.update(ser_proteins)
286                    transaction.savepoint()
287                    for pep_id, pep in ser_peptides.items():
288                        if not pep_id in peptides:
289                            peptides[str(pep_id)] = pep
290                    transaction.savepoint()
291                    for spec_id, spec in ser_spectra.items():
292                        if spec_id not in spectra:
293                            spectra[str(spec_id)] = spec
294                    transaction.savepoint()
295                    for psm_id, psm in ser_psms.items():
296                        if psm not in psms:
297                            psms[str(psm_id)] = psm
298                    transaction.savepoint()
299
300def main():
301    usage = "usage: %prog [options] mascotresult1.dat [mascotresult2.dat ...]"
302    optparser = OptionParser(usage)
303    optparser.add_option("-o",
304                      "--output",
305                      dest="out_fn",
306                      help="Output file [default: stdout]",
307                      metavar="FILE")
308    optparser.add_option("-s",
309                      "--scans-in-title",
310                      dest="scans_in_title",
311                      type="int",
312                      default="1",
313                      help="1 to parse the scan number from the TITLE field, "
314                           "0 to parse the scan number from the SCAN field "
315                           "[default: %default]")
316    optparser.add_option("-d",
317                      "--decoy-str",
318                      dest="decoy_str",
319                      default=r'^IPI:REV_:IPI',
320                      help="1 to parse the scan number from the TITLE field, "
321                           "0 to parse the scan number from the SCAN field "
322                           "[default: %default]")
323    (options, args) = optparser.parse_args()
324
325    if len(args) == 0:
326        optparser.print_usage()
327
328    else:
329        for arg in args:
330            if not arg.endswith(".dat"):
331                optparser.error("Provide only mascot dat files")
332            else:
333                if options.out_fn:
334                    out_fn = options.out_file
335                    out_file = open(out_fn, 'wb')
336                else:
337                    out_file = sys.stdout
338                for arg in args:
339                    with open(arg) as in_file:
340                        dat_parser = DatParser(options.decoy_str,
341                                scans_in_title = options.scans_in_title)
342                        psms = dat_parser.parse(in_file)
343                        write_csv(psms, out_file)
344
345if __name__ == '__main__':
346    main()
Note: See TracBrowser for help on using the repository browser.