1# Time-stamp: <2019-09-25 10:04:08 taoliu> 2 3"""Description: Naive call peaks from a single bedGraph track for 4scores. 5 6This code is free software; you can redistribute it and/or modify it 7under the terms of the BSD License (see the file LICENSE included with 8the distribution). 9""" 10 11# ------------------------------------ 12# python modules 13# ------------------------------------ 14import sys 15import os 16import logging 17from MACS2.IO import BedGraphIO 18# ------------------------------------ 19# constants 20# ------------------------------------ 21logging.basicConfig(level=20, 22 format='%(levelname)-5s @ %(asctime)s: %(message)s ', 23 datefmt='%a, %d %b %Y %H:%M:%S', 24 stream=sys.stderr, 25 filemode="w" 26 ) 27 28# ------------------------------------ 29# Misc functions 30# ------------------------------------ 31error = logging.critical # function alias 32warn = logging.warning 33debug = logging.debug 34info = logging.info 35# ------------------------------------ 36# Classes 37# ------------------------------------ 38 39# ------------------------------------ 40# Main function 41# ------------------------------------ 42def run( options ): 43 info("Read and build bedGraph...") 44 bio = BedGraphIO.bedGraphIO(options.ifile) 45 btrack = bio.build_bdgtrack(baseline_value=0) 46 47 if options.cutoff_analysis: 48 info("Analyze cutoff vs number of peaks/total length of peaks/average length of peak") 49 cutoff_analysis_result = btrack.cutoff_analysis( int(options.maxgap), int(options.minlen), 50 ) 50 info("Write report...") 51 if options.ofile: 52 fhd = open( os.path.join( options.outdir, options.ofile ), 'w' ) 53 else: 54 fhd = open ( os.path.join( options.outdir, "%s_l%d_g%d_cutoff_analysis.txt" % (options.oprefix,options.minlen,options.maxgap)), "w" ) 55 fhd.write( cutoff_analysis_result ) 56 info("Done") 57 else: 58 info("Call peaks from bedGraph...") 59 peaks = btrack.call_peaks(cutoff=float(options.cutoff),min_length=int(options.minlen),max_gap=int(options.maxgap),call_summits=options.call_summits) 60 61 info("Write peaks...") 62 if options.ofile: 63 options.oprefix = options.ofile 64 nf = open( os.path.join( options.outdir, options.ofile ), 'w' ) 65 else: 66 nf = open ( os.path.join( options.outdir, "%s_c%.1f_l%d_g%d_peaks.narrowPeak" % (options.oprefix,options.cutoff,options.minlen,options.maxgap)), "w" ) 67 peaks.write_to_narrowPeak(nf, name=options.oprefix.encode(), name_prefix=(options.oprefix+"_narrowPeak").encode(), score_column="score", trackline=options.trackline) 68 info("Done") 69 70 71 72