getSegmentsCoordinates.py

import subprocess
from Graph_gff import write_line

# creates a dictionnary with the length of the segments
def get_seg_len(gfa):
    # get the lines that start with "S"
    command="grep ^S "+gfa+" > seg_coord/segments.txt"
    subprocess.run(command,shell=True,timeout=None)

    # build a dictionnary with the segment sizes
    segments_size={}
    with open('seg_coord/segments.txt','r') as seg_file:
        line=seg_file.readline()
        while line:
            line=line.split()
            seg_id='s'+line[1]
            seg_size=len(line[2])
            segments_size[seg_id]=seg_size
            line=seg_file.readline()
    return segments_size

# check that a walk in "walk_names" has "name" in it (check that "name" is a valid walk from the list)
def check_walk_name(walk_names,name):
    name_found=False
    for walk in walk_names:
        if walk in name:
            name_found=True
            break
    return name_found

# outputs the coordinates of the segments on the walks in bed files
def seg_coord(gfa,walk_names):

    # create directory to store output files
    command="mkdir seg_coord"
    subprocess.run(command,shell=True)

    segments_size=get_seg_len(gfa)

    # get the lines that start with "W"
    command="grep ^W "+gfa+" | sed 's/>/,>/g' | sed 's/</,</g' > seg_coord/walks.txt"
    subprocess.run(command,shell=True,timeout=None)

    # on these lines, get the name of the genome to name the output bed file
    file_names=list()
    with open('seg_coord/walks.txt','r') as walk_file:
        line=walk_file.readline()
        while line:
            line=line.split()
            name=line[1]+"_"+line[3]

            if (check_walk_name(walk_names,name)) or ((len(walk_names)==1) and ("MINIGRAPH" not in name)): # len=1 if there is only the source genome.
                path_start=int(line[4])
                seq_name=line[3]

                file_name="seg_coord/"+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)
                    open_mode="w"
                else :
                    open_mode="a"

                path=line[6].split(',')
                position=path_start

                with open(file_name,open_mode) as output_bed:
                    for i in range(1, len(path)): # for each segment in the path, write the position of the segment in the output bed file

                        seg_start=position
                        seg_name='s'+path[i][1:]
                        seg_stop=position+segments_size[seg_name]

                        out_line=seq_name+'\t'+str(seg_start)+'\t'+str(seg_stop)+'\t'+path[i][0:1]+seg_name+'\n'
                        output_bed.write(out_line)

                        position=seg_stop
            line=walk_file.readline()