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