MACS2/bdgpeakcall_cmd.py

   1 # Time-stamp: <2019-09-25 10:04:08 taoliu>
   2
   3 """Description: Naive call peaks from a single bedGraph track for
   4 scores.
   5
   6 This code is free software; you can redistribute it and/or modify it
   7 under the terms of the BSD License (see the file LICENSE included with
   8 the distribution).
   9 """
  10
  11 # ------------------------------------
  12 # python modules
  13 # ------------------------------------
  14 import sys
  15 import os
  16 import logging
  17 from MACS2.IO import BedGraphIO
  18 # ------------------------------------
  19 # constants
  20 # ------------------------------------
  21 logging.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 # ------------------------------------
  31 error   = logging.critical              # function alias
  32 warn    = logging.warning
  33 debug   = logging.debug
  34 info    = logging.info
  35 # ------------------------------------
  36 # Classes
  37 # ------------------------------------
  38
  39 # ------------------------------------
  40 # Main function
  41 # ------------------------------------
  42 def 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