Mercurial > repos > jjohnson > ensembl_cdna_translate

comparison bedutil.py @ 8:5c92d0be6514 draft

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Uploaded
author jjohnson
date Thu, 14 Dec 2017 13:32:00 -0500
parents
children
comparison
equal deleted inserted replaced
7:d59e3ce10e74 8:5c92d0be6514
1 #!/usr/bin/env python
2 """
3 #
4 #------------------------------------------------------------------------------
5 # University of Minnesota
6 # Copyright 2016, Regents of the University of Minnesota
7 #------------------------------------------------------------------------------
8 # Author:
9 #
10 # James E Johnson
11 #
12 #------------------------------------------------------------------------------
13 """
14
15 import sys
16 from Bio.Seq import reverse_complement, translate
17
18
19 def bed_from_line(line,ensembl=False):
20 fields = line.rstrip('\r\n').split('\t')
21 if len(fields) < 12:
22 return None
23 (chrom, chromStart, chromEnd, name, score, strand,
24 thickStart, thickEnd, itemRgb,
25 blockCount, blockSizes, blockStarts) = fields[0:12]
26 bed_entry = BedEntry(chrom=chrom, chromStart=chromStart, chromEnd=chromEnd,
27 name=name, score=score, strand=strand,
28 thickStart=thickStart, thickEnd=thickEnd,
29 itemRgb=itemRgb,
30 blockCount=blockCount,
31 blockSizes=blockSizes.rstrip(','),
32 blockStarts=blockStarts.rstrip(','))
33 if ensembl and len(fields) >= 20:
34 bed_entry.second_name = fields[12]
35 bed_entry.cds_start_status = fields[13]
36 bed_entry.cds_end_status = fields[14]
37 bed_entry.exon_frames = fields[15].rstrip(',')
38 bed_entry.biotype = fields[16]
39 bed_entry.gene_name = fields[17]
40 bed_entry.second_gene_name = fields[18]
41 bed_entry.gene_type = fields[19]
42 return bed_entry
43
44
45
46 class BedEntry( object ):
47 def __init__(self, chrom=None, chromStart=None, chromEnd=None,
48 name=None, score=None, strand=None,
49 thickStart=None, thickEnd=None, itemRgb=None,
50 blockCount=None, blockSizes=None, blockStarts=None):
51 self.chrom = chrom
52 self.chromStart = int(chromStart)
53 self.chromEnd = int(chromEnd)
54 self.name = name
55 self.score = int(score) if score is not None else 0
56 self.strand = '-' if str(strand).startswith('-') else '+'
57 self.thickStart = int(thickStart) if thickStart else self.chromStart
58 self.thickEnd = int(thickEnd) if thickEnd else self.chromEnd
59 self.itemRgb = str(itemRgb) if itemRgb is not None else r'100,100,100'
60 self.blockCount = int(blockCount)
61 if isinstance(blockSizes, str) or isinstance(blockSizes, unicode):
62 self.blockSizes = [int(x) for x in blockSizes.split(',')]
63 elif isinstance(blockSizes, list):
64 self.blockSizes = [int(x) for x in blockSizes]
65 else:
66 self.blockSizes = blockSizes
67 if isinstance(blockStarts, str) or isinstance(blockSizes, unicode):
68 self.blockStarts = [int(x) for x in blockStarts.split(',')]
69 elif isinstance(blockStarts, list):
70 self.blockStarts = [int(x) for x in blockStarts]
71 else:
72 self.blockStarts = blockStarts
73 self.second_name = None
74 self.cds_start_status = None
75 self.cds_end_status = None
76 self.exon_frames = None
77 self.biotype = None
78 self.gene_name = None
79 self.second_gene_name = None
80 self.gene_type = None
81 self.seq = None
82 self.cdna = None
83 self.pep = None
84 # T26C
85 self.aa_change = []
86 # p.Trp26Cys g.<pos><ref>><alt> # g.1304573A>G
87 self.variants = []
88
89
90 def __str__(self):
91 return '%s\t%d\t%d\t%s\t%d\t%s\t%d\t%d\t%s\t%d\t%s\t%s' % (
92 self.chrom, self.chromStart, self.chromEnd,
93 self.name, self.score, self.strand,
94 self.thickStart, self.thickEnd, str(self.itemRgb), self.blockCount,
95 ','.join([str(x) for x in self.blockSizes]),
96 ','.join([str(x) for x in self.blockStarts]))
97
98
99 def get_splice_junctions(self):
100 splice_juncs = []
101 for i in range(self.blockCount - 1):
102 splice_junc = "%s:%d_%d" % (self.chrom, self.chromStart + self.blockSizes[i], self.chromStart + self.blockStarts[i+1])
103 splice_juncs.append(splice_junc)
104 return splice_juncs
105
106
107 def get_exon_seqs(self):
108 if not self.seq:
109 return None
110 exons = []
111 for i in range(self.blockCount):
112 # splice_junc = "%s:%d_%d" % (self.chrom, self.chromStart + self.blockSizes[i], self.chromStart + self.blockStarts[i+1])
113 exons.append(self.seq[self.blockStarts[i]:self.blockStarts[i] + self.blockSizes[i]])
114 if self.strand == '-': #reverse complement
115 exons.reverse()
116 for i,s in enumerate(exons):
117 exons[i] = reverse_complement(s)
118 return exons
119
120
121 def get_spliced_seq(self,strand=None):
122 if not self.seq:
123 return None
124 seq = ''.join(self.get_exon_seqs())
125 if strand and self.strand != strand:
126 seq = reverse_complement(seq)
127 return seq
128
129
130 def get_cdna(self):
131 if not self.cdna:
132 self.cdna = self.get_spliced_seq()
133 return self.cdna
134
135
136 def get_cds(self):
137 cdna = self.get_cdna()
138 if cdna:
139 if self.chromStart == self.thickStart and self.chromEnd == self.thickEnd:
140 return cdna
141 pos = [self.offset_of_pos(self.thickStart),self.offset_of_pos(self.thickEnd)]
142 if 0 <= min(pos) <= max(pos) <= len(cdna):
143 return cdna[min(pos):max(pos)]
144 return None
145
146
147 def get_cigar(self):
148 cigar = ''
149 md = ''
150 r = range(self.blockCount)
151 rev = self.strand == '-'
152 ## if rev: r.reverse()
153 xl = None
154 for x in r:
155 if xl is not None:
156 intronSize = abs(self.blockStarts[x] - self.blockSizes[xl] - self.blockStarts[xl])
157 cigar += '%dN' % intronSize
158 cigar += '%dM' % self.blockSizes[x]
159 xl = x
160 return cigar
161
162
163 def get_cigar_md(self):
164 cigar = ''
165 md = ''
166 r = range(self.blockCount)
167 rev = self.strand == '-'
168 ## if rev: r.reverse()
169 xl = None
170 for x in r:
171 if xl is not None:
172 intronSize = abs(self.blockStarts[x] - self.blockSizes[xl] - self.blockStarts[xl])
173 cigar += '%dN' % intronSize
174 cigar += '%dM' % self.blockSizes[x]
175 xl = x
176 md = '%d' % sum(self.blockSizes)
177 return (cigar,md)
178
179
180 def get_translation(self,sequence=None):
181 translation = None
182 seq = sequence if sequence else self.get_spliced_seq()
183 if seq:
184 seqlen = len(seq) / 3 * 3;
185 if seqlen >= 3:
186 translation = translate(seq[:seqlen])
187 return translation
188
189
190 def get_translations(self):
191 translations = []
192 seq = self.get_spliced_seq()
193 if seq:
194 for i in range(3):
195 translation = self.get_translation(sequence=seq[i:])
196 if translation:
197 translations.append(translation)
198 return translations
199
200 def pos_of_na(self, offset):
201 if offset is not None and 0 <= offset < sum(self.blockSizes):
202 r = range(self.blockCount)
203 rev = self.strand == '-'
204 if rev:
205 r.reverse()
206 nlen = 0
207 for x in r:
208 if offset < nlen + self.blockSizes[x]:
209 if rev:
210 return self.chromStart + self.blockStarts[x] + self.blockSizes[x] - (offset - nlen)
211 else:
212 return self.chromStart + self.blockStarts[x] + (offset - nlen)
213 nlen += self.blockSizes[x]
214 return None
215
216 def offset_of_pos(self, pos):
217 if not self.chromStart <= pos < self.chromEnd:
218 return -1
219 r = range(self.blockCount)
220 rev = self.strand == '-'
221 if rev:
222 r.reverse()
223 nlen = 0
224 for x in r:
225 bStart = self.chromStart + self.blockStarts[x]
226 bEnd = bStart + self.blockSizes[x]
227 ## print >> sys.stdout, "offset_of_pos %d %d %d %d" % (bStart,pos,bEnd,nlen)
228 if bStart <= pos < bEnd:
229 return nlen + (bEnd - pos if rev else pos - bStart)
230 nlen += self.blockSizes[x]
231
232 def apply_variant(self,pos,ref,alt):
233 pos = int(pos)
234 if not ref or not alt:
235 print >> sys.stderr, "variant requires ref and alt sequences"
236 return
237 if not self.chromStart <= pos <= self.chromEnd:
238 print >> sys.stderr, "variant not in entry %s: %s %d < %d < %d" % (self.name,self.strand, self.chromStart,pos,self.chromEnd)
239 print >> sys.stderr, "%s" % str(self)
240 return
241 if len(ref) != len(alt):
242 print >> sys.stderr, "variant only works for snp: %s %s" % (ref,alt)
243 return
244 if not self.seq:
245 print >> sys.stderr, "variant entry %s has no seq" % self.name
246 return
247 """
248 if self.strand == '-':
249 ref = reverse_complement(ref)
250 alt = reverse_complement(alt)
251 """
252 bases = list(self.seq)
253 offset = pos - self.chromStart
254 for i in range(len(ref)):
255 # offset = self.offset_of_pos(pos+i)
256 if offset is not None:
257 bases[offset+i] = alt[i]
258 else:
259 print >> sys.stderr, "variant offset %s: %s %d < %d < %d" % (self.name,self.strand,self.chromStart,pos+1,self.chromEnd)
260 print >> sys.stderr, "%s" % str(self)
261 self.seq = ''.join(bases)
262 self.variants.append("g.%d%s>%s" % (pos+1,ref,alt))
263
264 def get_variant_bed(self,pos,ref,alt):
265 pos = int(pos)
266 if not ref or not alt:
267 print >> sys.stderr, "variant requires ref and alt sequences"
268 return None
269 if not self.chromStart <= pos <= self.chromEnd:
270 print >> sys.stderr, "variant not in entry %s: %s %d < %d < %d" % (self.name,self.strand, self.chromStart,pos,self.chromEnd)
271 print >> sys.stderr, "%s" % str(self)
272 return None
273 if not self.seq:
274 print >> sys.stderr, "variant entry %s has no seq" % self.name
275 return None
276 tbed = BedEntry(chrom = self.chrom, chromStart = self.chromStart, chromEnd = self.chromEnd,
277 name = self.name, score = self.score, strand = self.strand,
278 thickStart = self.chromStart, thickEnd = self.chromEnd, itemRgb = self.itemRgb,
279 blockCount = self.blockCount, blockSizes = self.blockSizes, blockStarts = self.blockStarts)
280 bases = list(self.seq)
281 offset = pos - self.chromStart
282 tbed.seq = ''.join(bases[:offset] + list(alt) + bases[offset+len(ref):])
283 if len(ref) != len(alt):
284 diff = len(alt) - len(ref)
285 rEnd = pos + len(ref)
286 aEnd = pos + len(alt)
287 #need to adjust blocks
288 # change spans blocks,
289 for x in range(tbed.blockCount):
290 bStart = tbed.chromStart + tbed.blockStarts[x]
291 bEnd = bStart + tbed.blockSizes[x]
292 # change within a block or extends (last block), adjust blocksize
293 # seq: GGGcatGGG
294 # ref c alt tag: GGGtagatGGG
295 # ref cat alt a: GGGaGGG
296 if bStart <= pos < rEnd < bEnd:
297 tbed.blockSizes[x] += diff
298 return tbed
299
300 ## (start,end)
301 def get_subrange(self,tstart,tstop,debug=False):
302 chromStart = self.chromStart
303 chromEnd = self.chromEnd
304 if debug:
305 print >> sys.stderr, "%s" % (str(self))
306 r = range(self.blockCount)
307 if self.strand == '-':
308 r.reverse()
309 bStart = 0
310 bEnd = 0
311 for x in r:
312 bEnd = bStart + self.blockSizes[x]
313 if bStart <= tstart < bEnd:
314 if self.strand == '+':
315 chromStart = self.chromStart + self.blockStarts[x] +\
316 (tstart - bStart)
317 else:
318 chromEnd = self.chromStart + self.blockStarts[x] +\
319 self.blockSizes[x] - (tstart - bStart)
320 if bStart <= tstop < bEnd:
321 if self.strand == '+':
322 chromEnd = self.chromStart + self.blockStarts[x] +\
323 (tstop - bStart)
324 else:
325 chromStart = self.chromStart + self.blockStarts[x] +\
326 self.blockSizes[x] - (tstop - bStart)
327 if debug:
328 print >> sys.stderr,\
329 "%3d %s\t%d\t%d\t%d\t%d\t%d\t%d"\
330 % (x, self.strand, bStart, bEnd,
331 tstart, tstop, chromStart, chromEnd)
332 bStart += self.blockSizes[x]
333 return(chromStart, chromEnd)
334
335
336 # get the blocks for sub range
337 def get_blocks(self, chromStart, chromEnd):
338 tblockCount = 0
339 tblockSizes = []
340 tblockStarts = []
341 for x in range(self.blockCount):
342 bStart = self.chromStart + self.blockStarts[x]
343 bEnd = bStart + self.blockSizes[x]
344 if bStart > chromEnd:
345 break
346 if bEnd < chromStart:
347 continue
348 cStart = max(chromStart, bStart)
349 tblockStarts.append(cStart - chromStart)
350 tblockSizes.append(min(chromEnd, bEnd) - cStart)
351 tblockCount += 1
352 return (tblockCount, tblockSizes, tblockStarts)
353
354
355 def trim(self, tstart, tstop, debug=False):
356 (tchromStart, tchromEnd) =\
357 self.get_subrange(tstart, tstop, debug=debug)
358 (tblockCount, tblockSizes, tblockStarts) =\
359 self.get_blocks(tchromStart, tchromEnd)
360 tbed = BedEntry(
361 chrom=self.chrom, chromStart=tchromStart, chromEnd=tchromEnd,
362 name=self.name, score=self.score, strand=self.strand,
363 thickStart=tchromStart, thickEnd=tchromEnd, itemRgb=self.itemRgb,
364 blockCount=tblockCount,
365 blockSizes=tblockSizes, blockStarts=tblockStarts)
366 if self.seq:
367 ts = tchromStart-self.chromStart
368 te = tchromEnd - tchromStart + ts
369 tbed.seq = self.seq[ts:te]
370 return tbed
371
372
373 def get_filtered_translations(self,untrimmed=False,filtering=True,ignore_left_bp=0,ignore_right_bp=0,debug=False):
374 translations = [None,None,None]
375 seq = self.get_spliced_seq()
376 ignore = (ignore_left_bp if self.strand == '+' else ignore_right_bp) / 3
377 block_sum = sum(self.blockSizes)
378 exon_sizes = [x for x in self.blockSizes]
379 if self.strand == '-':
380 exon_sizes.reverse()
381 splice_sites = [sum(exon_sizes[:x]) / 3 for x in range(1,len(exon_sizes))]
382 if debug:
383 print >> sys.stderr, "splice_sites: %s" % splice_sites
384 junc = splice_sites[0] if len(splice_sites) > 0 else exon_sizes[0]
385 if seq:
386 for i in range(3):
387 translation = self.get_translation(sequence=seq[i:])
388 if translation:
389 tstart = 0
390 tstop = len(translation)
391 offset = (block_sum - i) % 3
392 if debug:
393 print >> sys.stderr, "frame: %d\ttstart: %d tstop: %d offset: %d\t%s" % (i,tstart,tstop,offset,translation)
394 if not untrimmed:
395 tstart = translation.rfind('*',0,junc) + 1
396 stop = translation.find('*',junc)
397 tstop = stop if stop >= 0 else len(translation)
398 offset = (block_sum - i) % 3
399 trimmed = translation[tstart:tstop]
400 if debug:
401 print >> sys.stderr, "frame: %d\ttstart: %d tstop: %d offset: %d\t%s" % (i,tstart,tstop,offset,trimmed)
402 if filtering and tstart > ignore:
403 continue
404 #get genomic locations for start and end
405 if self.strand == '+':
406 chromStart = self.chromStart + i + (tstart * 3)
407 chromEnd = self.chromEnd - offset - (len(translation) - tstop) * 3
408 else:
409 chromStart = self.chromStart + offset + (len(translation) - tstop) * 3
410 chromEnd = self.chromEnd - i - (tstart * 3)
411 #get the blocks for this translation
412 (tblockCount,tblockSizes,tblockStarts) = self.get_blocks(chromStart,chromEnd)
413 translations[i] = (chromStart,chromEnd,trimmed,tblockCount,tblockSizes,tblockStarts)
414 if debug:
415 print >> sys.stderr, "tblockCount: %d tblockStarts: %s tblockSizes: %s" % (tblockCount,tblockStarts,tblockSizes)
416 # translations[i] = (chromStart,chromEnd,trimmed,tblockCount,tblockSizes,tblockStarts)
417 return translations
418
419
420 def get_seq_id(self,seqtype='unk:unk',reference='',frame=None):
421 ## Ensembl fasta ID format
422 # >ID SEQTYPE:STATUS LOCATION GENE TRANSCRIPT
423 # >ENSP00000328693 pep:splice chromosome:NCBI35:1:904515:910768:1 gene:ENSG00000158815:transcript:ENST00000328693 gene_biotype:protein_coding transcript_biotype:protein_coding
424 frame_name = ''
425 chromStart = self.chromStart
426 chromEnd = self.chromEnd
427 strand = 1 if self.strand == '+' else -1
428 if frame != None:
429 block_sum = sum(self.blockSizes)
430 offset = (block_sum - frame) % 3
431 frame_name = '_' + str(frame + 1)
432 if self.strand == '+':
433 chromStart += frame
434 chromEnd -= offset
435 else:
436 chromStart += offset
437 chromEnd -= frame
438 location = "chromosome:%s:%s:%s:%s:%s" % (reference,self.chrom,chromStart,chromEnd,strand)
439 seq_id = "%s%s %s %s" % (self.name,frame_name,seqtype,location)
440 return seq_id
441
442
443 def get_line(self, start_offset = 0, end_offset = 0):
444 if start_offset or end_offset:
445 s_offset = start_offset if start_offset else 0
446 e_offset = end_offset if end_offset else 0
447 if s_offset > self.chromStart:
448 s_offset = self.chromStart
449 chrStart = self.chromStart - s_offset
450 chrEnd = self.chromEnd + e_offset
451 blkSizes = self.blockSizes
452 blkSizes[0] += s_offset
453 blkSizes[-1] += e_offset
454 blkStarts = self.blockStarts
455 for i in range(1,self.blockCount):
456 blkStarts[i] += s_offset
457 items = [str(x) for x in [self.chrom,chrStart,chrEnd,self.name,self.score,self.strand,self.thickStart,self.thickEnd,self.itemRgb,self.blockCount,','.join([str(x) for x in blkSizes]),','.join([str(x) for x in blkStarts])]]
458 return '\t'.join(items) + '\n'
459 return self.line