/usr/bin/bdgdiff is in macs 2.0.9.1-1.
This file is owned by root:root, with mode 0o755.
The actual contents of the file can be viewed below.
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 | #! /usr/bin/python
# Time-stamp: <2011-08-23 12:22:55 Tao Liu>
"""Description: Naive call differential peaks from 4 bedGraph tracks for scores.
Copyright (c) 2011 Tao Liu <taoliu@jimmy.harvard.edu>
This code is free software; you can redistribute it and/or modify it
under the terms of the BSD License (see the file COPYING included with
the distribution).
@status: experimental
@version: $Revision$
@author: Tao Liu
@contact: taoliu@jimmy.harvard.edu
"""
# ------------------------------------
# python modules
# ------------------------------------
import os
import sys
import re
import logging
from optparse import OptionParser
from MACS2.IO import bedGraphIO
from MACS2.IO.cCompositeScoreTrack import *
# ------------------------------------
# constants
# ------------------------------------
logging.basicConfig(level=20,
format='%(levelname)-5s @ %(asctime)s: %(message)s ',
datefmt='%a, %d %b %Y %H:%M:%S',
stream=sys.stderr,
filemode="w"
)
# ------------------------------------
# Misc functions
# ------------------------------------
error = logging.critical # function alias
warn = logging.warning
debug = logging.debug
info = logging.info
# ------------------------------------
# Classes
# ------------------------------------
# ------------------------------------
# Main function
# ------------------------------------
def main():
usage = "usage: %prog <--s1 bedGraph> <--s2 bedGraph> <--s12 bedGraph> <--s21 bedGraph> [-c CUTOFF] [-l MIN] [-g MAX] [-o PREFIX]"
description = "Call differential regions from four bedGraph files of pair-wise p/qvalue, with customized settings. Please read the descriptions of options carefully."
optparser = OptionParser(version="%prog 0.1",description=description,usage=usage,add_help_option=False)
optparser.add_option("-h","--help",action="help",help="Show this help message and exit.")
optparser.add_option("--s1",dest="st1c1bdg",type="string",
help="MACS p/qvalue score bedGraph, condition 1 ChIP treatment against condition 1 control input.")
optparser.add_option("--s2",dest="st2c2bdg",type="string",
help="MACS p/qvalue score bedGraph, condition 2 ChIP treatment against condition 2 control input.")
optparser.add_option("--s12",dest="st1t2bdg",type="string",
help="MACS p/qvalue score bedGraph, condition 1 ChIP treatment against condition 2 ChIP treatment.")
optparser.add_option("--s21",dest="st2t1bdg",type="string",
help="MACS p/qvalue score bedGraph, condition 2 ChIP treatment against condition 1 ChIP treatment.")
optparser.add_option("-c","--cutoff",dest="cutoff",type="float",
help="Cutoff depends on which method you used for score track. If the file contains p/qvalue scores from MACS2, score 5 means p/qvalue 1e-5. DEFAULT: 5",default=5)
optparser.add_option("-l","--min-length",dest="minlen",type="int",
help="minimum length of peak, better to set it as d value. DEFAULT: 200",default=200)
optparser.add_option("-g","--max-gap",dest="maxgap",type="int",
help="maximum gap between significant points in a peak, better to set it as tag size. DEFAULT: 30",default=30)
optparser.add_option("-o","--o-prefix",dest="oprefix",default="peak",type="string",
help="output file prefix, DEFAULT: peak")
(options,args) = optparser.parse_args()
if not (options.st1c1bdg and options.st2c2bdg and options.st1t2bdg and options.st2t1bdg):
optparser.print_help()
sys.exit()
info("Read and build bedGraph...")
info("Score of condition 1 treatment vs condition 1 control...")
t1c1_bio = bedGraphIO.bedGraphIO(options.st1c1bdg)
t1c1_btrack = t1c1_bio.build_bdgtrack(baseline_value=0)
info("Score of condition 2 treatment vs condition 2 control...")
t2c2_bio = bedGraphIO.bedGraphIO(options.st2c2bdg)
t2c2_btrack = t2c2_bio.build_bdgtrack(baseline_value=0)
info("Score of condition 1 treatment vs condition 2 treatment...")
t1t2_bio = bedGraphIO.bedGraphIO(options.st1t2bdg)
t1t2_btrack = t1t2_bio.build_bdgtrack(baseline_value=0)
info("Score of condition 2 treatment vs condition 1 treatment...")
t2t1_bio = bedGraphIO.bedGraphIO(options.st2t1bdg)
t2t1_btrack = t2t1_bio.build_bdgtrack(baseline_value=0)
info("Combine four score tracks...")
comp_btrack = make_compositeScoreTrack(t1c1_btrack,t2c2_btrack,t1t2_btrack,t2t1_btrack)
info("Call differential regions ...")
( consistent_peaks, condition1_peaks, condition2_peaks ) = comp_btrack.call_diff_regions(cutoff=options.cutoff,min_length=options.minlen,max_gap=options.maxgap)
info("Write peaks...")
consistent_f = open ("%s_c%.0f_l%d_g%d_consistent_peaks.bed" % (options.oprefix,options.cutoff,options.minlen,options.maxgap),"w")
condition1_f = open ("%s_c%.0f_l%d_g%d_condition1_unique_peaks.bed" % (options.oprefix,options.cutoff,options.minlen,options.maxgap),"w")
condition2_f = open ("%s_c%.0f_l%d_g%d_condition2_unique_peaks.bed" % (options.oprefix,options.cutoff,options.minlen,options.maxgap),"w")
consistent_peaks.write_to_bed(consistent_f,name_prefix=options.oprefix+"_consistent_peak_", score_column="score")
condition1_peaks.write_to_bed(condition1_f,name_prefix=options.oprefix+"_condition1_unique_peak_", score_column="score")
condition2_peaks.write_to_bed(condition2_f,name_prefix=options.oprefix+"_condition2_unique_peak_", score_column="score")
info("Done")
if __name__ == '__main__':
try:
main()
except KeyboardInterrupt:
sys.stderr.write("User interrupt me! ;-) See you!\n")
sys.exit(0)
|