MACS2/predictd_cmd.py

   1 # Time-stamp: <2019-09-20 11:59:26 taoliu>
   2
   3 """Description: predict fragment size.
   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 import sys
  16 import logging
  17
  18 # ------------------------------------
  19 # own python modules
  20 # ------------------------------------
  21 from MACS2.OptValidator import opt_validate_predictd as opt_validate
  22 from MACS2.PeakModel import PeakModel,NotEnoughPairsException
  23 from MACS2.Prob import binomial_cdf_inv
  24 from MACS2.OutputWriter import model2r_script
  25 from MACS2.Constants import *
  26 # ------------------------------------
  27 # Main function
  28 # ------------------------------------
  29 def run( o_options ):
  30     """The Main function/pipeline for duplication filter.
  31
  32     """
  33     # Parse options...
  34     options = opt_validate( o_options )
  35     # end of parsing commandline options
  36     info = options.info
  37     warn = options.warn
  38     debug = options.debug
  39     error = options.error
  40     #0 output arguments
  41     assert options.format != 'BAMPE', "Pair-end data with BAMPE option doesn't work with predictd command. You can pretend your data to be single-end with -f BAM. Please try again!"
  42
  43     #1 Read tag files
  44     info("# read alignment files...")
  45     treat = load_tag_files_options  (options)
  46
  47     info("# tag size = %d", options.tsize)
  48
  49     t0 = treat.total
  50     info("# total tags in alignment file: %d", t0)
  51
  52     #2 Build Model
  53     info("# Build Peak Model...")
  54
  55     try:
  56         peakmodel = PeakModel(treatment = treat,
  57                               max_pairnum = MAX_PAIRNUM,
  58                               opt = options
  59                               )
  60         info("# finished!")
  61         debug("#  Summary Model:")
  62         debug("#   min_tags: %d" % (peakmodel.min_tags))
  63         debug("#   d: %d" % (peakmodel.d))
  64         info("# predicted fragment length is %d bps" % peakmodel.d)
  65         info("# alternative fragment length(s) may be %s bps" % ','.join(map(str,peakmodel.alternative_d)))
  66         info("# Generate R script for model : %s" % (options.modelR))
  67         model2r_script(peakmodel,options.modelR, options.rfile )
  68         options.d = peakmodel.d
  69
  70     except NotEnoughPairsException:
  71         warn("# Can't find enough pairs of symmetric peaks to build model!")
  72
  73 def load_tag_files_options ( options ):
  74     """From the options, load alignment tags.
  75
  76     """
  77     options.info("# read treatment tags...")
  78     tp = options.parser(options.ifile[0], buffer_size=options.buffer_size)
  79     if not options.tsize:           # override tsize if user specified --tsize
  80         ttsize = tp.tsize()
  81         options.tsize = ttsize
  82     treat = tp.build_fwtrack()
  83     #treat.sort()
  84     if len(options.ifile) > 1:
  85         # multiple input
  86         for tfile in options.ifile[1:]:
  87             tp = options.parser(tfile, buffer_size=options.buffer_size)
  88             treat = tp.append_fwtrack( treat )
  89             #treat.sort()
  90     treat.finalize()
  91
  92     options.info("tag size is determined as %d bps" % options.tsize)
  93     return treat