adapted the code to keep the orientation of the segments from the walks. added...

adapted the code to keep the orientation of the segments from the walks. added functions to handle oriented segments

Graph_gff.py

@@ -23,6 +23,25 @@ def write_line(line,output,force):
         output[1]=""
 
 
+def invert_seg(seg):
+    if seg[0]==">":
+        inv_seg="<"+seg[1:]
+    elif seg[0]=="<":
+        inv_seg=">"+seg[1:]
+    else:
+        print(seg," not invertable")
+        return seg
+    return inv_seg
+
+
+def search_segment(segment):
+    if segment not in Segments:
+        if invert_seg(segment) in Segments:
+            return invert_seg(segment)
+        return False
+    return segment
+
+
 # functions to create the segments and the features
 
 # create a segment with its first feature
@@ -73,8 +92,8 @@ def add_seg(feat_id,new_seg_oriented):
 
 # add the position of the feature on its first and last segment
 def add_pos(feature):
-    feature.pos_start=get_feature_start_on_segment(feature.segments_list_source[0][1:],feature.id)
-    feature.pos_stop=get_feature_stop_on_segment(feature.segments_list_source[-1][1:],feature.id)
+    feature.pos_start=get_feature_start_on_segment(feature.segments_list_source[0],feature.id)
+    feature.pos_stop=get_feature_stop_on_segment(feature.segments_list_source[-1],feature.id)
 
 def load_feature_intersect(line,feature_id,seg_oriented):
     if feature_id not in Features: # if the feature doesn't exist, create it and add the current segment to its seg list
@@ -123,7 +142,7 @@ def load_intersect(intersect_path):
         line=line.split()
         # get the ids for the dictionnaries' keys
         feature_id=line[12].split(';')[0].split("=")[1].replace(".","_").replace(":","_")
-        segment_id=line[3][1:]
+        segment_id=line[3]
         strand=line[10]
         segment_oriented=line[3]
         feature_stranded=strand+feature_id
@@ -146,7 +165,7 @@ def load_intersect(intersect_path):
 
 # get the feature's start position on the segment
 def get_feature_start_on_segment(seg_id,feat_id):
-    s=Segments[seg_id]
+    s=Segments[search_segment(seg_id)]
     f=Features[feat_id]
     if s.start>=f.start:
         result=1
@@ -156,7 +175,7 @@ def get_feature_start_on_segment(seg_id,feat_id):
 
 # get the feature's stop position on the segment
 def get_feature_stop_on_segment(seg_id,feat_id):
-    s=Segments[seg_id]
+    s=Segments[search_segment(seg_id)]
     f=Features[feat_id]
     if s.stop<=f.stop:
         result=s.size
@@ -165,30 +184,24 @@ def get_feature_stop_on_segment(seg_id,feat_id):
     return result
 
 # generates the annotation for the gff of the graph, with the rank of the segment in the feature
-def make_annot_graph(feature,segment_id):
-
-    if ">"+segment_id in feature.segments_list_source:
-        segment_oriented=">"+segment_id
-    else:
-        segment_oriented="<"+segment_id
-        # what if the segment is present twice in the list, with the two orientation ??
+def make_annot_graph(feature,segment):
 
     # get the rank and the total number of ocurrences for the feature
-    rank=str(feature.segments_list_source.index(segment_oriented)+1)
+    rank=str(feature.segments_list_source.index(segment)+1)
     total=str(len(feature.segments_list_source))
     # create the annotation with the rank information
     annotation=feature.annot+";Rank_occurrence="+rank+";Total_occurrences="+total
     return annotation
 
-def create_line_graph_gff(feature_stranded,segment_id):
+def create_line_graph_gff(feature_stranded,segment):
     [strand,feature_id]=[feature_stranded[0:1],feature_stranded[1:]]
     feature=Features[feature_id]
 
-    annotation=make_annot_graph(feature,segment_id)
+    annotation=make_annot_graph(feature,segment)
     type=feature.type
-    start=get_feature_start_on_segment(segment_id,feature_id)
-    stop=get_feature_stop_on_segment(segment_id,feature_id)
-    line=segment_id+"\tGrAnnot\t"+type+"\t"+str(start)+"\t"+str(stop)+"\t.\t"+strand+"\t.\t"+annotation+"\n"
+    start=get_feature_start_on_segment(segment,feature_id)
+    stop=get_feature_stop_on_segment(segment,feature_id)
+    line=segment+"\tGrAnnot\t"+type+"\t"+str(start)+"\t"+str(stop)+"\t.\t"+strand+"\t.\t"+annotation+"\n"
     return line
 
 # go through all the segments in Segments and prints the gff, with one line for each segment/feature intersection
@@ -197,10 +210,10 @@ def graph_gff(graph_gff_file):
     graph_gff_file = open(graph_gff_file, 'w')
     output_graph_gff[2]=graph_gff_file
 
-    for segment_id in Segments:
-        features_on_seg=Segments[segment_id].features # get the list of the features on the segment
+    for segment_stranded in Segments:
+        features_on_seg=Segments[segment_stranded].features # get the list of the features on the segment
         for feature_stranded in features_on_seg: # go through all the segment's features, and print the gff line for each
-            line = create_line_graph_gff(feature_stranded,segment_id)
+            line = create_line_graph_gff(feature_stranded,segment_stranded)
             write_line(line,output_graph_gff,False)
     write_line("",output_graph_gff,True)