[Bio] / Babel / bin / source2ach.py Repository:
ViewVC logotype

Annotation of /Babel/bin/source2ach.py

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1.1 - (view) (download) (as text)

1 : tharriso 1.1 #!/usr/bin/env python
2 :    
3 :     import os, sys, re, hashlib
4 :     from optparse import OptionParser
5 :     from Bio import SeqIO
6 :     from Bio.Seq import Seq
7 :     from Bio.SeqRecord import SeqRecord
8 :     from Bio.Alphabet import IUPAC
9 :    
10 :    
11 :     __doc__ = """
12 :    
13 :     Script to take a repository file from the following formats:
14 :     genbank, kegg, swiss, fasta
15 :     and output the following tab seperated files:
16 :    
17 :     SOURCE.md52id2func:
18 :     md5, id, function, organism, source, beg_pos*, end_pos*, strand*, contig_id*, contig_desc*, contig_length*
19 :     SOURCE.md52seq:
20 :     md5, sequence
21 :     SOURCE.id2xref
22 :     id, external refrence 1, [external refrence 2, external refrence 3, ...]
23 :    
24 :     Values with * are optional, only genbank files contain that data."""
25 :    
26 :     desc_re = re.compile('RecName: Full=(.*?);')
27 :     entry_re = re.compile('^ENTRY\s+(\S+)\s+CDS\s+(.*)$', re.MULTILINE)
28 :     name_re = re.compile('^NAME\s+(.*)$', re.MULTILINE)
29 :     func_re = re.compile('^DEFINITION\s+(.*?)^\S', re.DOTALL | re.MULTILINE)
30 :     dbref_re = re.compile('^DBLINKS\s+(.*?)^\S', re.DOTALL | re.MULTILINE)
31 :     seq_re = re.compile('^AASEQ\s+\d+\s+(.*?)^\S', re.DOTALL | re.MULTILINE)
32 :     go_re = re.compile('^GO:(.*)')
33 :    
34 :     def get_org_map(filename):
35 :     omap = {}
36 :     if filename is None:
37 :     return omap
38 :     for text in kegg_iter(filename):
39 :     rec = get_kegg_rec(text)
40 :     if (not rec.name) or (not rec.description): continue
41 :     names = rec.name.split(', ')
42 :     if len(names) > 1:
43 :     omap[ names[1] ] = rec.description
44 :     return omap
45 :    
46 :     def kegg_iter(filename):
47 :     fhdl = file(filename, 'rb')
48 :     recs = ['']
49 :     for data in iter(lambda: fhdl.read(100000), ''):
50 :     recs = (recs[-1] + data).split("\n///\n")
51 :     for rec in recs[:-1]:
52 :     yield rec
53 :    
54 :     def get_kegg_rec(text):
55 :     record = SeqRecord(Seq("",IUPAC.protein), id='', name='', description='', dbxrefs=[], annotations={'organism':''})
56 :     entry_s = entry_re.search(text)
57 :     name_s = name_re.search(text)
58 :     func_s = func_re.search(text)
59 :     dbref_s = dbref_re.search(text)
60 :     seq_s = seq_re.search(text)
61 :     if entry_s:
62 :     record.id = entry_s.group(1)
63 :     record.annotations['organism'] = entry_s.group(2).strip()
64 :     if name_s:
65 :     record.name = name_s.group(1).strip()
66 :     if func_s:
67 :     record.description = re.sub('\s+', ' ', func_s.group(1).strip())
68 :     if dbref_s:
69 :     parts = dbref_s.group(1).split()
70 :     if (len(parts) % 2) == 0:
71 :     for i in range(0, len(parts)-1, 2):
72 :     record.dbxrefs.append(parts[i].replace('NCBI-', '') + parts[i+1])
73 :     if seq_s:
74 :     record.seq = Seq( re.sub('\s+', '', seq_s.group(1)), IUPAC.protein )
75 :     return record
76 :    
77 :     def format_factory(data_f, seq_f, ref_f, source, form, org_desc, get_go):
78 :     if form == 'genbank':
79 :     def parse_genbank(rec):
80 :     if (not rec) or (not rec.id) or (not rec.seq): return
81 :     cid = rec.id
82 :     cdesc = rec.description.rstrip('.')
83 :     clen = len(rec.seq)
84 :     org = rec.annotations['organism']
85 :     for feat in rec.features:
86 :     if feat.type == 'CDS':
87 :     if ( (not feat.qualifiers) or (not feat.location) or
88 :     ('translation' not in feat.qualifiers) or
89 :     ('protein_id' not in feat.qualifiers) or
90 :     ('product' not in feat.qualifiers) ): continue
91 :     seq = feat.qualifiers['translation'][0]
92 :     md5 = hashlib.md5(seq).hexdigest()
93 :     fid = feat.qualifiers['protein_id'][0]
94 :     func = feat.qualifiers['product'][0]
95 :     beg = feat.location.start.position
96 :     end = feat.location.end.position
97 :     if feat.strand:
98 :     strd = feat.strand
99 :     else:
100 :     strd = 1
101 :     if (beg > end):
102 :     strd = -1
103 :     beg = feat.location.end.position
104 :     end = feat.location.start.position
105 :     seq_f.write("%s\t%s\n" %(md5, seq))
106 :     data_f.write("\t".join([md5, fid, func, org, source, str(beg), str(end), str(strd), cid, cdesc, str(clen)]) + "\n")
107 :     if ('db_xref' in feat.qualifiers) and (len(feat.qualifiers['db_xref']) > 0):
108 :     ref_f.write(fid + "\t" + "\t".join(feat.qualifiers['db_xref']) + "\n")
109 :     return parse_genbank
110 :    
111 :     elif form == 'swiss':
112 :     def parse_swiss(rec):
113 :     if (not rec) or (not rec.id) or (not rec.seq): return
114 :     md5 = hashlib.md5(str(rec.seq)).hexdigest()
115 :     func = rec.description
116 :     srch = desc_re.search(func)
117 :     if srch: func = srch.group(1)
118 :     seq_f.write("%s\t%s\n" %(md5, str(rec.seq)))
119 :     data_f.write("\t".join([md5, rec.id, func, rec.annotations['organism'], source]) + "\n")
120 :     if len(rec.dbxrefs) > 0:
121 :     ref_str = rec.id
122 :     for ref in rec.dbxrefs:
123 :     go_m = go_re.match(ref)
124 :     if get_go and go_m:
125 :     # create seperate entry for GO ids
126 :     data_f.write("\t".join([md5, go_m.group(1), func, rec.annotations['organism'], "GO"]) + "\n")
127 :     ref_f.write(go_m.group(1) + "\t" + source + ":" + rec.id + "\n")
128 :     else:
129 :     ref_str += "\t" + ref
130 :     ref_f.write(ref_str + "\n")
131 :     return parse_swiss
132 :    
133 :     elif form == 'fasta':
134 :     def parse_fasta(rec):
135 :     if (not rec) or (not rec.id) or (not rec.seq): return
136 :     func = rec.description
137 :     md5 = hashlib.md5(str(rec.seq)).hexdigest()
138 :     seq_f.write("%s\t%s\n" %(md5, str(rec.seq)))
139 :     if func.startswith("("): func = rec.description.strip(')(')
140 :     if org_desc:
141 :     data_f.write("\t".join([md5, rec.id, '', func, source]) + "\n")
142 :     else:
143 :     data_f.write("\t".join([md5, rec.id, func, '', source]) + "\n")
144 :     return parse_fasta
145 :    
146 :     elif form == 'kegg':
147 :     def parse_kegg(text, orgmap):
148 :     rec = get_kegg_rec(text)
149 :     if (not rec.id) or (not rec.description) or (not rec.seq): return
150 :     md5 = hashlib.md5(str(rec.seq)).hexdigest()
151 :     org = rec.annotations['organism']
152 :     if org in orgmap:
153 :     org = orgmap[org]
154 :     seq_f.write("%s\t%s\n" %(md5, str(rec.seq)))
155 :     data_f.write("\t".join([md5, rec.id, rec.description, org, source]) + "\n")
156 :     if len(rec.dbxrefs) > 0:
157 :     ref_f.write(rec.id + "\t" + "\t".join(rec.dbxrefs) + "\n")
158 :     return parse_kegg
159 :    
160 :     else:
161 :     return None
162 :    
163 :    
164 :     usage = "usage: %prog [options] source input_file1 [input_file2 input_file3 ...]" + __doc__
165 :    
166 :     def main(args):
167 :     parser = OptionParser(usage=usage)
168 :     parser.add_option("-f", "--format", dest="format", metavar="FORMAT", default="fasta",
169 :     help="FORMAT inputed file is in. Must be one of following: genbank, kegg, swiss, fasta [default is fasta]")
170 :     parser.add_option("-m", "--org_map", dest="orgmap", metavar="FILE", default=None, help="Optional KEGG genome FILE for org name mapping")
171 :     parser.add_option("-d", "--out_dir", dest="outdir", metavar="DIR", default=None, help="DIR to write output [default is current dir]")
172 :     parser.add_option("-o", "--org_header", dest="orgheader", action="store_true", default=False, help="For fasta files, header description is organism [default is function]")
173 :     parser.add_option("-g", "--go_entry", dest="goentry", action="store_true", default=False, help="For swiss files, output for each GO alias in entry [default is off]")
174 :     parser.add_option("-v", "--verbose", dest="verbose", action="store_true", default=False, help="Wordy [default is off]")
175 :    
176 :     (opts, args) = parser.parse_args()
177 :     if len(args) < 2:
178 :     parser.error("Incorrect number of arguments")
179 :     source = args[0]
180 :     files = args[1:]
181 :     orgmap = get_org_map(opts.orgmap)
182 :     outdir = ''
183 :     if opts.outdir: outdir = opts.outdir
184 :    
185 :     data_file = open( os.path.join(outdir, source + ".md52id2func"), 'w' )
186 :     seq_file = open( os.path.join(outdir, source + ".md52seq"), 'w' )
187 :     ref_file = open( os.path.join(outdir, source + ".id2xref"), 'w' )
188 :    
189 :     parse_format = format_factory( data_file, seq_file, ref_file, source, opts.format, opts.orgheader, opts.goentry )
190 :     if not parse_format:
191 :     parser.error("Incorrect format entered")
192 :    
193 :     for f in files:
194 :     if opts.verbose:
195 :     sys.stdout.write("Parsing file %s ...\n" %f)
196 :     if (opts.format == 'kegg'):
197 :     for rec in kegg_iter(f):
198 :     parse_format(rec, orgmap)
199 :     else:
200 :     for rec in SeqIO.parse(f, opts.format):
201 :     parse_format(rec)
202 :    
203 :     data_file.close()
204 :     seq_file.close()
205 :     ref_file.close()
206 :     if opts.verbose:
207 :     sys.stdout.write("All Done\n")
208 :    
209 :     if __name__ == "__main__":
210 :     sys.exit(main(sys.argv))

MCS Webmaster
ViewVC Help
Powered by ViewVC 1.0.3