@files( [( ["%s%s.bed" % (x,y) for y in REPLICATES], x % "%s.reproducibility", x) for x in TRACKS_EXPERIMENTS ] )
    def makeReproducibility( infile, outfile, track ):
        '''compute overlap between intervals.

        Note the exon percentages are approxmations assuming that there are
        not more than one intervals in one set overlapping one in the other set.
        '''

        dbhandle = sqlite3.connect( PARAMS["database"] )

        data = []
        for replicate in REPLICATES:
            cc = dbhandle.cursor()
            statement = "SELECT contig, start, end, peakval FROM %(track)s%(replicate)s_intervals" % locals()
            cc.execute( statement )
            data.append( cc.fetchall() )

        ma = int(max( [ x[3] for x in itertools.chain( *data ) ] ) + 1)

        nexons1, nexons2, nexons_overlapping = [0] * ma, [0] * ma, [0] * ma
        nbases1, nbases2, nbases_overlapping = [0] * ma, [0] * ma, [0] * ma

        idx = IndexedGenome.IndexedGenome()
        for contig, start, end, peakval in data[0]:
            peakval = int(peakval)
            for x in range( 0, peakval):
                nexons1[x] += 1
                nbases1[x] += end - start
            idx.add( contig, start, end, peakval )

        for contig, start, end, peakval in data[1]:
            peakval = int(peakval)
            for x in range( 0, peakval):
                nexons2[x] += 1
                nbases2[x] += end - start

            try:
                intervals = list(idx.get( contig, start, end ))
            except KeyError:
                continue

            if len(intervals) == 0: continue

            for other_start,other_end,other_peakval in intervals:
                ovl = min(end, other_end ) - max(start, other_start )
                for x in range( 0, min( other_peakval, peakval) ):
                    nexons_overlapping[x] += 1
                    nbases_overlapping[x] += ovl

        outs = open( outfile, "w" )
        outs.write("peakval\tnexons1\tnexons2\tnexons_union\tnexons_ovl\tpexons_ovl\tpexons_union\tnbases1\tnbases2\tnbases_union\tnbases_ovl\tpbases_ovl\tpbases_union\n" )
        total_bases = nbases1[0] + nbases2[0] - nbases_overlapping[0]
        for x in range(0, ma):
            bases_union = nbases1[x] + nbases2[x] - nbases_overlapping[x]
            exons_union = nexons1[x] + nexons2[x] - nexons_overlapping[x]
            outs.write( "\t".join( (str(x),
                                    str(nexons1[x]),
                                    str(nexons2[x]),
                                    str(exons_union),
                                    str(nexons_overlapping[x]),
                                    IOTools.prettyPercent( nbases_overlapping[x], bases_union),
                                    IOTools.prettyPercent( bases_union, total_bases),
                                    str(nbases1[x]),
                                    str(nbases2[x]),
                                    str(bases_union),
                                    str(nbases_overlapping[x]),
                                    IOTools.prettyPercent( nbases_overlapping[x], bases_union),
                                    IOTools.prettyPercent( bases_union, total_bases),
                                    )) + "\n" )

        outs.close()
dataset/ETHPy150Open CGATOxford/cgat/obsolete/pipeline_vitaminD.py/makeReproducibility/swapped operands of `%`