# 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 'W' 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]) file_name=name+'_'+chromosome+'.bed' # si on écrit dans le fichier pour la première fois, l'écraser. sinon l'append. # certains chr d'un même individu sont fragmentés. les coordonnées des segments de tous les fragments sont inscrites dans le même 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') path=line[6].split(',') 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 chr='Chr'+chromosome[len(chromosome)-2:] 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=chr+'\t'+str(seg_start)+'\t'+str(seg_stop)+'\t'+orientation+'\t'+seg_name+'\n' out_bed.write(out_line) position+=segments_size[seg_name] out_bed.close()