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]

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