8
|
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
|