@transform( getLongIntervalGeneGTF, suffix(".gtf.gz"), ".H3K27Me3.log" )    
def longGeneChromatinIntersection(infile, outfile):
    chromatin = P.asList(PARAMS["bed_h3k27me3"])
    track = P.snip(os.path.basename(infile), ".gtf.gz")
    for bed in chromatin:
        chromatin_track = P.snip(os.path.basename(bed), ".bed")
        statement = '''zcat %(infile)s | python %(scriptsdir)s/gff2bed.py --is-gtf > long_intervals/%(track)s.bed;
                       cat long_intervals/%(track)s.bed %(bed)s | awk 'OFS="\\t" {print $1,$2,$3}'
                       | mergeBed -i stdin | awk 'OFS="\\t" {print $1,$2,$3,"merged"NR}' > long_intervals/%(track)s_%(chromatin_track)s.merged.bed;
                       echo "Track" > long_intervals/%(track)s_%(chromatin_track)s.counts; 
                       echo "%(track)s" >> long_intervals/%(track)s_%(chromatin_track)s.counts; 
                       echo "Chromatin_track" >> long_intervals/%(track)s_%(chromatin_track)s.counts; 
                       echo "%(chromatin_track)s" >> long_intervals/%(track)s_%(chromatin_track)s.counts; 
                       echo "Total_merged_intervals" >> long_intervals/%(track)s_%(chromatin_track)s.counts; 
                       cat long_intervals/%(track)s_%(chromatin_track)s.merged.bed | wc -l >> long_intervals/%(track)s_%(chromatin_track)s.counts; 
                       echo "track_and_chromatin_track" >> long_intervals/%(track)s_%(chromatin_track)s.counts; 
                       intersectBed -a long_intervals/%(track)s_%(chromatin_track)s.merged.bed -b long_intervals/%(track)s.bed -u | intersectBed -a stdin -b %(bed)s -u > long_intervals/%(track)s_%(chromatin_track)s.shared.bed; 
                       cat long_intervals/%(track)s_%(chromatin_track)s.shared.bed | wc -l >> long_intervals/%(track)s_%(chromatin_track)s.counts; 
                       echo "chromatin_track_only" >> long_intervals/%(track)s_%(chromatin_track)s.counts; 
                       intersectBed -a long_intervals/%(track)s_%(chromatin_track)s.merged.bed -b long_intervals/%(track)s.bed -v > long_intervals/%(track)s_%(chromatin_track)s.%(chromatin_track)s.unique.bed; 
                       cat long_intervals/%(track)s_%(chromatin_track)s.%(chromatin_track)s.unique.bed | wc -l >> long_intervals/%(track)s_%(chromatin_track)s.counts; 
                       echo "track_only" >> long_intervals/%(track)s_%(chromatin_track)s.counts; 
                       intersectBed -a long_intervals/%(track)s_%(chromatin_track)s.merged.bed -b %(bed)s -v > long_intervals/%(track)s_%(chromatin_track)s.%(track)s.unique.bed; 
                       cat long_intervals/%(track)s_%(chromatin_track)s.%(track)s.unique.bed | wc -l >> long_intervals/%(track)s_%(chromatin_track)s.counts;                   
                       sed -i '{N;s/\\n/\\t/g}' long_intervals/%(track)s_%(chromatin_track)s.counts;
                       touch %(outfile)s '''
        P.run()

############################################################
calculate intersection of chromatin marks and CAPSeq overlapped genes
dataset/ETHPy150Open CGATOxford/cgat/obsolete/pipeline_proj007.py/longGeneChromatinIntersection