-
Notifications
You must be signed in to change notification settings - Fork 2
/
Copy pathtranscript.py
134 lines (107 loc) · 4.4 KB
/
transcript.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
class Transcript:
'''
'''
def __init__(self, chrom, start, end, cov, name):
self.chrom = chrom
self.start = start
self.end = end
self.cov = cov
self.exons = []
self.name = name
def scoreTranscript(self, transcript, threshold):
''' Compare this transcript to the given transcript and return a score representing their closeness
0 = no match, 1 = perfect match
'''
if (not self.chrom == transcript.chrom) or (self.start > transcript.end) or (self.end < transcript.start):
#if ((self.start > transcript.end) or (self.end < transcript.start)):
return 0
thisMatches = []
scores = []
for i, e1 in enumerate(self.exons):
closest = 0
closestScore = 0
for j in range(len(transcript.exons)):
e2 = transcript.exons[j]
currScore = self.scoreExons(e1, e2, threshold)
if currScore > closestScore:
#logging.info('scoreExons returned %0.2f' % currScore)
closestScore = currScore
closest = j
thisMatches.append(closest)
scores.append(closestScore)
otherMatches = []
for i in range(len(transcript.exons)):
e1 = transcript.exons[i]
closest = 0
closestScore = 0
for j in range(len(self.exons)):
e2 = self.exons[j]
currScore = self.scoreExons(e1, e2, threshold)
if currScore > closestScore:
#logging.info('scoreExons returned %0.2f' % currScore)
closestScore = currScore
closest = j
otherMatches.append(closest)
totalScore = 0.0
count = len(thisMatches)
for i in range(len(thisMatches)):
if otherMatches[thisMatches[i]] == i:
totalScore += scores[i]
for i in range(len(otherMatches)):
if not thisMatches[otherMatches[i]] == i:
count += 1
return totalScore / float(count)
def scoreTranscript2(self, transcript, threshold):
if (not self.chrom == transcript.chrom) or ((self.start > transcript.end) or (self.end < transcript.start)):
return 0
if not len(self.exons) == len(transcript.exons):
return 0
score = 1
for i in range(len(self.exons)):
s = self.scoreExons(self.exons[i], transcript.exons[i], threshold)
if s == 0:
return 0
else:
score *= s
return score
def scoreTranscript3(self, transcript, threshold):
#Same as scoreTranscript() above but calculate score as proportion of overlap
if (not self.chrom == transcript.chrom) or ((self.start > transcript.end) or (self.end < transcript.start)):
return 0
overlap = 0
i = 0
for e1 in self.exons:
while i < len(transcript.exons) and transcript.exons[i][1] <= e1[0]:
i += 1
if i == len(transcript.exons):
break
while i < len(transcript.exons) and transcript.exons[i][0] < e1[1]:
overlap += min(transcript.exons[i][1], e1[1]) - max(transcript.exons[i][0], e1[0])
if transcript.exons[i][0] < e1[1]:
i += 1
len1 = 0
for e1 in transcript.exons:
len1 += e1[1] - e1[0]
len2 = 0
for e2 in transcript.exons:
len2 += e2[1] - e2[0]
if overlap > len1 or overlap > len2:
print(self.exons)
print(transcript.exons)
print('Length 1 = %d' % len1)
print('Length 2 = %d' % len2)
print('Overlap = %d' % overlap)
exit()
return float(overlap * overlap) / float(len1 * len2)
def scoreExons(self, exon1, exon2, threshold):
''' Returns a value between 0 and 1, where 1 indicates that the 2 exons are a perfect match and 0 indicates no match
'''
if threshold == 0:
if exon1[0] == exon2[0] and exon1[1] == exon2[1]:
return 1
else:
return 0
else:
dx = min(abs(exon1[0] - exon2[0]), threshold)
dy = min(abs(exon1[1] - exon2[1]), threshold)
return 1 - dx / (2.0*threshold) - dy / (2.0*threshold)