-
nina.marthe_ird.fr authorednina.marthe_ird.fr authored
getSegmentsCoordinates.py 2.71 KiB
import subprocess
from Graph_gff import write_line
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[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)
out_bed = open(file_name, 'w')
else :
out_bed = open(file_name, 'a')
output_bed=[0,"",out_bed]
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
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'
write_line(out_line,output_bed,False)
position=seg_stop
write_line("",output_bed,True)
out_bed.close()