/usr/bin/bdgbroadcall 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 | #! /usr/bin/python
# Time-stamp: <2011-09-08 11:57:20 Tao Liu>
"""Description: Fine-tuning script to call broad peaks from a single bedGraph track 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 sys
import logging
from optparse import OptionParser
from MACS2.IO import bedGraphIO
# ------------------------------------
# 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 <-i bedGraph> [-c CUTOFF1] [-C CUTOFF2] [-l MIN] [-g MAX1] [-G MAX2] [-o PREFIX]"
description = "Call broad peaks from MACS pvalue or qscore score bedGraph output, with customized settings. Output two files for narrow peaks in encodePeak format, and one for broad peaks in bed12 format."
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("-i","--ifile",dest="ifile",type="string",
help="MACS pvalue score bedGraph")
optparser.add_option("-c","--cutoff-peak",dest="cutoffpeak",type="float",
help="Cutoff for peaks depending on which method you used for score track. If the file contains qvalue scores from MACS2, score 2 means qvalue 0.01. DEFAULT: 2",default=2)
optparser.add_option("-C","--cutoff-link",dest="cutofflink",type="float",
help="Cutoff for linking regions/low abundance regions depending on which method you used for score track. If the file contains qvalue scores from MACS2, score 1 means qvalue 0.1, and score 0.3 means qvalue 0.5. DEFAULT: 5",default=1)
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","--lvl1-max-gap",dest="lvl1maxgap",type="int",
help="maximum gap between significant peaks, better to set it as tag size. DEFAULT: 30",default=30)
optparser.add_option("-G","--lvl2-max-gap",dest="lvl2maxgap",type="int",
help="maximum linking between significant peaks, better to set it as 4 times of d value. DEFAULT: 800",default=800)
optparser.add_option("-o","--o-prefix",dest="oprefix",default="peak",
help="output file prefix, DEFAULT: peak")
(options,args) = optparser.parse_args()
if not options.ifile:
optparser.print_help()
sys.exit()
info("Read and build bedGraph...")
bio = bedGraphIO.bedGraphIO(options.ifile)
btrack = bio.build_bdgtrack(baseline_value=0)
info("Call peaks from bedGraph...")
(peaks,bpeaks) = btrack.call_broadpeaks (lvl1_cutoff=options.cutoffpeak, lvl2_cutoff=options.cutofflink, min_length=options.minlen, lvl1_max_gap=options.lvl1maxgap, lvl2_max_gap=options.lvl2maxgap)
info("Write peaks...")
nf = open ("%s_c%.0f_l%d_g%d_peaks.encodePeak" % (options.oprefix,options.cutoffpeak,options.minlen,options.lvl1maxgap),"w")
bf = open ("%s_c%.0f_C%.0f_l%d_g%d_G%d_broad.bed" % (options.oprefix,options.cutoffpeak,options.cutofflink,options.minlen,options.lvl1maxgap,options.lvl2maxgap),"w")
peaks.write_to_narrowPeak(nf, name_prefix=options.oprefix+"_encodePeak", score_column="score")
bpeaks.write_to_gappedPeak(bf, name_prefix=options.oprefix+"_broadRegion")
info("Done")
if __name__ == '__main__':
try:
main()
except KeyboardInterrupt:
sys.stderr.write("User interrupt me! ;-) See you!\n")
sys.exit(0)
|