MACS2/bdgbroadcall_cmd.py

   1 # Time-stamp: <2019-09-25 10:04:48 taoliu>
   2
   3 """Description: Fine-tuning script to call broad peaks from a single
   4 bedGraph track for 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 # ------------------------------------
  13 # python modules
  14 # ------------------------------------
  15
  16 import sys
  17 import os
  18 import logging
  19 from MACS2.IO import BedGraphIO
  20 # ------------------------------------
  21 # constants
  22 # ------------------------------------
  23 logging.basicConfig(level=20,
  24                     format='%(levelname)-5s @ %(asctime)s: %(message)s ',
  25                     datefmt='%a, %d %b %Y %H:%M:%S',
  26                     stream=sys.stderr,
  27                     filemode="w"
  28                     )
  29
  30 # ------------------------------------
  31 # Misc functions
  32 # ------------------------------------
  33 error   = logging.critical              # function alias
  34 warn    = logging.warning
  35 debug   = logging.debug
  36 info    = logging.info
  37 # ------------------------------------
  38 # Classes
  39 # ------------------------------------
  40
  41 # ------------------------------------
  42 # Main function
  43 # ------------------------------------
  44 def run( options ):
  45     info("Read and build bedGraph...")
  46     bio = BedGraphIO.bedGraphIO(options.ifile)
  47     btrack = bio.build_bdgtrack(baseline_value=0)
  48
  49     info("Call peaks from bedGraph...")
  50
  51     bpeaks = btrack.call_broadpeaks (lvl1_cutoff=options.cutoffpeak, lvl2_cutoff=options.cutofflink, min_length=options.minlen, lvl1_max_gap=options.lvl1maxgap, lvl2_max_gap=options.lvl2maxgap)
  52
  53     info("Write peaks...")
  54
  55     if options.ofile:
  56         bf = open( os.path.join( options.outdir, options.ofile ), "w" )
  57         options.oprefix = options.ofile
  58     else:
  59         bf = open ( os.path.join( options.outdir, "%s_c%.1f_C%.2f_l%d_g%d_G%d_broad.bed12" % (options.oprefix,options.cutoffpeak,options.cutofflink,options.minlen,options.lvl1maxgap,options.lvl2maxgap)), "w" )
  60     bpeaks.write_to_gappedPeak(bf, name_prefix=(options.oprefix+"_broadRegion").encode(), score_column="score", trackline=options.trackline)
  61     info("Done")