1# Time-stamp: <2019-09-25 10:04:08 taoliu>
2
3"""Description: Naive call peaks from a single bedGraph track for
4scores.
5
6This code is free software; you can redistribute it and/or modify it
7under the terms of the BSD License (see the file LICENSE included with
8the distribution).
9"""
10
11# ------------------------------------
12# python modules
13# ------------------------------------
14import sys
15import os
16import logging
17from MACS2.IO import BedGraphIO
18# ------------------------------------
19# constants
20# ------------------------------------
21logging.basicConfig(level=20,
22                    format='%(levelname)-5s @ %(asctime)s: %(message)s ',
23                    datefmt='%a, %d %b %Y %H:%M:%S',
24                    stream=sys.stderr,
25                    filemode="w"
26                    )
27
28# ------------------------------------
29# Misc functions
30# ------------------------------------
31error   = logging.critical		# function alias
32warn    = logging.warning
33debug   = logging.debug
34info    = logging.info
35# ------------------------------------
36# Classes
37# ------------------------------------
38
39# ------------------------------------
40# Main function
41# ------------------------------------
42def run( options ):
43    info("Read and build bedGraph...")
44    bio = BedGraphIO.bedGraphIO(options.ifile)
45    btrack = bio.build_bdgtrack(baseline_value=0)
46
47    if options.cutoff_analysis:
48        info("Analyze cutoff vs number of peaks/total length of peaks/average length of peak")
49        cutoff_analysis_result = btrack.cutoff_analysis( int(options.maxgap), int(options.minlen), 50 )
50        info("Write report...")
51        if options.ofile:
52            fhd = open( os.path.join( options.outdir, options.ofile ), 'w' )
53        else:
54            fhd = open ( os.path.join( options.outdir, "%s_l%d_g%d_cutoff_analysis.txt" % (options.oprefix,options.minlen,options.maxgap)), "w" )
55        fhd.write( cutoff_analysis_result )
56        info("Done")
57    else:
58        info("Call peaks from bedGraph...")
59        peaks = btrack.call_peaks(cutoff=float(options.cutoff),min_length=int(options.minlen),max_gap=int(options.maxgap),call_summits=options.call_summits)
60
61        info("Write peaks...")
62        if options.ofile:
63            options.oprefix = options.ofile
64            nf = open( os.path.join( options.outdir, options.ofile ), 'w' )
65        else:
66            nf = open ( os.path.join( options.outdir, "%s_c%.1f_l%d_g%d_peaks.narrowPeak" % (options.oprefix,options.cutoff,options.minlen,options.maxgap)), "w" )
67        peaks.write_to_narrowPeak(nf, name=options.oprefix.encode(), name_prefix=(options.oprefix+"_narrowPeak").encode(), score_column="score", trackline=options.trackline)
68        info("Done")
69
70
71
72