Newer
Older
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
# input : gfa file with paths.
gfa="platGB2_r426_chr10_cactus.gfa"
# output : bed files giving the coordinates of the segments on the genomes (/minigraph segments)
import subprocess
# récupérer les lignes qui commencent par 'S'
command="grep ^S "+gfa+" > segments.txt"
subprocess.run(command,shell=True)
segments = open('segments.txt','r')
lines=segments.readlines()
segments.close()
# faire un dictionnaire avec les tailles des segments
segments_size={}
for line in lines:
line=line.split()
seg_id='s'+line[1]
seg_size=len(line[2])
segments_size[seg_id]=seg_size
#print(segments_size['1'])
# récupérer les lignes qui commencent par 'P'
command="grep ^W "+gfa+" | sed 's/>/,>/g' | sed 's/</,</g' > walks.txt"
subprocess.run(command,shell=True)
walks=open('walks.txt','r')
lines=walks.readlines()
walks.close()
# sur ces lignes, récupérer le nom du génome pour nommer le bed
file_names=list()
for line in lines :
line=line.split()
name=line[1]
chr_field=line[3].split('_')
chromosome=chr_field[len(chr_field)-1]
path_start=int(line[4])
path_stop=int(line[5])
file_name=name+'_'+chromosome+'.bed'
if file_name not in file_names:
file_names.append(file_name)
out_bed = open(file_name, 'w')
else :
out_bed = open(file_name, 'a')
# si le fichier est créé pour la premiere fois
path=line[6].split(',')
seg_cpt=1
position=path_start
for i in range(1, len(path)): # pour chaque segment du path, print dans le bed la position du segment.
# calcul coordonnées : start=position, stop=position+taille_segment-1, ensuite position+=taille_segment
seg_start=position
seg_name='s'+path[i][1:]
seg_stop=position+segments_size[seg_name]-1
if path[i][0:1]=='>':
orientation='+'
elif path[i][0:1]=='<':
orientation='-'
else:
orientation='.'
out_line=seg_name+' '+str(seg_start)+' '+str(seg_stop)+' '+orientation+'\n'
out_bed.write(out_line)
position+=segments_size[seg_name]
out_bed.close()