script pour récupérer les coordonnées des segments sur les génomes de...

script pour récupérer les coordonnées des segments sur les génomes de référence à partir d'un gfa avec walks

seg_coord/getSegmentsCoordinates.py

+# 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()
+ 
+
+