@transform( thresholdFoldChange, suffix(".foldchange"), ".shared.foldchange")
def sharedIntervalsFoldChangeThreshold(infile, outfile):
    '''identify shared intervals between datasets at different foldchange thresholds'''

    # Open foldchange file and read
    fc_file = open(infile, "r")
    foldchange = []
    for line in fc_file:
        threshold, interval_count = line.split()
        foldchange.append(threshold)

    in_track = P.snip( os.path.basename( infile ), ".foldchange" )
    in_dir = os.path.dirname(infile)
    out_dir = in_dir.replace("macs","fc")
    try: os.mkdir( out_dir )
    except __HOLE__: pass

    if os.path.exists(outfile):
        statement = '''rm %(outfile)s'''
        P.run()

    # for each foldchange 
    for fc in foldchange:

        in_bed = "foldchange/" + in_track.replace("-","_") + ".fc" + str(fc) + ".bed"

        # For each track
        for track in TRACKS:
            if (str(track) != in_track):
                compare_bed = "foldchange/" + str(track).replace("-","_") + ".fc" + str(fc) + ".bed"
                statement = '''echo %(track)s %(fc)s >> %(outfile)s; intersectBed -a %(in_bed)s -b %(compare_bed)s -u | wc -l >> %(outfile)s; ''' 
                P.run()

    statement = '''sed -i '{N;s/\\n/\\t/}' %(outfile)s; sed -i '{s/ /\\t/g}' %(outfile)s; '''
    P.run()

############################################################
dataset/ETHPy150Open CGATOxford/cgat/obsolete/pipeline_proj012_chipseq.py/sharedIntervalsFoldChangeThreshold