def combcall2draw ( cvsfile, write_func, *args ):
    argv = args[0]
    if len(argv) < 6:
        sys.stderr.write("Need 6 extra arguments for 'combcall2draw', options <loc column> <score column1[,score column2,...]> <cutoff1[,cutoff2,cutoff3]> <min length> <max gap> <pdf filename>\ne.g. command: <0> <1,2,3> <0.5,0.6,0.7> <10000> <2000> <a.pdf>, means to use the first column as genome coordinations to call enriched regions from the combinition of #1, #2 and #3, the thresholds to call enriched region are 0.5 for column 1, 0.6 for column 2 and 0.7 for column 3, the minimum length of region is 10k, and the maximum gap to link two nearby regions is 2k. Then the figure will be saved in a.pdf.\n")
        sys.exit()        
    cor_column = cvsfile.fieldnames[int(argv[0])]
    var_columns = map(lambda x:cvsfile.fieldnames[int(x)],argv[1].split(","))
    cutoffs = map(float,argv[2].split(","))
    
    min_len = int(argv[3])
    max_gap = int(argv[4])
    wtrack = WigTrackI()                 # combined track containing 1 if either of track is above cutoff
    add_func = wtrack.add_loc
    
    for l in cvsfile:
        cor = l.setdefault(cor_column,None)
        if not cor or cor =="NA":
            continue
        
        for i in range(len(var_columns)):
            var_column = var_columns[i]
            cutoff = cutoffs[i]
            var = l.setdefault(var_column,None)
            if var and var != "NA" and float(var) > cutoff:
                (chrom,start,end) = cor.split(".")
                add_func(chrom,int(start),1.1)
                break
            
    wtrack.span = int(end)-int(start)
    bpeaks = wtrack.call_peaks(cutoff=1.0,min_length=min_len,max_gap=max_gap)
    #f = argv[5]
    fhd = open(argv[5].replace("pdf","bed"),"w")
    fhd.write(bpeaks.tobed())
    
    from Bio.Graphics import BasicChromosome
    from reportlab.lib.colors import gray, black, white
    entries = [("chrI", 15072419),
               ("chrII", 15279316),
               ("chrIII", 13783681),
               ("chrIV", 17493784),
               ("chrV", 20919398),
               ("chrX", 17718852)]
    max_length = max([x[1] for x in entries])
    chr_diagram = BasicChromosome.Organism()
    for name, length in entries:
        cur_chromosome = BasicChromosome.Chromosome(name)
        #Set the length, adding and extra 20 percent for the tolomeres:
        cur_chromosome.scale_num = max_length * 1.1
        # Add an opening telomere
        start = BasicChromosome.TelomereSegment()
        start.scale = 0.05 * max_length
        start.fill_color=gray
        cur_chromosome.add(start)
        #Add a body - using bp as the scale length here.
        try:
            cpeaks = bpeaks.peaks[name]
        except:
            cpeaks = []
        body_regions = []
        last_pos = 0
        for p in cpeaks:
            body_regions.append( (p[0]-last_pos,white) ) # outside regions
            body_regions.append( (p[1]-p[0],black) ) # enriched regions
            last_pos = p[1]
            assert p[1] < length
        body_regions.append( (length-last_pos,white) ) # last part

        for b,c in body_regions:
            body = BasicChromosome.ChromosomeSegment()
            body.fill_color= c
            body.scale = b
            cur_chromosome.add(body)
        
        #Add a closing telomere
        end = BasicChromosome.TelomereSegment(inverted=True)
        end.scale = 0.05 * max_length
        end.fill_color=gray
        cur_chromosome.add(end)
        #This chromosome is done
        chr_diagram.add(cur_chromosome)
        
    chr_diagram.draw(argv[5], "Highlight regions in Caenorhabditis elegans" )
User specifies several columns to consider, this tool will call
dataset/ETHPy150Open taoliu/taolib/Scripts/ce_histone_matrix.py/combcall2draw