#NEWLINE# def combcall2draw( cvsfile , write_func , * args ) : #NEWLINE# #INDENT# 'User specifies several columns to consider, this tool will call\n    regions where either of the column is above its threshold.\n    \n    ' #NEWLINE# argv = args [ 0 ] #NEWLINE# if ( len ( argv ) < 6 ) : #NEWLINE# #INDENT# 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" ) #NEWLINE# sys . exit ( ) #NEWLINE# #UNINDENT# cor_column = cvsfile . fieldnames [ int ( argv [ 0 ] ) ] #NEWLINE# var_columns = map ( ( lambda x : cvsfile . fieldnames [ int ( x ) ] ) , argv [ 1 ] . split ( ',' ) ) #NEWLINE# cutoffs = map ( float , argv [ 2 ] . split ( ',' ) ) #NEWLINE# min_len = int ( argv [ 3 ] ) #NEWLINE# max_gap = int ( argv [ 4 ] ) #NEWLINE# wtrack = WigTrackI ( ) #NEWLINE# add_func = wtrack . add_loc #NEWLINE# for l in cvsfile : #NEWLINE# #INDENT# cor = l . setdefault ( cor_column , None ) #NEWLINE# if ( ( not cor ) or ( cor == 'NA' ) ) : #NEWLINE# #INDENT# continue #NEWLINE# #UNINDENT# for i in range ( len ( var_columns ) ) : #NEWLINE# #INDENT# var_column = var_columns [ i ] #NEWLINE# cutoff = cvsfile [ i ] #NEWLINE# var = l . setdefault ( var_column , None ) #NEWLINE# if ( var and ( var != 'NA' ) and ( float ( var ) > cutoff ) ) : #NEWLINE# #INDENT# ( chrom , start , end ) = cor . split ( '.' ) #NEWLINE# add_func ( chrom , int ( start ) , 1.1 ) #NEWLINE# break #NEWLINE# #UNINDENT# #UNINDENT# #UNINDENT# wtrack . span = ( int ( end ) - int ( start ) ) #NEWLINE# bpeaks = wtrack . call_peaks ( cutoff = 1.0 , min_length = min_len , max_gap = max_gap ) #NEWLINE# fhd = open ( argv [ 5 ] . replace ( 'pdf' , 'bed' ) , 'w' ) #NEWLINE# fhd . write ( bpeaks . tobed ( ) ) #NEWLINE# from Bio.Graphics import BasicChromosome #NEWLINE# from reportlab.lib.colors import gray , black , white #NEWLINE# entries = [ ( 'chrI' , 15072419 ) , ( 'chrII' , 15279316 ) , ( 'chrIII' , 13783681 ) , ( 'chrIV' , 17493784 ) , ( 'chrV' , 20919398 ) , ( 'chrX' , 17718852 ) ] #NEWLINE# max_length = max ( [ x [ 1 ] for x in entries ] ) #NEWLINE# chr_diagram = BasicChromosome . Organism ( ) #NEWLINE# for ( name , length ) in entries : #NEWLINE# #INDENT# cur_chromosome = BasicChromosome . Chromosome ( name ) #NEWLINE# cur_chromosome . scale_num = ( max_length * 1.1 ) #NEWLINE# start = BasicChromosome . TelomereSegment ( ) #NEWLINE# start . scale = ( 0.05 * max_length ) #NEWLINE# start . fill_color = gray #NEWLINE# cur_chromosome . add ( start ) #NEWLINE# try : #NEWLINE# #INDENT# cpeaks = bpeaks . peaks [ name ] #NEWLINE# #UNINDENT# except : #NEWLINE# #INDENT# cpeaks = [ ] #NEWLINE# #UNINDENT# body_regions = [ ] #NEWLINE# last_pos = 0 #NEWLINE# for p in cpeaks : #NEWLINE# #INDENT# body_regions . append ( ( ( p [ 0 ] - last_pos ) , white ) ) #NEWLINE# body_regions . append ( ( ( p [ 1 ] - p [ 0 ] ) , black ) ) #NEWLINE# last_pos = p [ 1 ] #NEWLINE# assert ( p [ 1 ] < length ) #NEWLINE# #UNINDENT# body_regions . append ( ( ( length - last_pos ) , white ) ) #NEWLINE# for ( b , c ) in body_regions : #NEWLINE# #INDENT# body = BasicChromosome . ChromosomeSegment ( ) #NEWLINE# body . fill_color = c #NEWLINE# body . scale = b #NEWLINE# cur_chromosome . add ( body ) #NEWLINE# #UNINDENT# end = BasicChromosome . TelomereSegment ( inverted = True ) #NEWLINE# end . scale = ( 0.05 * max_length ) #NEWLINE# end . fill_color = gray #NEWLINE# cur_chromosome . add ( end ) #NEWLINE# chr_diagram . add ( cur_chromosome ) #NEWLINE# #UNINDENT# chr_diagram . draw ( argv [ 5 ] , 'Highlight regions in Caenorhabditis elegans' )