-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy patheliot_main.py~
287 lines (218 loc) · 9.03 KB
/
eliot_main.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
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
from typing import *
class ORF:
"""Open Reading Frame defined by start and stop nucleotide"""
def __init__(self, source: str, start: int, stop: int, code_table_id: int, protein: str, strand: str):
##########################
# À corriger / Améliorer #
##########################
self.source = source
self.strand = strand
if self.strand == "+":
self.start = start
self.stop = stop
else:
self.start = len(source) - start
self.stop = len(source) - stop
self.code_table = code_table_id
self.protein = protein
self.length = abs(self.stop - self.start)
self.product = None # A rajouter
self.name = None # A rajouter
def __len__(self):
return self.length
def __repr__(self):
return "{} : {}..{}".format(self.strand, self.start, self.stop)
def as_dict(self):
return {"strand": self.strand,
"start": self.start,
"stop": self.stop,
"transl_table": self.code_table,
"protein": self.protein
}
def find_orf(seq: str, threshold: int, code_table_id: int) -> [ORF]:
"""Give a list of all ORF in the sequence if they are grater than the threshold
This function is written by Eliot Ragueneau.
Args:
seq: Sequence source
threshold: Minimum size of the ORF in the list
code_table_id: NCBI identifier of the translation table used on this sequence
Returns:
list of ORF
"""
##########################
# À corriger / Améliorer #
##########################
transl_table, start_table = get_genetic_code(code_table_id)
strands = {"+": seq, "-": reversed_complement(seq)}
orf_list = []
for strand in strands:
dict_init = {}
for i in range(len(strands[strand])):
if strands[strand][i:i + 3] in start_table:
dict_init[i] = None
list_stop = [] # Seulement pour n'avoir que les ORFs maximum
for init in dict_init:
prot = ""
for i in range(init + 3, len(strands[strand]), 3):
codon = strands[strand][i - 3: i]
if len(codon) == 3:
aa = transl_table[codon]
if aa == "*":
if i - init > threshold and i not in list_stop: # Seulement pour n'avoir que les ORFs maximum
list_stop.append(i)
orf_list.append(ORF(seq,
start=init,
stop=i,
code_table_id=code_table_id,
protein=prot,
strand=strand))
break
prot += aa
else:
break
return orf_list
def get_genetic_code(ncbi_id: int) -> Tuple[dict, dict]:
"""Give the initiation codons and the translation table by its id
This function is written by Eliot Ragueneau.
Args:
ncbi_id: ncbi translation table identifier
Returns:
transl_table: translating dictionary
start_table: list of start codons
"""
with open("Translation_tables/{}.txt".format(ncbi_id), "r") as file:
lines = [line.strip() for line in file.readlines()]
transl_table = {lines[2][i] + lines[3][i] + lines[4][i]: lines[0][i]
for i in range(len(lines[0]))}
start_table = {lines[2][i] + lines[3][i] + lines[4][i]: lines[1][i]
for i in range(len(lines[0])) if lines[1][i] == "M"}
return transl_table, start_table
def get_lengths(orf_list: list) -> [int]:
"""Give the list of orf lengths from a list of orf
This function is written by Eliot Ragueneau.
Args:
orf_list: list of ORF.
Returns:
list of lengths of ORF
"""
return [len(orf) for orf in orf_list]
def get_longest_orf(orf_list: List[ORF]) -> ORF:
"""Give the longest orf from a list of orf
This function is written by Eliot Ragueneau.
Args:
orf_list: list of ORF.
Returns:
longest ORF
"""
return max(orf_list, key=lambda orf: len(orf))
def get_top_longest_orf(orf_list: List[ORF], value: float) -> [ORF]:
"""Return the value% top longest orfs from a list of orf
This function is written by Eliot Ragueneau.
Args:
orf_list: list of ORF.
value: 0 > float > 1 : % of the top longest orf to show
Returns:
list of top ORF
"""
orf_list.sort(key=lambda orf: len(orf))
return orf_list[int(- value * len(orf_list)) - 1:]
def reversed_complement(dna_seq: str):
"""Return the reversed complement of the given DNA sequence
This function is written by Eliot Ragueneau.
Args:
dna_seq: DNA sequence to be reversed.
Returns:
reversed complement DNA sequence
"""
complement_dict = {"A": "T", "T": "A", "C": "G", "G": "C"}
return "".join([complement_dict[i] for i in dna_seq[::-1]])
def read_csv(filename: str, separator: str = ";") -> [dict]:
with open(filename, 'r') as file:
list_lines = [line.strip() for line in file.readlines()]
keys = list_lines[0].split(separator)
return [{key: element for key in keys for element in line.split(separator)} for line in list_lines[1:]]
def write_csv(filename: str, data: list, separator: str = ";"):
if isinstance(data[0], ORF):
data = [orf.as_dict() for orf in data]
filename += ".csv" if filename[-4:] != ".csv" else ""
with open(filename, "w") as file:
file.write("{}\n".format(separator.join(data[0]))) # Write only the keys
for dictionary in data:
file.write("{}\n".format(separator.join([str(dictionary[key]) for key in dictionary])))
class GenBank:
"""GenBank data structure (More interesting than dictionaries)"""
def __init__(self, filename: str):
"""Parse a GenBank file
This function is written by .
Args:
filename: .gb file to parse
Returns:
list of ORF (either as ORF or as dict)
"""
text = self.read_flat_file(filename)
features = self.get_features(text)
self.genes = self.get_genes(features)
self.definition = None
self.type = None # dna, rna, or protein
self.data = None # only if available otherwise set to ‘xxx’. if seq too large, the entry does not contain data.
self.id = None
self.length = None
self.gbtype = None
self.organism = None
self.code_table_id = None
@staticmethod
def read_flat_file(filename: str) -> str:
"""Load a file in memory by returning a string
This function is written by .
Args:
filename: file to open
Returns:
string of the whole file (with \n)
"""
# Voir with statements sur Sam et Max
pass
@staticmethod
def get_features(txt: str) -> str:
"""Extract features lines from flat text and return them
This function is written by .
Args:
txt: flat text with features to extract
Returns:
string of features
"""
# Voir package builtin re : peut être utile ici
pass
@staticmethod
def get_genes(features: str) -> [ORF]:
"""Extract gene and CDS data from their section inside features table
This function is written by .
Args:
features: text with gene to extract
Returns:
list of ORF (either as ORF or as dict)
"""
# Voir package builtin re : peut être utile ici
pass
@staticmethod
def read(filename):
"""Parse a GenBank file
This function is written by .
Args:
filename: .gb file to parse
Returns:
list of ORF (either as ORF or as dict)
"""
return GenBank(filename)
def read_fasta(filename: str) -> str:
dna = ""
with open(filename, 'r') as fasta:
for fasta_line in fasta:
if fasta_line[0] != ">":
dna += fasta_line.strip()
if __name__ == '__main__':
dna = read_fasta("influenza.fasta")
list_of_orf = find_orf(dna[:50001], 300, 11)
list_of_orf.sort(key=lambda orf: orf.start, reverse=True)
print(get_lengths(list_of_orf))
print(get_longest_orf(list_of_orf))
print(get_top_longest_orf(list_of_orf, 0.2))