MACS2/bdgcmp_cmd.py

   1 # Time-stamp: <2019-09-25 10:04:37 taoliu>
   2
   3 """Description: compare bdg files
   4
   5 This code is free software; you can redistribute it and/or modify it
   6 under the terms of the BSD License (see the file LICENSE included with
   7 the distribution).
   8 """
   9
  10 import sys
  11 import os
  12 import logging
  13
  14 from MACS2.IO import BedGraphIO
  15 from MACS2.OptValidator import opt_validate_bdgcmp as opt_validate
  16
  17 from math import log as mlog
  18
  19 # ------------------------------------
  20 # constants
  21 # ------------------------------------
  22 logging.basicConfig(level=20,
  23                     format='%(levelname)-5s @ %(asctime)s: %(message)s ',
  24                     datefmt='%a, %d %b %Y %H:%M:%S',
  25                     stream=sys.stderr,
  26                     filemode="w"
  27                     )
  28
  29 # ------------------------------------
  30 # Misc functions
  31 # ------------------------------------
  32 error   = logging.critical              # function alias
  33 warn    = logging.warning
  34 debug   = logging.debug
  35 info    = logging.info
  36 # ------------------------------------
  37 # Main function
  38 # ------------------------------------
  39
  40 def run( options ):
  41     options = opt_validate( options )
  42     scaling_factor = options.sfactor
  43     pseudo_depth = 1.0/scaling_factor   # not an actual depth, but its reciprocal, a trick to override SPMR while necessary.
  44
  45     info("Read and build treatment bedGraph...")
  46     tbio = BedGraphIO.bedGraphIO(options.tfile)
  47     tbtrack = tbio.build_bdgtrack()
  48
  49     info("Read and build control bedGraph...")
  50     cbio = BedGraphIO.bedGraphIO(options.cfile)
  51     cbtrack = cbio.build_bdgtrack()
  52
  53     info("Build scoreTrackII...")
  54     sbtrack = tbtrack.make_scoreTrackII_for_macs( cbtrack, depth1 = pseudo_depth, depth2 = pseudo_depth )
  55     if abs(scaling_factor-1) > 1e-6:
  56         # Only for the case while your input is SPMR from MACS2 callpeak; Let's override SPMR.
  57         info("Values in your input bedGraph files will be multiplied by %f ..." % scaling_factor)
  58         sbtrack.change_normalization_method( ord('M') ) # a hack to override SPMR
  59     sbtrack.set_pseudocount( options.pseudocount )
  60
  61     already_processed_method_list = []
  62     for (i, method) in enumerate(options.method):
  63         if method in already_processed_method_list:
  64             continue
  65         else:
  66             already_processed_method_list.append( method )
  67
  68         info("Calculate scores comparing treatment and control by '%s'..." % method)
  69         if options.ofile:
  70             ofile = os.path.join( options.outdir, options.ofile[ i ] )
  71         else:
  72             ofile = os.path.join( options.outdir, options.oprefix + "_" + method + ".bdg" )
  73         # build score track
  74         if method == 'ppois':
  75             sbtrack.change_score_method( ord('p') )
  76         elif method == 'qpois':
  77             sbtrack.change_score_method( ord('q') )
  78         elif method == 'subtract':
  79             sbtrack.change_score_method( ord('d') )
  80         elif method == 'logFE':
  81             sbtrack.change_score_method( ord('f') )
  82         elif method == 'FE':
  83             sbtrack.change_score_method( ord('F') )
  84         elif method == 'logLR':             # log likelihood
  85             sbtrack.change_score_method( ord('l') )
  86         elif method == 'slogLR':             # log likelihood
  87             sbtrack.change_score_method( ord('s') )
  88         elif method == 'max':
  89             sbtrack.change_score_method( ord('M') )
  90         else:
  91             raise Exception("Can't reach here!")
  92
  93         info("Write bedGraph of scores...")
  94         ofhd = open(ofile,"w")
  95         sbtrack.write_bedGraph(ofhd,name="%s_Scores" % (method.upper()),description="Scores calculated by %s" % (method.upper()), column = 3)
  96         info("Finished '%s'! Please check '%s'!" % (method, ofile))