Sort our progress

def annotated_reads(txome,sorted_read_stream,output,args):
   if not os.path.exists(args.tempdir+'/'+output):
      os.makedirs(args.tempdir+'/'+output)
   ts = OrderedStream(iter(txome.transcript_stream()))
   rs = OrderedStream(GPDStream(sorted_read_stream))
   mls = MultiLocusStream([ts,rs])
   
   aof = gzip.open(args.tempdir+'/'+output+'/annotated.txt.gz','w')
   of1 = gzip.open(args.tempdir+'/'+output+'/unannotated_multiexon.gpd.gz','w')
   of2 = gzip.open(args.tempdir+'/'+output+'/partial_annotated_multiexon.gpd.gz','w')
   of3 = gzip.open(args.tempdir+'/'+output+'/unannotated_singleexon.gpd.gz','w')
   z = 0
   for ml in mls:
      refs, reads = ml.payload
      if len(refs) == 0: continue
      """Check and see if we have single exons annotated"""
      single_exon_refs = [x for x in refs if x.get_exon_count()==1]
      single_exon_reads = [x for x in reads if x.get_exon_count()==1]
      #print str(len(single_exon_refs))+" "+str(len(single_exon_reads))
      unannotated_single_exon_reads = []
      for seread in single_exon_reads:
         ovs  = [(x.exons[0].length,seread.exons[0].length,x.exons[0].overlap_size(seread.exons[0]),x) for x in single_exon_refs if x.exons[0].overlaps(seread.exons[0])]
         #check for minimum overlap
         ovs = [x for x in ovs if x[2] >= args.single_exon_minimum_overlap]
         ovs = [x for x in ovs if float(x[2])/float(max(x[0],x[1])) >= args.single_exon_mutual_overlap]
         ovs = sorted(ovs,key=lambda x: float(x[2])/float(max(x[0],x[1])),reverse=True)
         if len(ovs) == 0: 
            unannotated_single_exon_reads.append(seread)
            continue
         best_ref = ovs[0][3]
         aof.write(seread.name+"\t"+best_ref.name+"\tSE"+"\n")

      """May want an optional check for any better matches among exons"""
      #now we can look for best matches among multi-exon transcripts
      reads = [x for x in reads if x.get_exon_count() > 1]
      multiexon_refs = [x for x in refs if x.get_exon_count() > 1]
      unannotated_multi_exon_reads = []
      partial_annotated_multi_exon_reads = []
      for read in reads:
         # we dont' need to have multiple exons matched. one is enough to call
         ovs = [y for y in [(x,x.exon_overlap(read,
                                              multi_minover=args.multi_exon_minimum_overlap,
                                              multi_endfrac=args.multi_exon_end_frac,
                                              multi_midfrac=args.multi_exon_mid_frac,
                                              multi_consec=False)) for x in multiexon_refs] if y[1]]
         for o in ovs: o[1].analyze_overs()
         full =  sorted([x for x in ovs if x[1].is_subset()==1],
                        key = lambda y: float(y[1].overlap_size())/float(max(y[1].tx_obj1.length,y[1].tx_obj2.length)),
                        reverse=True
                       )
         if len(full) > 0:
            aof.write(read.name+"\t"+full[0][0].name+"\tfull"+"\n")
            continue
         subset =  sorted([x for x in ovs if x[1].is_subset()==2],
                          key = lambda y: (y[1].match_exon_count(),
                                           y[1].min_overlap_fraction()),
                          reverse = True
                         )
         if len(subset) > 0:
            aof.write(read.name+"\t"+subset[0][0].name+"\tpartial"+"\n")
            continue
         #check for supersets
         superset = sorted([x for x in ovs if x[1].is_subset()==3],
                          key = lambda y: (y[1].match_exon_count(),
                                           y[1].min_overlap_fraction()),
                          reverse = True
                         )
         if len(superset) > 0:
            partial_annotated_multi_exon_reads.append((read,superset[0][0]))
            #print read.name+"\t"+superset[0][0].name+"\tsuper"
            continue
         #check for noncompatible overlaps
         overset = sorted([x for x in ovs if x[1].match_exon_count > 0],
                          key = lambda y: (y[1].consecutive_exon_count(),
                                           y[1].min_overlap_fraction()),
                          reverse = True
                         )
         #print [(x[1].consecutive_exon_count(), x[1].min_overlap_fraction()) for x in overset]
         if len(overset) > 0:
            partial_annotated_multi_exon_reads.append((read,overset[0][0]))
            #print read.name+"\t"+overset[0][0].name+"\tover"
            continue
         unannotated_multi_exon_reads.append(read)
      """Now we have partially annotated multi and unannotated multi and unannotated single"""
      if len(unannotated_multi_exon_reads) > 0:
         sys.stderr.write(str(z)+" "+str(len(unannotated_multi_exon_reads))+"   \r")
         d = Deconvolution(downsample(unannotated_multi_exon_reads,args.downsample_locus))
         groups = d.parse(tolerance=20,downsample=args.downsample)
         for tx in groups:
            z+=1
            of1.write(tx.get_gpd_line()+"\n")
      if len(partial_annotated_multi_exon_reads) > 0:
         sys.stderr.write(str(z)+" "+str(len(partial_annotated_multi_exon_reads))+"   \r")
         ### set the direction of the transcript
         for v in partial_annotated_multi_exon_reads:
            v[0].set_strand(v[1].direction)
         ### set the gene name of the transcript
         for v in partial_annotated_multi_exon_reads:
            v[0].set_gene_name(v[1].gene_name)         
         d = Deconvolution(downsample([x[0] for x in partial_annotated_multi_exon_reads],args.downsample_locus))
         groups = d.parse(tolerance=20,downsample=args.downsample,use_gene_names=True)
         for tx in groups:
            z += 1
            of2.write(tx.get_gpd_line()+"\n")
      """Do the unannotated single exon reads"""
      g = Graph()
      for r in unannotated_single_exon_reads:
         if len([x for x in partial_annotated_multi_exon_reads if x[0].overlaps(r)]) > 0: continue
         if len([x for x in unannotated_multi_exon_reads if x.overlaps(r)]) > 0: continue
         if len([x for x in txome.transcripts if x.overlaps(r)]) > 0: continue
         g.add_node(Node(r))
      for i in range(0,len(g.nodes)):
         for j in range(0,len(g.nodes)):
            if i == j: continue
            if g.nodes[i].payload.overlaps(g.nodes[j].payload):
               g.add_edge(Edge(g.nodes[i],g.nodes[j]))
      g.merge_cycles()
      for r in g.roots:
         se = []
         se += r.payload_list
         children = g.get_children(r)
         for child in children:  se += child.payload_list
         rng = GenomicRange(se[0].exons[0].chr,
                            min([x.exons[0].start for x in se]),
                            max([x.exons[0].end for x in se]))
         tx = Transcript([rng],Transcript.Options(direction='+'))
         of3.write(tx.get_gpd_line()+"\n")

   of1.close()
   of2.close()
   of3.close()
   sys.stderr.write("\n")      
   """Sort our progress"""
   sys.stderr.write("sort transcriptome\n")
   sort_gpd(args.tempdir+'/'+output+'/partial_annotated_multiexon.gpd.gz',args.tempdir+'/'+output+'/partial_annotated_multiexon.sorted.gpd.gz',args)
   sort_gpd(args.tempdir+'/'+output+'/unannotated_multiexon.gpd.gz',args.tempdir+'/'+output+'/unannotated_multiexon.sorted.gpd.gz',args)
   sort_gpd(args.tempdir+'/'+output+'/unannotated_singleexon.gpd.gz',args.tempdir+'/'+output+'/unannotated_singleexon.sorted.gpd.gz',args)
   """We still have the unannotated single exon reads to deal with"""