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