#NEWLINE# def main( ) : #NEWLINE# #INDENT# usage = 'usage: %prog [options]' #NEWLINE# description = 'Call peaks from Wiggle file.' #NEWLINE# optparser = OptionParser ( version = '%prog 0.1' , description = description , usage = usage , add_help_option = False ) #NEWLINE# optparser . add_option ( '-h' , '--help' , action = 'help' , help = 'Show this help message and exit.' ) #NEWLINE# optparser . add_option ( '-i' , '--wfile' , dest = 'wfile' , type = 'string' , help = 'input wiggle file. *REQUIRED' ) #NEWLINE# optparser . add_option ( '-o' , '--bfile' , dest = 'bfile' , type = 'string' , help = 'output bed file. *REQUIRED' ) #NEWLINE# optparser . add_option ( '-c' , '--cutoff' , dest = 'cutoff' , type = 'float' , help = 'pvalue cutoff. default: 1e-5' , default = 1e-05 ) #NEWLINE# optparser . add_option ( '-l' , '--min-length' , dest = 'minlen' , type = 'int' , help = 'minimum length of peak, default: 300' , default = 300 ) #NEWLINE# optparser . add_option ( '-g' , '--maxgap' , dest = 'maxgap' , type = 'int' , help = 'maximum gap between significant points in a peak, default: 50' , default = 50 ) #NEWLINE# optparser . add_option ( '--skip-normal' , dest = 'normalize' , action = 'store_false' , help = 'Skip normalization the data before calling peaks' , default = True ) #NEWLINE# optparser . add_option ( '--skip-null' , dest = 'nullmodel' , action = 'store_false' , help = 'Skip sampling the data to use null (background) data to calculate mean and std' , default = True ) #NEWLINE# optparser . add_option ( '--sample-percent' , dest = 'samplepercent' , type = 'int' , default = 10 , help = 'Percentage to sample the data to build null model. Default: 10' ) #NEWLINE# ( options , args ) = optparser . parse_args ( ) #NEWLINE# if ( ( not options . wfile ) or ( not options . bfile ) or ( not options . cutoff ) ) : #NEWLINE# #INDENT# optparser . print_help ( ) #NEWLINE# sys . exit ( ) #NEWLINE# #UNINDENT# f = options . wfile #NEWLINE# if ( not os . path . isfile ( f ) ) : #NEWLINE# #INDENT# error ( ( '%s is not valid!' % f ) ) #NEWLINE# sys . exit ( 1 ) #NEWLINE# #UNINDENT# try : #NEWLINE# #INDENT# fhd = open ( f ) #NEWLINE# #UNINDENT# except : #NEWLINE# #INDENT# error ( ( "Can't read %s" % f ) ) #NEWLINE# sys . exit ( 1 ) #NEWLINE# #UNINDENT# try : #NEWLINE# #INDENT# bfhd = open ( options . bfile , 'w' ) #NEWLINE# #UNINDENT# except : #NEWLINE# #INDENT# error ( ( "Can't open %s to write" % options . bfile ) ) #NEWLINE# sys . exit ( 1 ) #NEWLINE# #UNINDENT# info ( 'open wiggle file...' ) #NEWLINE# wio = WiggleIO . WiggleIO ( fhd ) #NEWLINE# info ( 'construct wiggle track object...' ) #NEWLINE# wtrack = wio . build_wigtrack ( ) #NEWLINE# if options . normalize : #NEWLINE# #INDENT# info ( 'Normalizing...' ) #NEWLINE# ( sum_v , max_v , min_v , mean_v , std_v ) = wtrack . normalize ( null = options . nullmodel , sample_percent = options . samplepercent ) #NEWLINE# if options . nullmodel : #NEWLINE# #INDENT# info ( 'before normalization (NULL model):' ) #NEWLINE# info ( ( 'mean of data: %.2f' % mean_v ) ) #NEWLINE# info ( ( 'std of data: %.2f' % std_v ) ) #NEWLINE# info ( 'after normalization (NULL model):' ) #NEWLINE# #UNINDENT# else : #NEWLINE# #INDENT# info ( 'before normalization:' ) #NEWLINE# info ( ( 'mean of data: %.2f' % mean_v ) ) #NEWLINE# info ( ( 'std of data: %.2f' % std_v ) ) #NEWLINE# info ( 'after normalization' ) #NEWLINE# #UNINDENT# #UNINDENT# if options . nullmodel : #NEWLINE# #INDENT# ( sum_v , max_v , min_v , mean_v , std_v ) = wtrack . null_model_summary ( sample = options . samplepercent ) #NEWLINE# #UNINDENT# else : #NEWLINE# #INDENT# ( sum_v , max_v , min_v , mean_v , std_v ) = wtrack . summary ( ) #NEWLINE# #UNINDENT# info ( ( 'mean of data: %.2f' % mean_v ) ) #NEWLINE# info ( ( 'std of data: %.2f' % std_v ) ) #NEWLINE# info ( 'Calculate score cutoff using pvalue cutoff...' ) #NEWLINE# scorecutoff = normal_cdf_inv ( options . cutoff , mu = mean_v , sigma2 = std_v , lower = False ) #NEWLINE# info ( ( 'scorecutoff is %f' % scorecutoff ) ) #NEWLINE# info ( 'call peaks...' ) #NEWLINE# wpeaks = wtrack . call_peaks ( cutoff = scorecutoff , min_length = options . minlen , max_gap = options . maxgap ) #NEWLINE# info ( 'write to bed file...' ) #NEWLINE# bfhd . write ( wpeaks . tobed ( ) ) #NEWLINE# bfhd . close ( ) #NEWLINE# info ( 'finished' )