Mercurial > repos > vipints > fml_gff3togtf

comparison GFFtools-GX/helper.py @ 3:ff2c2e6f4ab3

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Uploaded version 2.0.0 of gfftools ready to import to local instance
author vipints
date Wed, 11 Jun 2014 16:29:25 -0400
parents
children
comparison
equal deleted inserted replaced
2:db3c67b03d55 3:ff2c2e6f4ab3
1 #!/usr/bin/env python
2 """
3 Common utility functions
4 """
5
6 import os
7 import re
8 import sys
9 import gzip
10 import bz2
11 import numpy
12
13 def init_gene():
14 """
15 Initializing the gene structure
16 """
17
18 gene_det = [('id', 'f8'),
19 ('anno_id', numpy.dtype),
20 ('confgenes_id', numpy.dtype),
21 ('name', 'S25'),
22 ('source', 'S25'),
23 ('gene_info', numpy.dtype),
24 ('alias', 'S15'),
25 ('name2', numpy.dtype),
26 ('strand', 'S2'),
27 ('score', 'S15'),
28 ('chr', 'S15'),
29 ('chr_num', numpy.dtype),
30 ('paralogs', numpy.dtype),
31 ('start', 'f8'),
32 ('stop', 'f8'),
33 ('transcripts', numpy.dtype),
34 ('transcript_type', numpy.dtype),
35 ('transcript_info', numpy.dtype),
36 ('transcript_status', numpy.dtype),
37 ('transcript_valid', numpy.dtype),
38 ('exons', numpy.dtype),
39 ('exons_confirmed', numpy.dtype),
40 ('cds_exons', numpy.dtype),
41 ('utr5_exons', numpy.dtype),
42 ('utr3_exons', numpy.dtype),
43 ('tis', numpy.dtype),
44 ('tis_conf', numpy.dtype),
45 ('tis_info', numpy.dtype),
46 ('cdsStop', numpy.dtype),
47 ('cdsStop_conf', numpy.dtype),
48 ('cdsStop_info', numpy.dtype),
49 ('tss', numpy.dtype),
50 ('tss_info', numpy.dtype),
51 ('tss_conf', numpy.dtype),
52 ('cleave', numpy.dtype),
53 ('cleave_info', numpy.dtype),
54 ('cleave_conf', numpy.dtype),
55 ('polya', numpy.dtype),
56 ('polya_info', numpy.dtype),
57 ('polya_conf', numpy.dtype),
58 ('is_alt', 'f8'),
59 ('is_alt_spliced', 'f8'),
60 ('is_valid', numpy.dtype),
61 ('transcript_complete', numpy.dtype),
62 ('is_complete', numpy.dtype),
63 ('is_correctly_gff3_referenced', 'S5'),
64 ('splicegraph', numpy.dtype) ]
65
66 return gene_det
67
68 def open_file(fname):
69 """
70 Open the file (supports .gz .bz2) and returns the handler
71
72 @args fname: input file name for reading
73 @type fname: str
74 """
75
76 try:
77 if os.path.splitext(fname)[1] == ".gz":
78 FH = gzip.open(fname, 'rb')
79 elif os.path.splitext(fname)[1] == ".bz2":
80 FH = bz2.BZ2File(fname, 'rb')
81 else:
82 FH = open(fname, 'rU')
83 except Exception as error:
84 sys.exit(error)
85
86 return FH
87
88 def add_CDS_phase(strand, cds):
89 """
90 Calculate CDS phase and add to the CDS exons
91
92 @args strand: feature strand information
93 @type strand: +/-
94 @args cds: coding exon coordinates
95 @type cds: numpy array [[int, int, int]]
96 """
97
98 cds_region, cds_flag = [], 0
99 if strand == '+':
100 for cdspos in cds:
101 if cds_flag == 0:
102 cdspos = (cdspos[0], cdspos[1], 0)
103 diff = (cdspos[1]-(cdspos[0]-1))%3
104 else:
105 xy = 0
106 if diff == 0:
107 cdspos = (cdspos[0], cdspos[1], 0)
108 elif diff == 1:
109 cdspos = (cdspos[0], cdspos[1], 2)
110 xy = 2
111 elif diff == 2:
112 cdspos = (cdspos[0], cdspos[1], 1)
113 xy = 1
114 diff = ((cdspos[1]-(cdspos[0]-1))-xy)%3
115 cds_region.append(cdspos)
116 cds_flag = 1
117 elif strand == '-':
118 cds.reverse()
119 for cdspos in cds:
120 if cds_flag == 0:
121 cdspos = (cdspos[0], cdspos[1], 0)
122 diff = (cdspos[1]-(cdspos[0]-1))%3
123 else:
124 xy = 0
125 if diff == 0:
126 cdspos = (cdspos[0], cdspos[1], 0)
127 elif diff == 1:
128 cdspos = (cdspos[0], cdspos[1], 2)
129 xy = 2
130 elif diff == 2:
131 cdspos = (cdspos[0], cdspos[1], 1)
132 xy = 1
133 diff = ((cdspos[1]-(cdspos[0]-1))-xy)%3
134 cds_region.append(cdspos)
135 cds_flag = 1
136 cds_region.reverse()
137 return cds_region
138
139 def buildUTR(cc, ec, strand):
140 """
141 Build UTR regions from a given set of CDS and exon coordiantes of a gene
142
143 @args cc: coding exon coordinates
144 @type cc: numpy array [[int, int, int]]
145 @args ec: exon coordinates
146 @type ec: numpy array [[int, int]]
147 @args strand: feature strand information
148 @type strand: +/-
149 """
150
151 utr5 = []
152 utr3 = []
153 if strand == '+':
154 cds_s = cc[0][0]
155 for ex in ec:
156 if ex[0] <= cds_s and cds_s <= ex[1]:
157 if ex[0] != cds_s:utr5.append((ex[0], cds_s-1))
158 break
159 else:
160 utr5.append(ex)
161 cds_e = cc[-1][1]
162 for i in range(len(ec)):
163 i += 1
164 if ec[-i][0] <= cds_e and cds_e <= ec[-i][1]:
165 if ec[-i][1] != cds_e:utr3.append((cds_e +1, ec[-i][1]))
166 break
167 else:
168 utr3.append(ec[-i])
169 utr3.reverse()
170 elif strand == '-':
171 cds_s = cc[-1][1]
172 for i in range(len(ec)):
173 i += 1
174 if ec[-i][0] <= cds_s and cds_s <= ec[-i][1]:
175 if ec[-i][1] != cds_s:utr5.append((cds_s+1, ec[-i][1]))
176 break
177 else:
178 utr5.append(ec[-i])
179 utr5.reverse()
180 cds_e = cc[0][0]
181 for ex in ec:
182 if ex[0] <= cds_e and cds_e <= ex[1]:
183 if ex[0] != cds_e:utr3.append((ex[0], cds_e-1))
184 break
185 else:
186 utr3.append(ex)
187 return utr5, utr3
188
189 def make_Exon_cod(strand_p, five_p_utr, cds_cod, three_p_utr):
190 """
191 Create exon cordinates from UTR's and CDS region
192
193 @args strand_p: feature strand information
194 @type strand_p: +/-
195 @args five_p_utr: five prime utr exon coordinates
196 @type five_p_utr: numpy array [[int, int]]
197 @args cds_cod: coding exon coordinates
198 @type cds_cod: numpy array [[int, int, int]]
199 @args three_p_utr: three prime utr exon coordinates
200 @type three_p_utr: numpy array [[int, int]]
201 """
202
203 exon_pos = []
204 if strand_p == '+':
205 utr5_start, utr5_end = 0, 0
206 if five_p_utr != []:
207 utr5_start, utr5_end = five_p_utr[-1][0], five_p_utr[-1][1]
208 cds_5start, cds_5end = cds_cod[0][0], cds_cod[0][1]
209 jun_exon = []
210 if cds_5start-utr5_end == 0 or cds_5start-utr5_end == 1:
211 jun_exon = [utr5_start, cds_5end]
212 if len(cds_cod) == 1:
213 five_prime_flag = 0
214 if jun_exon != []:
215 five_p_utr = five_p_utr[:-1]
216 five_prime_flag = 1
217 for utr5 in five_p_utr:
218 exon_pos.append(utr5)
219 jun_exon = []
220 utr3_start, utr3_end = 0, 0
221 if three_p_utr != []:
222 utr3_start = three_p_utr[0][0]
223 utr3_end = three_p_utr[0][1]
224 if utr3_start-cds_5end == 0 or utr3_start-cds_5end == 1:
225 jun_exon = [cds_5start, utr3_end]
226 three_prime_flag = 0
227 if jun_exon != []:
228 cds_cod = cds_cod[:-1]
229 three_p_utr = three_p_utr[1:]
230 three_prime_flag = 1
231 if five_prime_flag == 1 and three_prime_flag == 1:
232 exon_pos.append([utr5_start, utr3_end])
233 if five_prime_flag == 1 and three_prime_flag == 0:
234 exon_pos.append([utr5_start, cds_5end])
235 cds_cod = cds_cod[:-1]
236 if five_prime_flag == 0 and three_prime_flag == 1:
237 exon_pos.append([cds_5start, utr3_end])
238 for cds in cds_cod:
239 exon_pos.append(cds)
240 for utr3 in three_p_utr:
241 exon_pos.append(utr3)
242 else:
243 if jun_exon != []:
244 five_p_utr = five_p_utr[:-1]
245 cds_cod = cds_cod[1:]
246 for utr5 in five_p_utr:
247 exon_pos.append(utr5)
248 exon_pos.append(jun_exon) if jun_exon != [] else ''
249 jun_exon = []
250 utr3_start, utr3_end = 0, 0
251 if three_p_utr != []:
252 utr3_start = three_p_utr[0][0]
253 utr3_end = three_p_utr[0][1]
254 cds_3start = cds_cod[-1][0]
255 cds_3end = cds_cod[-1][1]
256 if utr3_start-cds_3end == 0 or utr3_start-cds_3end == 1:
257 jun_exon = [cds_3start, utr3_end]
258 if jun_exon != []:
259 cds_cod = cds_cod[:-1]
260 three_p_utr = three_p_utr[1:]
261 for cds in cds_cod:
262 exon_pos.append(cds)
263 exon_pos.append(jun_exon) if jun_exon != [] else ''
264 for utr3 in three_p_utr:
265 exon_pos.append(utr3)
266 elif strand_p == '-':
267 utr3_start, utr3_end = 0, 0
268 if three_p_utr != []:
269 utr3_start = three_p_utr[-1][0]
270 utr3_end = three_p_utr[-1][1]
271 cds_3start = cds_cod[0][0]
272 cds_3end = cds_cod[0][1]
273 jun_exon = []
274 if cds_3start-utr3_end == 0 or cds_3start-utr3_end == 1:
275 jun_exon = [utr3_start, cds_3end]
276 if len(cds_cod) == 1:
277 three_prime_flag = 0
278 if jun_exon != []:
279 three_p_utr = three_p_utr[:-1]
280 three_prime_flag = 1
281 for utr3 in three_p_utr:
282 exon_pos.append(utr3)
283 jun_exon = []
284 (utr5_start, utr5_end) = (0, 0)
285 if five_p_utr != []:
286 utr5_start = five_p_utr[0][0]
287 utr5_end = five_p_utr[0][1]
288 if utr5_start-cds_3end == 0 or utr5_start-cds_3end == 1:
289 jun_exon = [cds_3start, utr5_end]
290 five_prime_flag = 0
291 if jun_exon != []:
292 cds_cod = cds_cod[:-1]
293 five_p_utr = five_p_utr[1:]
294 five_prime_flag = 1
295 if three_prime_flag == 1 and five_prime_flag == 1:
296 exon_pos.append([utr3_start, utr5_end])
297 if three_prime_flag == 1 and five_prime_flag == 0:
298 exon_pos.append([utr3_start, cds_3end])
299 cds_cod = cds_cod[:-1]
300 if three_prime_flag == 0 and five_prime_flag == 1:
301 exon_pos.append([cds_3start, utr5_end])
302 for cds in cds_cod:
303 exon_pos.append(cds)
304 for utr5 in five_p_utr:
305 exon_pos.append(utr5)
306 else:
307 if jun_exon != []:
308 three_p_utr = three_p_utr[:-1]
309 cds_cod = cds_cod[1:]
310 for utr3 in three_p_utr:
311 exon_pos.append(utr3)
312 if jun_exon != []:
313 exon_pos.append(jun_exon)
314 jun_exon = []
315 (utr5_start, utr5_end) = (0, 0)
316 if five_p_utr != []:
317 utr5_start = five_p_utr[0][0]
318 utr5_end = five_p_utr[0][1]
319 cds_5start = cds_cod[-1][0]
320 cds_5end = cds_cod[-1][1]
321 if utr5_start-cds_5end == 0 or utr5_start-cds_5end == 1:
322 jun_exon = [cds_5start, utr5_end]
323 if jun_exon != []:
324 cds_cod = cds_cod[:-1]
325 five_p_utr = five_p_utr[1:]
326 for cds in cds_cod:
327 exon_pos.append(cds)
328 if jun_exon != []:
329 exon_pos.append(jun_exon)
330 for utr5 in five_p_utr:
331 exon_pos.append(utr5)
332 return exon_pos