docs: def estimate_shift(signal, genome=None, windowsize=5000, thresh=None,
                     nwindows=1000, maxlag=500, array_kwargs=None,
                     verbose=False):
      """
      Experimental: cross-correlation to estimate the shift width of ChIP-seq
      data
  
      This can be interpreted as the binding site footprint.
  
      For ChIP-seq, the plus and minus strand reads tend to be shifted in the 5'
      direction away from each other.  Various ChIP-seq peak-callers estimate
      this distance; this function provides a quick, tunable way to do so using
      cross-correlation.  The resulting shift can then be incorporated into
      subsequent calls to `array` by adding the shift_width kwarg.
  
  
      :param signal: genomic_signal object
      :param genome: String assembly for constructing windows
      :param nwindows: Number of windows to compute cross-correlation on
      :param windowsize: Size of each window to compute cross-correlation on.
      :param thresh: Threshold read coverage to run cross-correlation on.  This
          is likely to be a function of the fragment size provided in
          `array_kwargs` `windowsize`.  If `thresh` is small, then the cross
          correlation can be noisy.
      :param maxlag: Max shift to look for
      :param array_kwargs: Kwargs passed directly to genomic_signal.array, with
          the default of `bins=windowsize` for single-bp resolution, and
          `read_strand` will be overwritten.
      :param verbose: Be verbose.
  
      Returns lags and a `maxlag*2+1` x `nwindows` matrix of cross-correlations.
  
      You can then plot the average cross-correlation function with::
  
          plt.plot(lags, shift.mean(axis=0))
  
      and get the distance to shift with::
  
          d = lags[np.argmax(shift.mean(axis=0))]
  
      and then plot that with::
  
          plt.axvline(d, color='k', linestyle='--')
  
      The number of windows with at least `thresh` coverage is::
  
          shift.shape[0]
      """
      if thresh is None:
          thresh = 0
  
      if genome is None:
          genome = signal.genome()
  
      if array_kwargs is None:
          array_kwargs = {}
  
      array_kwargs.pop('read_strand', None)
  
      if 'bins' not in array_kwargs:
          array_kwargs['bins'] = windowsize
  
      def add_strand(f, strand):
          fields = f.fields[:]
          while len(fields) < 5:
              fields.append('.')
          fields.append(strand)
          return pybedtools.create_interval_from_list(fields)
  
      windows = pybedtools.BedTool()\
          .window_maker(genome=genome, w=windowsize)
  
      random_subset = pybedtools.BedTool(windows[:nwindows])\
          .shuffle(genome=genome).saveas()
  
      if verbose:
          sys.stderr.write("Getting plus-strand signal for %s regions...\n"
                           % nwindows)
          sys.stderr.flush()
  
      plus = signal.array(
          features=random_subset,
          read_strand="+",
          **array_kwargs).astype(float)
  
      if verbose:
          sys.stderr.write("Getting minus-strand signal for %s regions...\n"
                           % nwindows)
          sys.stderr.flush()
  
      minus = signal.array(
          features=random_subset,
          read_strand="-",
          **array_kwargs).astype(float)
  
      # only do cross-correlation if you have enough reads to do so
      enough = ((plus.sum(axis=1) / windowsize) > thresh) \
          & ((minus.sum(axis=1) / windowsize) > thresh)
  
      if verbose:
          sys.stderr.write(
              "Running cross-correlation on %s regions that passed "
              "threshold\n" % sum(enough))
      results = np.zeros((sum(enough), 2 * maxlag + 1))
      for i, xy in enumerate(izip(plus[enough], minus[enough])):
          x, y = xy
          results[i] = xcorr(x, y, maxlag)
  
      lags = np.arange(-maxlag, maxlag + 1)
  
      return lags, results