MACS2/bdgdiff_cmd.py

   1 # Time-stamp: <2019-09-25 12:27:01 taoliu>
   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 sys
  15 import os
  16 import logging
  17 from MACS2.IO import BedGraphIO
  18 from MACS2.IO import ScoreTrack
  19
  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     if options.maxgap >= options.minlen:
  46         error("MAXGAP should be smaller than MINLEN! Your input is MAXGAP = %d and MINLEN = %d" % (options.maxgap, options.minlen))
  47
  48     LLR_cutoff = options.cutoff
  49     ofile_prefix = options.oprefix
  50
  51     info("Read and build treatment 1 bedGraph...")
  52     t1bio = BedGraphIO.bedGraphIO(options.t1bdg)
  53     t1btrack = t1bio.build_bdgtrack()
  54
  55     info("Read and build control 1 bedGraph...")
  56     c1bio = BedGraphIO.bedGraphIO(options.c1bdg)
  57     c1btrack = c1bio.build_bdgtrack()
  58
  59     info("Read and build treatment 2 bedGraph...")
  60     t2bio = BedGraphIO.bedGraphIO(options.t2bdg)
  61     t2btrack = t2bio.build_bdgtrack()
  62
  63     info("Read and build control 2 bedGraph...")
  64     c2bio = BedGraphIO.bedGraphIO(options.c2bdg)
  65     c2btrack = c2bio.build_bdgtrack()
  66
  67     depth1 = options.depth1
  68     depth2 = options.depth2
  69
  70     if depth1 > depth2:         # scale down condition 1 to size of condition 2
  71         depth1 = depth2 / depth1
  72         depth2 = 1.0
  73     elif depth1 < depth2:       # scale down condition 2 to size of condition 1
  74         depth2 = depth1/ depth2
  75         depth1 = 1.0
  76     else:                       # no need to scale down any
  77         depth1 = 1.0
  78         depth2 = 1.0
  79
  80     twoconditionscore = ScoreTrack.TwoConditionScores( t1btrack,
  81                                                         c1btrack,
  82                                                         t2btrack,
  83                                                         c2btrack,
  84                                                         depth1,
  85                                                         depth2 )
  86     twoconditionscore.build()
  87     twoconditionscore.finalize()
  88     (cat1,cat2,cat3) = twoconditionscore.call_peaks(min_length=options.minlen, max_gap=options.maxgap, cutoff=options.cutoff)
  89
  90     info("Write peaks...")
  91
  92     if options.ofile:
  93         ofiles = [os.path.join( options.outdir, x ) for x in options.ofile]
  94         name_prefix = [ x.encode() for x in options.ofile ]
  95     else:
  96         ofiles = [ os.path.join( options.outdir, "%s_c%.1f_cond1.bed" % (options.oprefix,options.cutoff)),
  97                    os.path.join( options.outdir, "%s_c%.1f_cond2.bed" % (options.oprefix,options.cutoff)),
  98                    os.path.join( options.outdir, "%s_c%.1f_common.bed" % (options.oprefix,options.cutoff))
  99                    ]
 100         name_prefix = [ x.encode() for x in [ options.oprefix+"_cond1_", options.oprefix+"_cond2_", options.oprefix+"_common_" ]]
 101
 102     nf = open( ofiles[ 0 ], 'w' )
 103     cat1.write_to_bed(nf, name_prefix=name_prefix[ 0 ], name=b"condition 1", description=b"unique regions in condition 1", score_column="score")
 104     nf.close()
 105
 106     nf = open( ofiles[ 1 ], 'w' )
 107     cat2.write_to_bed(nf, name_prefix=name_prefix[ 1 ], name=b"condition 2", description=b"unique regions in condition 2", score_column="score")
 108     nf.close()
 109
 110     nf = open( ofiles[ 2 ], 'w' )
 111     cat3.write_to_bed(nf, name_prefix=name_prefix[ 2 ], name=b"common", description=b"common regions in both conditions", score_column="score")
 112     nf.close()
 113     info("Done")
 114