MACS3/Commands/bdgcmp_cmd.py

   1 # Time-stamp: <2024-10-02 16:06:33 Tao Liu>
   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 os
  11 from MACS3.IO import BedGraphIO
  12 from MACS3.Utilities.OptValidator import opt_validate_bdgcmp
  13
  14 # ------------------------------------
  15 # constants
  16 # ------------------------------------
  17
  18 # ------------------------------------
  19 # Misc functions
  20 # ------------------------------------
  21
  22 # ------------------------------------
  23 # Main function
  24 # ------------------------------------
  25
  26
  27 def run(options):
  28     options = opt_validate_bdgcmp(options)
  29     info = options.info
  30     # warn = options.warn
  31     # debug = options.debug
  32     # error = options.error
  33
  34     scaling_factor = options.sfactor
  35
  36     # not an actual depth, but its reciprocal, a trick to override
  37     # SPMR while necessary.
  38     pseudo_depth = 1.0/scaling_factor
  39
  40     info("Read and build treatment bedGraph...")
  41     tbio = BedGraphIO.bedGraphIO(options.tfile)
  42     tbtrack = tbio.read_bedGraph()
  43
  44     info("Read and build control bedGraph...")
  45     cbio = BedGraphIO.bedGraphIO(options.cfile)
  46     cbtrack = cbio.read_bedGraph()
  47
  48     info("Build ScoreTrackII...")
  49     sbtrack = tbtrack.make_ScoreTrackII_for_macs(cbtrack, depth1=pseudo_depth, depth2=pseudo_depth)
  50     if abs(scaling_factor-1) > 1e-6:
  51         # Only for the case while your input is SPMR from MACS3 callpeak; Let's override SPMR.
  52         info("Values in your input bedGraph files will be multiplied by %f ..." % scaling_factor)
  53         sbtrack.change_normalization_method(ord('M'))  # a hack to override SPMR
  54     sbtrack.set_pseudocount(options.pseudocount)
  55
  56     already_processed_method_list = []
  57     for (i, method) in enumerate(options.method):
  58         if method in already_processed_method_list:
  59             continue
  60         else:
  61             already_processed_method_list.append(method)
  62
  63         info("Calculate scores comparing treatment and control by '%s'..." % method)
  64         if options.ofile:
  65             ofile = os.path.join(options.outdir, options.ofile[i])
  66         else:
  67             ofile = os.path.join(options.outdir, options.oprefix + "_" + method + ".bdg")
  68         # build score track
  69         if method == 'ppois':
  70             sbtrack.change_score_method(ord('p'))
  71         elif method == 'qpois':
  72             sbtrack.change_score_method(ord('q'))
  73         elif method == 'subtract':
  74             sbtrack.change_score_method(ord('d'))
  75         elif method == 'logFE':
  76             sbtrack.change_score_method(ord('f'))
  77         elif method == 'FE':
  78             sbtrack.change_score_method(ord('F'))
  79         elif method == 'logLR':             # log likelihood
  80             sbtrack.change_score_method(ord('l'))
  81         elif method == 'slogLR':             # log likelihood
  82             sbtrack.change_score_method(ord('s'))
  83         elif method == 'max':
  84             sbtrack.change_score_method(ord('M'))
  85         else:
  86             raise Exception("Can't reach here!")
  87
  88         info("Write bedGraph of scores...")
  89         ofhd = open(ofile, "w")
  90         # write_bedGraph function for ScoreTrack
  91         sbtrack.write_bedGraph(ofhd, name="%s_Scores" % (method.upper()), description="Scores calculated by %s" % (method.upper()), column=3)
  92         info("Finished '%s'! Please check '%s'!" % (method, ofile))