getSegmentsCoordinates.py

import subprocess

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

def getSeqNumber(name_field):
    seq_number=""
    if has_numbers(name_field[-3:]):
        for char in reversed(name_field): # take the last digits of the field
            if not char.isdigit():
                break
            else:
                seq_number+=char
    else:
        for char in reversed(name_field): # take the last uppercase chars of the fied
            if not char.isupper():
                break
            else:
                seq_number+=char
    return seq_number[::-1]
    
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)
    segments = open('seg_coord/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

    return segments_size

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 # rename name_found var

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)
    walks=open('seg_coord/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]

        if check_walk_name(walk_names,name) | (len(walk_names)==1): # len=1 if there is only the source genome.

            path_start=int(line[4])
            seq_name=name.split('_')[-1]
            if (seq_name[0:10]=="chromosome") | (seq_name[0:10]=="Chromosome"):
                seq_name="Chr"+seq_name[10:]

            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)
                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
                
                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'
                out_bed.write(out_line)
                
                position+=segments_size[seg_name]
            out_bed.close()
    
    command="rm seg_coord/segments.txt && rm seg_coord/walks.txt"
    subprocess.run(command,shell=True,timeout=None)

    command="if ls test/*_MINIGRAPH_* 1> /dev/null 2>&1; then mkdir seg_coord/minigraph_segments && mv *_MINIGRAPH_* seg_coord/minigraph_segments/; fi"
    subprocess.run(command,shell=True,timeout=None)