getSegmentsCoordinates.py

import subprocess
import sys

def has_numbers(inputString):
    return any(char.isdigit() for char in inputString)

def getChrName(chromosome_field):
    chromosome_id=""
    if has_numbers(chromosome_field):
        for char in reversed(chromosome_field): # take the last digits of the field
            if not char.isdigit():
                break
            else:
                chromosome_id+=char
        chromosome_id="Chr"+chromosome_id[::-1]
    else:
        for char in reversed(chromosome_field): # take the last uppercase chars of the fied
            if not char.isupper():
                break
            else:
                chromosome_id+=char
        chromosome_id="Chr"+chromosome_id[::-1]
    return chromosome_id

if not(len(sys.argv)==2) :
        print("expected input : gfa file with walks.")
        print("output : bed files giving the coordinates of the segments on the genomes (or on minigraph segments).")
        sys.exit(1)
elif (sys.argv[1]=="-h") :
    print("expected input : gfa file with walks.")
    print("output : bed files giving the coordinates of the segments on the genomes (or on minigraph segments).")
    sys.exit(1)

gfa=sys.argv[1]

# get the lines that start with "S"
command="grep ^S "+gfa+" > segments.txt"
subprocess.run(command,shell=True)
segments = open('segments.txt','r')
lines=segments.readlines()
segments.close()

# build a dictionnary with the segment sizes
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

# get the lines that start with "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()

# on these lines, get the name of the genome to name the output bed file
file_names=list()
for line in lines :
    line=line.split()
    name=line[3]
    path_start=int(line[4])
    chromosome_field=line[3]
    chromosome_id=getChrName(chromosome_field)

    file_name=name+'.bed'

    # if we are writing in the file for the first time, overwrite it. else, append it
    # this is because chromosomes can be fragmented. the coordinates of all the fragments from the same chromosome will be written in the same bed file.
    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)): # for each segment in the path, write the position of the segment in the output bed file
        # coordinates calculation : start=position, stop=position+segment_size-1, then position+=segment_size
        
        chr='Chr'+chromosome_id[len(chromosome_id)-2:]
        
        seg_start=position
        seg_name='s'+path[i][1:]
        seg_stop=position+segments_size[seg_name]

        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'+path[i][0:1]+seg_name+'\n'
        out_bed.write(out_line)
        
        position+=segments_size[seg_name]
    out_bed.close()
 
command="rm segments.txt && rm walks.txt"
subprocess.run(command,shell=True)


command="if [ 0 -lt $(ls _MINIGRAPH_* 2>/dev/null | wc -w) ]; then mkdir minigraph_segments && mv _MINIGRAPH_* minigraph_segments; fi"
subprocess.run(command,shell=True)