Newer
Older
nina.marthe_ird.fr
committed
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
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
from Graph_gff import Features,load_intersect
from Functions import get_segment_sequence,convert_strand
target_genome_name="genome3_chr10"
intersect_path='/home/nina/annotpangenome/test_data/input_data_inf/intersect.bed'
load_intersect(intersect_path)
gfa="test_data/input_data_inf/graph_test.gfa"
def get_segments_sequence_and_paths(gfa):
file_gfa=open(gfa,'r')
lines_gfa=file_gfa.readlines()
file_gfa.close()
seg_seq={}
paths={}
for line in lines_gfa:
line=line.split()
if (line[0]=="S"): # get the sequence of the segment
seg_id='s'+line[1]
seg_seq[seg_id]=line[2]
if (line[0]=="W") & (line[1]!="_MINIGRAPH_"): # get the walk of the genome
path=line[6].replace(">",";>")
path=path.replace("<",";<").split(';')
list_path=[]
for segment in path:
if segment[0:1]=='>':
list_path.append('+s'+segment[1:])
elif segment[0:1]=='<':
list_path.append('-s'+segment[1:])
paths[line[3]]=list_path
return [paths,seg_seq]
[paths,seg_seq]=get_segments_sequence_and_paths(gfa)
segments_on_target_genome={}
pos_seg="test_data/input_data_inf/genome3_chr10.bed"
def get_segments_positions_on_genome(pos_seg):
bed=open(pos_seg,'r')
lines=bed.readlines() # read line by line ?
bed.close()
for line in lines:
line=line.split()
[seg,chrom,start,stop,strand]=[line[3][1:],line[0],line[1],line[2],line[3][0:1]]
segments_on_target_genome[seg]=[chrom,start,stop,strand]
get_segments_positions_on_genome(pos_seg)
def add_feature_sequence(feature,seg_seq):
feature_sequence=""
for segment in feature.segments_list:
if segment==feature.segments_list[0]:
feature_sequence+=get_segment_sequence(seg_seq,segment)[feature.pos_start-1:] # revérifier les +/- 1 pour la position, avec de vraies données
elif segment==feature.segments_list[-1]:
feature_sequence+=get_segment_sequence(seg_seq,segment)[0:feature.pos_stop] # revérifier les +/- 1 pour la position, avec de vraies données
else:
feature_sequence+=get_segment_sequence(seg_seq,segment)
feature.sequence=feature_sequence
def get_first_seg(list_seg):
first_seg_found=''
for segment in list_seg:
if segment[1:] in segments_on_target_genome:
first_seg_found=segment[1:]
break
return first_seg_found
def get_feature_path(paths,first_seg,last_seg):
# find the path in azucena.
first_strand=convert_strand(segments_on_target_genome[first_seg][3])
first_seg_stranded=first_strand+first_seg
last_strand=convert_strand(segments_on_target_genome[last_seg][3])
last_seg_stranded=last_strand+last_seg
index_first_seg=int(paths[target_genome_name].index(first_seg_stranded))
index_last_seg=int(paths[target_genome_name].index(last_seg_stranded))
first_index=min(index_first_seg,index_last_seg)
last_index=max(index_last_seg,index_first_seg)
list_segfeat_azu=paths[target_genome_name][first_index:last_index+1]
list_segfeat_azu_corrected=[convert_strand(segment_stranded[0])+segment_stranded[1:] for segment_stranded in list_segfeat_azu]
return list_segfeat_azu_corrected
def get_rna(dna_sequence):
return dna_sequence.replace("T","U")
# penser à transcrire la séquence codante du gène !!
def get_aa(codon):
match codon[0:2]:
case "UU":
if (codon[2]=="U") | (codon[2]=="C"):
return "Phe"
else:
return "Leu"
case "UC":
return "Ser"
case "UA":
if (codon[2]=="U") | (codon[2]=="C"):
return "Tyr"
else:
return "*"
case "UG":
if (codon[2]=="U") | (codon[2]=="C"):
return "Cys"
elif codon[2]=="A":
return "*"
else:
return "Trp"
case "CU":
return "Leu"
case "CC":
return "Pro"
case "CA":
if (codon[2]=="U") | (codon[2]=="C"):
return "His"
else:
return "Gln"
case "CG":
return "Arg"
case "AU":
if codon[2]=="G":
return "Met"
else:
return "Ile"
case "AC":
return "Thr"
case "AA":
if (codon[2]=="U") | (codon[2]=="C"):
return "Asn"
else:
return "Lys"
case "AG":
if (codon[2]=="U") | (codon[2]=="C"):
return "Ser"
else:
return "Arg"
case "GU":
return "Val"
case "GC":
return "Ala"
case "GA":
if (codon[2]=="U") | (codon[2]=="C"):
return "Asp"
else:
return "Glu"
case "GG":
return "Gly"
def traduction(sequence_arn):
list_codons=decoupe_codon(sequence_arn)
prot=list()
for codon in list_codons:
prot.append(get_aa(codon))
return prot
from textwrap import wrap
def decoupe_codon(sequence):
return wrap(sequence,3)
var=open("test_data/variations.txt",'r')
lines=var.readlines()
var.close()
# dict cds-var
cds_var={}
for line in lines:
line=line.split()
if line[1]=="CDS":
cds_id=line[0].replace('.','_').replace(':','_')
if cds_id not in cds_var.keys():
cds_var[cds_id]=list()
cds_var[cds_id].append(line)
def get_sequence_before(first_seg,seg_seq,n,paths,feat):
first_strand=convert_strand(first_seg[0])
first_seg_stranded=first_strand+first_seg[1:]
index_first_seg=int(paths[target_genome_name].index(first_seg_stranded))
sequence_before=seg_seq[first_seg[1:]][0:feat.pos_start-1] # sequence left on the segment on which the cds start (can be empty)
current_index=index_first_seg-1
while (len(sequence_before)<n) & (current_index>=0):
segment=paths[target_genome_name][current_index]
sequence_before=seg_seq[segment[1:]]+sequence_before
current_index-=1
return sequence_before[0:99]
def get_sequence_after(last_seg,seg_seq,n,paths,feat):
last_strand=convert_strand(last_seg[0])
last_seg_stranded=last_strand+last_seg[1:]
index_last_seg=int(paths[target_genome_name].index(last_seg_stranded))
sequence_after=seg_seq[last_seg[1:]][feat.pos_stop:] # sequence left on the segment on which the cds ends (can be empty)
current_index=index_last_seg+1
while (len(sequence_after)<n) & (current_index>len(paths[target_genome_name])):
segment=paths[target_genome_name][current_index]
sequence_after=sequence_after+seg_seq[segment[1:]]
current_index+=1
return sequence_after[len(sequence_after)-100:]
'''
first_strand=convert_strand(segments_on_target_genome[first_seg][3])
first_seg_stranded=first_strand+first_seg
last_strand=convert_strand(segments_on_target_genome[last_seg][3])
last_seg_stranded=last_strand+last_seg
index_first_seg=int(paths[target_genome_name].index(first_seg_stranded))
index_last_seg=int(paths[target_genome_name].index(last_seg_stranded))
first_index=min(index_first_seg,index_last_seg)
last_index=max(index_last_seg,index_first_seg)
list_segfeat_azu=paths[target_genome_name][first_index:last_index+1]
list_segfeat_azu_corrected=[convert_strand(segment_stranded[0])+segment_stranded[1:] for segment_stranded in list_segfeat_azu]
'''
for cds_id in cds_var.keys(): # for a gene that has cds, concatenate all cds to make a prot. then detail var by cds.
cds=Features[cds_id]
print("analysing variations in cds",cds_id)
add_feature_sequence(cds,seg_seq)
cds_prot=traduction(get_rna(cds.sequence))
list_seg=Features[cds_id].segments_list
first_seg=get_first_seg(list_seg)
last_seg=get_first_seg(reversed(list_seg))
path_on_target=get_feature_path(paths,first_seg,last_seg)
new_sequence=""
for segment in path_on_target:
if segment==cds.segments_list[0]:
new_sequence+=get_segment_sequence(seg_seq,segment)[cds.pos_start-1:]
elif segment==cds.segments_list[-1]:
new_sequence+=get_segment_sequence(seg_seq,segment)[0:cds.pos_stop]
else:
new_sequence+=get_segment_sequence(seg_seq,segment)
new_prot=traduction(get_rna(new_sequence))
print("original prot = ", cds_prot)
print("new prot = ", new_prot) # print new version with new start and stop codons of deleted. (before and after the gene sequence)
if new_prot[0]!="Met":
print("no Met at the start of the new version of the protein -> start codon loss")
if "Met" in new_prot:
print("Met found at position", (new_prot.index("Met")+1),"-> possible later start codon") # print cette version de la prot
print("look for start codon before. if none, print 'no start codon, likely gene not active'")
# récupérer n pb avant, les traduire, chercher la dernière Met (list[::-1].index("Met")), donner cette version de la prot
sequence_before=get_sequence_before(path_on_target[0],seg_seq,100,paths,cds)
print(sequence_before)
first_stop_index=new_prot.index("*") if "*" in new_prot else "None"
if "*" not in new_prot:
print("no stop codon")
# récupérer n pb après, les traduire, chercher le premier *, donner cette version de la prot.
elif first_stop_index+1!=len(new_prot):
print("early stop codon at position",(first_stop_index+1),"instead of",len(new_prot))
else:
print("stop codon found at expected position")
sequence_after=get_sequence_after(path_on_target[-1],seg_seq,100,paths,cds)
print(sequence_after)
for var in cds_var[cds_id]:
#print("\n",var)
size_var=int(var[11])
type_var=var[8]
pos_var=int(var[12])-1
if type_var=="insertions":
sequence_var=var[10]
if size_var%3==0:
print("pas de décalage du cadre de lecture")
trad_seq_ins=traduction(get_rna(sequence_var))
if pos_var%3==0:
print("insertion entre deux codons")
if "*" in trad_seq_ins:
print("apparition d'un codon stop")
print(f'ancienne sequence : {", ".join(cds_prot)}')
print(f'nouvelle sequence : {", ".join(cds_prot[0:(pos_var//3)-1])}, *')
else:
print(f'ancienne sequence : {", ".join(cds_prot)}')
print(f'{type_var} de {size_var//3} acides amines {", ".join(trad_seq_ins)} à la position {pos_var//3}')
else:
print("insertion au milieu d'un codon, changement de certains acides amines")
elif type_var=="deletions":
sequence_var=line[9]
if size_var%3==0:
print("pas de décalage du cadre de lecture")
trad_seq_ins=traduction(get_rna(sequence_var))
if pos_var%3==0:
print("deletion de codons entiers")
if "*" in trad_seq_ins:
print("disparition d'un codon stop")
else:
print(f'ancienne sequence : {", ".join(cds_prot)}')
print(f'{type_var} de {size_var//3} acides amines {", ".join(trad_seq_ins)} à la position {pos_var//3}')
else:
print("deletion au milieu d'un codon, changement de certains acides amines")