MACS3/Commands/bdgdiff_cmd.py

   1 # Time-stamp: <2024-10-02 16:11:19 Tao Liu>
   2
   3 """Description: Naive call differential peaks from 4 bedGraph tracks for scores.
   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 # ------------------------------------
  11 # python modules
  12 # ------------------------------------
  13
  14 import os
  15 from MACS3.IO import BedGraphIO
  16 from MACS3.Signal import ScoreTrack
  17
  18 # ------------------------------------
  19 # constants
  20 # ------------------------------------
  21
  22 # ------------------------------------
  23 # Misc functions
  24 # ------------------------------------
  25 from MACS3.Utilities.Logger import logging
  26
  27 logger = logging.getLogger(__name__)
  28 debug = logger.debug
  29 info = logger.info
  30 error = logger.critical
  31 warn = logger.warning
  32 # ------------------------------------
  33 # Classes
  34 # ------------------------------------
  35
  36 # ------------------------------------
  37 # Main function
  38 # ------------------------------------
  39
  40
  41 def run(options):
  42     if options.maxgap >= options.minlen:
  43         error("MAXGAP should be smaller than MINLEN! Your input is MAXGAP = %d and MINLEN = %d" % (options.maxgap, options.minlen))
  44
  45     info("Read and build treatment 1 bedGraph...")
  46     t1bio = BedGraphIO.bedGraphIO(options.t1bdg)
  47     t1btrack = t1bio.read_bedGraph()
  48
  49     info("Read and build control 1 bedGraph...")
  50     c1bio = BedGraphIO.bedGraphIO(options.c1bdg)
  51     c1btrack = c1bio.read_bedGraph()
  52
  53     info("Read and build treatment 2 bedGraph...")
  54     t2bio = BedGraphIO.bedGraphIO(options.t2bdg)
  55     t2btrack = t2bio.read_bedGraph()
  56
  57     info("Read and build control 2 bedGraph...")
  58     c2bio = BedGraphIO.bedGraphIO(options.c2bdg)
  59     c2btrack = c2bio.read_bedGraph()
  60
  61     depth1 = options.depth1
  62     depth2 = options.depth2
  63
  64     if depth1 > depth2:         # scale down condition 1 to size of condition 2
  65         depth1 = depth2 / depth1
  66         depth2 = 1.0
  67     elif depth1 < depth2:       # scale down condition 2 to size of condition 1
  68         depth2 = depth1 / depth2
  69         depth1 = 1.0
  70     else:                       # no need to scale down any
  71         depth1 = 1.0
  72         depth2 = 1.0
  73
  74     twoconditionscore = ScoreTrack.TwoConditionScores(t1btrack,
  75                                                       c1btrack,
  76                                                       t2btrack,
  77                                                       c2btrack,
  78                                                       depth1,
  79                                                       depth2)
  80     twoconditionscore.build()
  81     twoconditionscore.finalize()
  82     (cat1, cat2, cat3) = twoconditionscore.call_peaks(min_length=options.minlen, max_gap=options.maxgap, cutoff=options.cutoff)
  83
  84     info("Write peaks...")
  85
  86     if options.ofile:
  87         ofiles = [os.path.join(options.outdir, x) for x in options.ofile]
  88         name_prefix = [x.encode() for x in options.ofile]
  89     else:
  90         ofiles = [os.path.join(options.outdir, "%s_c%.1f_cond1.bed" % (options.oprefix, options.cutoff)),
  91                   os.path.join(options.outdir, "%s_c%.1f_cond2.bed" % (options.oprefix, options.cutoff)),
  92                   os.path.join(options.outdir, "%s_c%.1f_common.bed" % (options.oprefix, options.cutoff))
  93                   ]
  94         name_prefix = [x.encode() for x in [options.oprefix+"_cond1_", options.oprefix+"_cond2_", options.oprefix+"_common_"]]
  95
  96     nf = open(ofiles[0], 'w')
  97     cat1.write_to_bed(nf, name_prefix=name_prefix[0], name=b"condition 1", description=b"unique regions in condition 1", score_column="score")
  98     nf.close()
  99
 100     nf = open(ofiles[1], 'w')
 101     cat2.write_to_bed(nf, name_prefix=name_prefix[1], name=b"condition 2", description=b"unique regions in condition 2", score_column="score")
 102     nf.close()
 103
 104     nf = open(ofiles[2], 'w')
 105     cat3.write_to_bed(nf, name_prefix=name_prefix[2], name=b"common", description=b"common regions in both conditions", score_column="score")
 106     nf.close()
 107     info("Done")