MACS2/pileup_cmd.py

   1 # Time-stamp: <2019-09-25 10:03:31 taoliu>
   2
   3 """Description: Pileup alignment files
   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_pileup as opt_validate
  22 from MACS2.Pileup import pileup_and_write, pileup_and_write_pe
  23 from MACS2.Constants import *
  24 # ------------------------------------
  25 # Main function
  26 # ------------------------------------
  27 def run( o_options ):
  28     """The Main function/pipeline for duplication filter.
  29
  30     """
  31     # Parse options...
  32     options = opt_validate( o_options )
  33     # end of parsing commandline options
  34     info = options.info
  35     warn = options.warn
  36     debug = options.debug
  37     error = options.error
  38     #0 output arguments
  39     options.PE_MODE = options.format in ('BAMPE','BEDPE')
  40     #assert options.format != 'BAMPE', "Pair-end data with BAMPE option currently doesn't work with pileup command. You can pretend your data to be single-end with -f BAM. Please try again!"
  41
  42
  43     #0 prepare output file
  44     outfile = os.path.join( options.outdir, options.outputfile ).encode()
  45     if os.path.isfile( outfile ):
  46         info("# Existing file %s will be replaced!" % outfile )
  47         os.unlink( outfile )
  48
  49     #1 Read tag files
  50     info("# read alignment files...")
  51     if options.PE_MODE:
  52         info("# read input file in Paired-end mode.")
  53         treat = load_frag_files_options ( options ) # return PETrackI object
  54         t0 = treat.total # total fragments
  55         info("# total fragments/pairs in alignment file: %d" % (t0) )
  56         info("# Pileup paired-end alignment file.")
  57         pileup_and_write_pe(treat, outfile )
  58
  59     else:
  60         (tsize, treat) = load_tag_files_options  (options)
  61
  62         info("# tag size = %d", tsize)
  63
  64         t0 = treat.total
  65         info("# total tags in alignment file: %d", t0)
  66
  67         if options.bothdirection:
  68             info("# Pileup alignment file, extend each read towards up/downstream direction with %d bps" % options.extsize)
  69             pileup_and_write(treat, outfile, options.extsize * 2, 1, directional=False, halfextension=False)
  70         else:
  71             info("# Pileup alignment file, extend each read towards downstream direction with %d bps" % options.extsize)
  72             pileup_and_write(treat, outfile, options.extsize, 1, directional=True, halfextension=False)
  73
  74     info("# Done! Check %s" % options.outputfile)
  75
  76 def load_tag_files_options ( options ):
  77     """From the options, load alignment tags.
  78
  79     """
  80     options.info("# read treatment tags...")
  81     tp = options.parser(options.ifile[0], buffer_size=options.buffer_size)
  82     tsize = tp.tsize()
  83     treat = tp.build_fwtrack()
  84     #treat.sort()
  85     if len(options.ifile) > 1:
  86         # multiple input
  87         for tfile in options.ifile[1:]:
  88             tp = options.parser(tfile, buffer_size=options.buffer_size)
  89             treat = tp.append_fwtrack( treat )
  90             #treat.sort()
  91     treat.finalize()
  92
  93     options.info("tag size is determined as %d bps" % tsize)
  94     return (tsize, treat)
  95
  96 def load_frag_files_options ( options ):
  97     """From the options, load treatment fragments and control fragments (if available).
  98
  99     """
 100     options.info("# read treatment fragments...")
 101
 102     tp = options.parser(options.ifile[0], buffer_size=options.buffer_size)
 103     treat = tp.build_petrack()
 104     #treat.sort()
 105     if len(options.ifile) > 1:
 106         # multiple input
 107         for tfile in options.ifile[1:]:
 108             tp = options.parser(tfile, buffer_size=options.buffer_size)
 109             treat = tp.append_petrack( treat )
 110             #treat.sort()
 111     treat.finalize()
 112     return treat