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