source: trunk/mzcms/parsers.py @ 55

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

Better variable naming to update data folders

File size: 13.4 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                    for prot_id, prot in ser_proteins.items():
286                        if prot_id not in proteins:
287                            proteins[str(prot_id)] = prot
288                    transaction.savepoint()
289                    for pep_id, pep in ser_peptides.items():
290                        if not pep_id in peptides:
291                            peptides[str(pep_id)] = pep
292                    transaction.savepoint()
293                    for spec_id, spec in ser_spectra.items():
294                        if spec_id not in spectra:
295                            spectra[str(spec_id)] = spec
296                    transaction.savepoint()
297                    for psm_id, psm in ser_psms.items():
298                        if psm not in psms:
299                            psms[str(psm_id)] = psm
300                    transaction.savepoint()
301
302def main():
303    usage = "usage: %prog [options] mascotresult1.dat [mascotresult2.dat ...]"
304    optparser = OptionParser(usage)
305    optparser.add_option("-o",
306                      "--output",
307                      dest="out_fn",
308                      help="Output file [default: stdout]",
309                      metavar="FILE")
310    optparser.add_option("-s",
311                      "--scans-in-title",
312                      dest="scans_in_title",
313                      type="int",
314                      default="1",
315                      help="1 to parse the scan number from the TITLE field, "
316                           "0 to parse the scan number from the SCAN field "
317                           "[default: %default]")
318    optparser.add_option("-d",
319                      "--decoy-str",
320                      dest="decoy_str",
321                      default=r'^IPI:REV_:IPI',
322                      help="1 to parse the scan number from the TITLE field, "
323                           "0 to parse the scan number from the SCAN field "
324                           "[default: %default]")
325    (options, args) = optparser.parse_args()
326
327    if len(args) == 0:
328        optparser.print_usage()
329
330    else:
331        for arg in args:
332            if not arg.endswith(".dat"):
333                optparser.error("Provide only mascot dat files")
334            else:
335                if options.out_fn:
336                    out_fn = options.out_file
337                    out_file = open(out_fn, 'wb')
338                else:
339                    out_file = sys.stdout
340                for arg in args:
341                    with open(arg) as in_file:
342                        dat_parser = DatParser(options.decoy_str,
343                                scans_in_title = options.scans_in_title)
344                        psms = dat_parser.parse(in_file)
345                        write_csv(psms, out_file)
346
347if __name__ == '__main__':
348    main()
Note: See TracBrowser for help on using the repository browser.