from __future__ import annotations
import tempfile
from textwrap import dedent
import shutil
import subprocess
import operator
import os
import sys
import random
import string
import pprint
from itertools import islice
from multiprocessing import Pool
import gzip
from typing import Any, Callable, Iterable, Iterator, Literal, Optional, TYPE_CHECKING, cast
import pysam
from warnings import warn
import pathlib
from pathlib import Path

from .helpers import (
    get_tempdir,
    _tags,
    call_bedtools,
    _flatten_list,
    _check_sequence_stderr,
    isBAM,
    isBGZIP,
    isGZIP,
    BEDToolsError,
    pybedtoolsError,
    _call_randomintersect,
    SplitOutput,
    FisherOutput,
)
from . import helpers
from .cbedtools import (
    IntervalFile,
    IntervalIterator,
    Interval,
    create_interval_from_list,
    BedToolsFileError,
)
import pybedtools
from . import settings
from . import filenames

if TYPE_CHECKING:
    import pandas as pd
    import matplotlib.colors as mcolors

_implicit_registry = {}
_other_registry = {}
_bam_registry = {}


def _jaccard_output_to_dict(s, **kwargs) -> dict:
    """
    jaccard method doesn't return an interval file, rather, it returns a short
    summary of results.  Here, we simply parse it into a dict for convenience.
    """
    if isinstance(s, str):
        _s = open(s).read()
    elif hasattr(s, "next") or hasattr(s, "__next__"):
        _s = "".join([i for i in s])
    else:
        raise ValueError("Unexpected object %r" % s)
    header, data = _s.splitlines()
    header = header.split()
    data = data.split()
    data[0] = int(data[0])
    data[1] = int(data[1])
    data[2] = float(data[2])
    data[3] = int(data[3])
    return dict(list(zip(header, data)))


def _reldist_output_handler(s, **kwargs):
    """
    reldist, if called with -detail, returns a valid BED file with the relative
    distance as the last field.  In that case, return the BedTool immediately.
    If not -detail, then the results are a table, in which case here we parse
    into a dict for convenience.
    """
    if "detail" in kwargs:
        return BedTool(s)
    if isinstance(s, str):
        iterable = open(s)
    if hasattr(s, "next"):
        iterable = s
    header = next(iterable).split()
    results = {}
    for h in header:
        results[h] = []
    for i in iterable:
        reldist, count, total, fraction = i.split()
        data = [float(reldist), int(count), int(total), float(fraction)]
        for h, d in zip(header, data):
            results[h].append(d)
    return results


def _wraps(
    prog: Optional[str] = None,
    implicit: Optional[str] = None,
    bam: Optional[str] = None,
    other: Optional[str] = None,
    uses_genome: bool = False,
    make_tempfile_for: Optional[str] = None,
    check_stderr:Optional[Callable]=None,
    add_to_bedtool:Optional[dict] = None,
    nonbam: Optional[Literal["ALL"] | str | list[str]] = None,
    force_bam: bool = False,
    genome_none_if: Optional[list[str]] =None,
    genome_if: Optional[list[str]] = None,
    genome_ok_if: Optional[list[str]] = None,
    does_not_return_bedtool: Optional[Callable] =None,
    arg_order: Optional[list[str]] = None,
):
    """
    Do-it-all wrapper, to be used as a decorator.

    *prog* is the name of the BEDTools program that will be called.  The help
    for this program will also be added to the decorated method's docstring.

    *implicit* is the BEDTools program arg that should be filled in
    automatically.

    *bam* will disable the implicit substitution if *bam* is in the kwargs.
    This is typically 'abam' or 'ibam' if the program accepts BAM input.

    *other* is the BEDTools program arg that is passed in as the second input,
    if supported.  Within the semantics of BEDTools, the typical case will be
    that if implicit='a' then other='b'; if implicit='i' then other=None.

    *uses_genome*, if True, will check for 'g' and/or 'genome' args and
    retrieve the corresponding genome files as needed.

    *make_tempfile_for* is used for the sequence methods and indicates which
    kwarg should have a tempfile made for it if it's not provided ('fo' for the
    sequence methods)

    *check_stderr*, if not None, is a function that accepts a string (which
    will be anything written to stdout when calling the wrapped program).  This
    function should return True if the string is OK, and False if it should
    truly be considered an error.  This is needed for wrapping fastaFromBed,
    which will report to stderr that it's creating an index file.

    *add_to_bedtool* is used for sequence methods.  It is a dictionary mapping
    kwargs to attributes to be created in the resulting BedTool.  Typically it
    is {'fo':'seqfn'} which will add the resulting sequence name to the
    BedTool's .seqfn attribute. If *add_to_bedtool* is not None, then the
    returned BedTool will be *self* with the added attribute.  If a key is
    "stdout" (e.g., {"stdout": attr_name}), then save the stdout of the command
    as a tempfile and store the tempfile's name in the attribute.  This is
    required for linksBed and bedToIgv.

    *nonbam* is a kwarg that even if the input file was a BAM, the output will
    *not* be BAM format.  For example, the `-bed` arg for intersectBed will
    cause the output to be in BED format, not BAM.  If not None, this can be a
    string, a list of strings, or the special string "ALL", which means that
    the wrapped program will never return BAM output.

    *force_bam*, if True, will force the output to be BAM.  This is used for
    bedToBam.

    *genome_none_if* is a list of arguments that will ignore the requirement
    for a genome.  This is needed for window_maker, where -b and -g are
    mutually exclusive.

    *genome_ok_if* is a list of arguments that, if they are in
    *genome_none_if*, are still OK to pass in.  This is needed for bedtool
    genomecov, where -g is not needed if -ibam is specified...but it's still OK
    if the user passes a genome arg.

    *genome_if* is a list of arguments that will trigger the requirement for
    a genome; otherwise no genome needs to be specified.

    *does_not_return_bedtool*, if not None, should be a function that handles
    the returned output.  Its signature should be ``func(output, kwargs)``,
    where `output` is the output from the [possibly streaming] call to BEDTools
    and `kwargs` are passed verbatim from the wrapped method call. Some
    examples of methods that use this are jaccard, reldist, fisher, and split
    methods.

    *arg_order*, if not None, is a sorted list of arguments. This is used by
    handle_kwargs() to deal with things like issues 81 and 345, where some
    BEDTools programs are sensitive to argument order.
    """

    # NOTE: We are calling each BEDTools program to get its help and adding
    # that to the docstring of each method. This is run at import time. However
    # if BEDTools is not on the path at import time, `not_implemented` is set
    # to True and isn't reset later until the module is reloaded.
    #
    # helpers.set_bedtools_path therefore will trigger a module reload.
    not_implemented = False

    # Call the program with -h to get help, which prints to stderr.
    try:
        p = subprocess.Popen(
            helpers._version_2_15_plus_names(prog) + ["-h"],
            stdout=subprocess.PIPE,
            stderr=subprocess.PIPE,
        )
        help_str = p.communicate()[1].decode()

        # underscores throw off ReStructuredText syntax of docstrings, so
        # replace 'em
        help_str = help_str.replace("_", "**")

        # indent
        help_str = help_str.split("\n")
        help_str = ["\n\n**Original BEDTools help:**::"] + ["\t" + i for i in help_str]
        help_str = "\n".join(help_str) + "\n"

    # If the program can't be found, then we'll eventually replace the method
    # with a version that does nothing but raise a NotImplementedError (plus
    # a helpful message).
    except OSError:
        help_str = (
            '"%s" does not appear to be installed '
            "or on the path, so this method is "
            "disabled.  Please install a more recent "
            "version of BEDTools and re-import to "
            "use this method." % prog
        )
        not_implemented = True

    def decorator(func):
        """
        Accepts a function to be wrapped; discards the original and returns a
        new, rebuilt-from-scratch function based on the kwargs passed to
        _wraps().
        """
        # Register the implicit (as well as bam and other) args in the global
        # registry.  BedTool.handle_kwargs() will access these at runtime.  The
        # registry is keyed by program name (like intersectBed).
        _implicit_registry[prog] = implicit
        if other is not None:
            _other_registry[prog] = other
        if bam is not None:
            _bam_registry[prog] = bam

        # Here's where we replace an unable-to-be-found program's method with
        # one that only returns a NotImplementedError
        if not_implemented:

            def not_implemented_func(*args, **kwargs):
                raise NotImplementedError(help_str)

            return not_implemented_func

        _add_doc = []
        if implicit:
            _add_doc.append(
                dedent(
                    """
                    For convenience, the file or stream this BedTool points to
                    is implicitly passed as the `-%s` argument to `%s`
                """
                    % (implicit, prog)
                )
            )

        if uses_genome:
            _add_doc.append(
                dedent(
                    """
                    There are two alternatives for supplying a genome.  Use
                    `g="genome.filename"` if you have a genome's chrom sizes
                    saved as a file. This is the what BEDTools expects when
                    using it from the command line. Alternatively, use the
                    `genome="assembly.name"` (for example, `genome="hg19"`) to
                    use chrom sizes for that assembly without having to manage
                    a separate file.  The `genome` argument triggers a call
                    `pybedtools.chromsizes`, so see that method for more
                    details.
                """
                )
            )

        def wrapped(self, *args, **kwargs):
            """
            A newly created function that will be returned by the _wraps()
            decorator
            """

            # Only one non-keyword argument is supported; this is then assumed
            # to be "other" (e.g., `-b` for intersectBed)
            if len(args) > 0:
                assert len(args) == 1
                kwargs[other] = args[0]

            # Add the implicit values to kwargs.  If the current BedTool is
            # BAM, it will automatically be passed to the appropriate
            # BAM-support arg (like `-abam`).  But this also allows the user to
            # explicitly specify the abam kwarg, which will override the
            # auto-substitution.
            # Note: here, `implicit` is something like "a"; `bam` is something
            # like "abam"
            if (
                (implicit not in kwargs)
                and (bam not in kwargs)
                and (implicit is not None)
            ):
                if not self._isbam:
                    kwargs[implicit] = self.fn
                else:
                    # It is a bam file.  If this program supports BAM as the
                    # first input, then we set it here
                    if bam is not None:
                        kwargs[bam] = self.fn

                    # Otherwise, BEDTools can't currently handle it, so raise
                    # an exception.
                    else:
                        raise pybedtoolsError(
                            '"%s" currently can\'t handle BAM '
                            "input, please use bam_to_bed() first." % prog
                        )

            # Should this function handle genome files?
            check_for_genome = uses_genome
            if uses_genome:
                if genome_none_if:
                    for i in genome_none_if:
                        if i in kwargs or i == implicit:
                            check_for_genome = False

                    # for genomecov, if -ibam then -g is optional.  So it's OK
                    # for the user to provide genome or g kwargs, even if
                    # -ibam.
                    if genome_ok_if:
                        for i in genome_ok_if:
                            if i in kwargs or i == implicit:
                                if ("g" in kwargs) or ("genome" in kwargs):
                                    check_for_genome = True
                if genome_if:
                    check_for_genome = False
                    for i in genome_if:
                        if (i in kwargs) or (i == implicit):
                            check_for_genome = True
            if check_for_genome:
                kwargs = self.check_genome(**kwargs)

            # For sequence methods, we may need to make a tempfile that will
            # hold the resulting sequence.  For example, fastaFromBed needs to
            # make a tempfile for 'fo' if no 'fo' was explicitly specified by
            # the user.
            if make_tempfile_for is not None:
                if make_tempfile_for not in kwargs:
                    kwargs[make_tempfile_for] = self._tmp()

            # At runtime, this will parse the kwargs, convert streams to
            # tempfiles if needed, and return all the goodies
            cmds, tmp, stdin = self.handle_kwargs(prog=prog,
                                                  arg_order=arg_order,
                                                  **kwargs)

            # Decide whether the output is BAM format or not.
            result_is_bam = False

            # By default, if the current BedTool is BAM, then the result should
            # be, too.
            if self._isbam:
                result_is_bam = True

            # If nonbam is "ALL", then this method will never return BAM
            # output.
            if nonbam == "ALL":
                result_is_bam = False

            # If any of the `nonbam` args are found in kwargs, then result is
            # not a BAM.  Side note: the _nonbam name mangling is necessary to
            # keep the nonbam arg passed into the original _wraps() decorator
            # in scope.
            if nonbam is not None and nonbam != "ALL":
                if isinstance(nonbam, str):
                    _nonbam = [nonbam]
                else:
                    _nonbam = nonbam
                for i in _nonbam:
                    if i in kwargs:
                        result_is_bam = False
                        break

            if force_bam:
                result_is_bam = True

            decode_output = not result_is_bam

            # Do the actual call
            stream = call_bedtools(
                cmds,
                tmp,
                stdin=stdin,
                check_stderr=check_stderr,
                decode_output=decode_output,
            )

            if does_not_return_bedtool:
                return does_not_return_bedtool(stream, **kwargs)

            # Post-hoc editing of the BedTool -- for example, this is used for
            # the sequence methods to add a `seqfn` attribute to the resulting
            # BedTool.
            if add_to_bedtool is not None:
                for kw, attr in list(add_to_bedtool.items()):
                    if kw == "stdout":
                        value = stream
                    else:
                        value = kwargs[kw]
                    setattr(self, attr, value)
                    result = self
            else:
                result = BedTool(stream)

            result._isbam = result_is_bam
            result._cmds = cmds
            del kwargs
            return result

        # Now add the edited docstring (original Python docstring plus BEDTools
        # help) to the newly created method above
        if func.__doc__ is None:
            orig = ""
        else:
            orig = func.__doc__

        wrapped.__doc__ = orig + "\n".join(_add_doc) + help_str

        # Add the original method's name to a new attribute so we can access it
        # when logging history
        wrapped._name = func.__name__  # type: ignore

        return wrapped

    return decorator


def _log_to_history(method: Callable):
    """
    Decorator to add a method and its kwargs to the history.

    Assumes that you only add this decorator to bedtool instances that
    return other bedtool instances
    """

    def decorated(self, *args, **kwargs):

        # this calls the actual method in the first place; *result* is
        # whatever you get back
        result = method(self, *args, **kwargs)

        # add appropriate tags
        parent_tag = self._tag
        result_tag = result._tag

        history_step = HistoryStep(
            method, args, kwargs, self, parent_tag, result_tag
        )

        # only add the current history to the new bedtool if there's
        # something to add
        if len(self.history) > 0:
            result.history.append(self.history)

        # but either way, add this history step to the result.
        result.history.append(history_step)

        return result

    decorated.__doc__ = method.__doc__
    return decorated


class BedTool(object):
    TEMPFILES = filenames.TEMPFILES

    def __init__(self, fn: Optional[Any] = None,
                 from_string: bool = False,
                 remote: bool = False):
        """
        Wrapper around Aaron Quinlan's ``BEDtools`` suite of programs
        (https://github.com/arq5x/bedtools); also contains many useful
        methods for more detailed work with BED files.

        *fn* is typically the name of a BED-like file, but can also be
        one of the following:

            * a string filename
            * another BedTool object
            * an iterable of Interval objects
            * an open file object
            * a "file contents" string (see below)

        If *from_string* is True, then you can pass a string that contains
        the contents of the BedTool you want to create.  This will treat all
        spaces as TABs and write to tempfile, treating whatever you pass as
        *fn* as the contents of the bed file.  This also strips empty lines.

        Typical usage is to point to an existing file::

            a = BedTool('a.bed')

        But you can also create one from scratch from a string::

            >>> s = '''
            ... chrX  1  100
            ... chrX 25  800
            ... '''
            >>> a = BedTool(s, from_string=True)

        Or use examples that come with pybedtools::

            >>> example_files = pybedtools.list_example_files()
            >>> assert 'a.bed' in example_files
            >>> a = pybedtools.example_bedtool('a.bed')

        """
        if remote:
            raise ValueError(
                "Remote BAM no longer supported (since BEDTools does not " "support it)"
            )
        self.remote = remote
        self._isbam = False
        self._bam_header = ""
        self._cmds = []
        if from_string:
            if fn is None or not isinstance(fn, str):
                raise ValueError("from_string=True requires a string to parse")
            bed_contents = fn
            fn = self._tmp()
            fout = open(fn, "w")
            for line in bed_contents.splitlines():
                if len(line.strip()) == 0:
                    continue
                line = "\t".join(line.split()) + "\n"
                fout.write(line)
            fout.close()

        else:
            # if fn is a Path object, we have to use its string representation
            if isinstance(fn, pathlib.PurePath):
                fn = str(fn)

            # our work is already done
            if isinstance(fn, BedTool):
                fn = fn.fn

            # from_string=False, so assume it's a filename
            elif isinstance(fn, str):
                if remote:
                    self._isbam = True
                else:
                    if not os.path.exists(fn):
                        msg = 'File "%s" does not exist' % fn
                        raise FileNotFoundError(msg)
                    self._isbam = isBAM(fn)

                # TODO: we dont' really need this, but it's added here for
                # compatibility with existing tests
                if self._isbam:
                    header = pysam.Samfile(fn).header.to_dict()
                    # For example:
                    # {
                    #     'HD': {'VN': '1.0', 'SO': 'coordinate'},
                    #     'SQ': [
                    #         {'LN': 23011544,
                    #          'SN': 'chr2L'},
                    #         {'LN': 21146708,
                    #          'SN': 'chr2R'},
                    #         {'LN': 24543557,
                    #          'SN': 'chr3L'},
                    #         {'LN': 27905053,
                    #          'SN': 'chr3R'},
                    #         {'LN': 1351857,
                    #          'SN': 'chr4'},
                    #         {'LN': 22422827,
                    #          'SN': 'chrX'}
                    #     ]
                    # }

                    txt_header = []
                    for k, v in header.items():
                        if isinstance(v, list):
                            for i in v:
                                if isinstance(i, dict):
                                    txt_header.append(
                                        "\t".join(
                                            ["@" + k]
                                            + [
                                                ":".join(map(str, j))
                                                for j in sorted(i.items(), reverse=True)
                                            ]
                                        )
                                    )
                                elif isinstance(i, str):
                                    txt_header.append(i)

                        elif isinstance(v, dict):
                            txt_header.append(
                                "\t".join(
                                    ["@" + k]
                                    + [
                                        ":".join(map(str, j))
                                        for j in sorted(v.items(), reverse=True)
                                    ]
                                )
                            )
                        else:
                            raise ValueError("unhandled type in BAM header")
                    self._bam_header = "\n".join(txt_header) + "\n"

            # If tuple or list, then save as file first
            # (fixes #73)
            elif isinstance(fn, (list, tuple)):
                fn = BedTool(iter(fn)).saveas().fn

            # Otherwise assume iterator, say an open file as from
            # subprocess.PIPE
            else:
                fn = fn

        self.fn = fn
        tag = "".join([random.choice(string.ascii_lowercase) for _ in range(8)])
        self._tag = tag
        _tags[tag] = self
        self._hascounts = False
        self._file_type = None
        self.seqfn = None
        self.fastq = None
        self.igv_script = None
        self.links_html = None
        self.history = History()

    @classmethod
    def from_dataframe(
        cls,
        df: pd.DataFrame,
        outfile: Optional[str] =None,
        sep: str ="\t",
        header: bool =False,
        na_rep:str =".",
        index:bool =False,
        **kwargs
    ) -> BedTool:
        """
        Creates a BedTool from a pandas.DataFrame.

        If `outfile` is None, a temporary file will be used. Otherwise it can
        be a specific filename or an open file handle. Additional kwargs will
        be passed to `pandas.DataFrame.to_csv`.

        The fields of the resulting BedTool will match the order of columns in
        the dataframe.
        """
        try:
            import pandas
        except ImportError:
            raise ImportError("pandas must be installed to use dataframes")
        if outfile is None:
            outfile = cls._tmp()
        default_kwargs = dict(sep=sep, header=header, na_rep=na_rep, index=index)
        default_kwargs.update(kwargs)
        df.to_csv(outfile, **default_kwargs)

        if isinstance(outfile, str):
            fn = outfile
        else:
            try:
                fn = outfile.name
            except AttributeError:
                raise ValueError(
                    "`outfile` is not a string and doesn't have a `name` attribute. "
                    "Unable to determine filename."
                )
        return BedTool(fn)

    def split(self, func: Callable, *args, **kwargs) -> BedTool:
        """
        Split each feature using a user-defined function.

        Calls the provided function `func` with each interval.  In contrast to
        `each` (which does something similar), this method expects `func` to
        return an *iterable* of Interval objects.

        args and kwargs are passed directly to `func`.

        Returns a new BedTool.
        """

        def generator():
            for orig_interval in self:
                for interval in func(orig_interval, *args, **kwargs):
                    yield interval

        return BedTool(generator())

    def truncate_to_chrom(self, genome: str | dict) -> BedTool:
        """
        Ensure all features fall within chromosome limits.

        Some peak-callers extend peaks such that the boundaries overstep
        chromosome coordinates.  Upon uploading such a file to a genome browser
        like UCSC, this results in an error like::

            Error line 101 of custom track: chromEnd larger than chrom chr2
            size

        Use this method to clean your file, truncating any out-of-bounds
        features to fit within the chromosome coordinates of `genome`.

        `genome` can be either an assembly name ('dm3') or a dictionary where
        keys are chrom and values are (start, stop) tuples.
        """
        if isinstance(genome, dict):
            chromdict = genome
        else:
            assert isinstance(genome, str)
            chromdict = helpers.chromsizes(genome)

        tmp = self._tmp()
        with open(tmp, "w") as fout:
            for chrom, coords in list(chromdict.items()):
                start, stop = coords
                start = str(start)
                stop = str(stop)
                fout.write("\t".join([chrom, start, stop]) + "\n")
        return self.intersect(tmp)

    def tabix_intervals(self, interval_or_string: Interval | str, check_coordinates: bool=False) -> BedTool:
        """
        Retrieve all intervals within coordinates from a "tabixed" BedTool.

        Given either a string in "chrom:start-stop" format, or an interval-like
        object with chrom, start, stop attributes, return a *streaming* BedTool
        of the features in this BedTool that overlap the provided interval.

        If the coordinates are invalid, an empty generator is returned unless
        `check_coordinates=True` in which case a ValueError will be raised.
        """
        if not self._tabixed():
            raise ValueError(
                "This BedTool has not been indexed for tabix "
                "-- please use the .tabix() method"
            )

        # tabix expects 1-based coords, but BEDTools works with
        # zero-based. pybedtools and pysam also work with zero-based. So we can
        # pass zero-based directly to the pysam tabix interface.
        tbx = pysam.TabixFile(self.fn)

        # If an interval is passed, use its coordinates directly
        if isinstance(interval_or_string, Interval):
            interval: Interval = interval_or_string
            chrom, start, end = interval.chrom, interval.start, interval.stop
        # Parse string directly instead of relying on Interval, in order to
        # permit full chromosome fetching
        else:
            match = helpers.coord_re.search(interval_or_string)
            # Assume string is contig if it doesn't fit chrom:start-end format
            if match is None:
                chrom = interval_or_string
                start, end = None, None
            # Otherwise parse the coordinates
            else:
                chrom, start, end = match.group(1, 2, 3)
                start, end = int(start), int(end)

        # Fetch results.
        try:
            results = tbx.fetch(str(chrom), start, end)
        except ValueError:
            if check_coordinates:
                raise
            else:
                results = []

        # pysam.ctabix.TabixIterator does not include newlines when yielding so
        # we need to add them.
        def gen():
            for i in results:
                yield i + "\n"

        # xref #190
        x = BedTool(gen()).saveas()
        tbx.close()
        return x

    def tabix_contigs(self):
        """
        Returns a list of contigs from the tabix index.
        """
        if not self._tabixed():
            raise ValueError(
                "This BedTool has not been indexed for tabix "
                "-- please use the .tabix() method"
            )

        tbx = pysam.TabixFile(self.fn)
        return tbx.contigs

    def tabix(self, in_place: bool = True, force: bool = False, is_sorted: bool = False) -> BedTool:
        """
        Prepare a BedTool for use with Tabix.

        Returns a new BedTool that has been BGZIP compressed
        and indexed by tabix.

        Parameters
        ----------

        in_place : bool
            If True (default), then assume the file is already sorted and
            replace the existing file with the BGZIPed version.

        force : bool
            If True (default is False), then overwrite both the index and the
            BGZIP file.

        is_sorted : bool
            If True (default is False), then assume the file is already sorted
            so that BedTool.bgzip() doesn't have to do that work.
        """
        # Return quickly if nothing to do
        if self._tabixed() and not force:
            return self

        # Make sure it's BGZIPed
        fn = self.bgzip(in_place=in_place, force=force, is_sorted=is_sorted)
        if self.file_type is not None and self.file_type not in ["bam", "empty"]:
            pysam.tabix_index(fn, force=force, preset=self.file_type)  # type: ignore
        return BedTool(fn)

    def _tabixed(self):
        """
        Verifies that we're working with a tabixed file: a string filename
        pointing to a BGZIPed file with a .tbi file in the same dir.
        """
        if (
            isinstance(self.fn, str)
            and isBGZIP(self.fn)
            and os.path.exists(self.fn + ".tbi")
        ):
            return True

    def bgzip(self, in_place: bool = True, force: bool = False, is_sorted: bool = False) -> str:
        """
        Helper function for more control over "tabixed" BedTools.

        Checks to see if we already have a BGZIP file; if not then prepare one.
        Always leaves the original file alone.  You can always just make a
        BedTool out of an already sorted and BGZIPed file to avoid this step.

        `in_place` will put the BGZIPed file in the same dir (possibly after
        sorting to tempfile).

        If `is_sorted`, then assume the file is already sorted. Otherwise call
        bedtools sort with the `-header` option.

        `force` will overwrite without asking.
        """
        # It may already be BGZIPed...
        if isinstance(self.fn, str) and not force:
            if isBGZIP(self.fn):
                return self.fn

        # If not in_place, then make a tempfile for the BGZIPed version
        if not in_place:
            # Get tempfile name, sorted or not
            if not is_sorted:
                fn = self.sort(header=True).fn
            else:
                fn = self._tmp()

            # Register for later deletion
            outfn = fn + ".gz"
            BedTool.TEMPFILES.append(outfn)

            # Creates tempfile.gz
            pysam.tabix_compress(fn, outfn, force=force)
            return outfn

        # Otherwise, make sure the BGZIPed version has a similar name to the
        # current BedTool's file
        if in_place:
            if not is_sorted:
                fn = self.sort(header=True).saveas().fn
            else:
                fn = self.fn
            outfn = self.fn + ".gz"
            pysam.tabix_compress(fn, outfn, force=force)
            return outfn

    def delete_temporary_history(self, ask: bool = True, raw_input_func=None):
        """
        Use at your own risk!  This method will delete temp files. You will be
        prompted for deletion of files unless you specify *ask=False*.

        Deletes all temporary files created during the history of this BedTool
        up to but not including the file this current BedTool points to.

        Any filenames that are in the history and have the following pattern
        will be deleted::

            <TEMP_DIR>/pybedtools.*.tmp

        (where <TEMP_DIR> is the result from get_tempdir() and is by default
        "/tmp")

        Any files that don't have this format will be left alone.

        (*raw_input_func* is used for testing)
        """
        flattened_history = _flatten_list(self.history)
        to_delete = []
        tempdir = get_tempdir()
        for i in flattened_history:
            fn = i.fn
            if fn.startswith(os.path.join(os.path.abspath(tempdir), "pybedtools")):
                if fn.endswith(".tmp"):
                    to_delete.append(fn)

        if raw_input_func is None:
            raw_input_func = input

        str_fns = "\n\t".join(to_delete)
        if ask:
            answer = raw_input_func("Delete these files?\n\t%s\n(y/N) " % str_fns)

            if not answer.lower()[0] == "y":
                print("OK, not deleting.")
                return
        for fn in to_delete:
            os.unlink(fn)
        return

    def filter(self, func: Callable, *args, **kwargs) -> BedTool:
        """
        Filter features by user-defined function.

        Takes a function *func* that is called for each feature in the
        `BedTool` object and returns only those for which the function returns
        True.

        *args and **kwargs are passed directly to *func*.

        Returns a streaming BedTool; if you want the filename then use the
        .saveas() method.

        >>> a = pybedtools.example_bedtool('a.bed')
        >>> subset = a.filter(lambda b: b.chrom == 'chr1' and b.start < 150)
        >>> len(a), len(subset)
        (4, 2)

        so it has extracted 2 records from the original 4.

        """
        return BedTool((f for f in self if func(f, *args, **kwargs)))

    def field_count(self, n:int=10) -> int:
        """
        Number of fields in each line of this BedTool (checks `n` lines)

        Return the number of fields in the features this file contains.  Checks
        the first *n* features.

        >>> a = pybedtools.example_bedtool('a.bed')
        >>> a.field_count()
        6
        """
        if self.file_type == "empty":
            return 0
        i = 0
        fields = set([])
        for feat in self:
            if i > n:
                break
            i += 1
            # TODO: make this more efficient.
            fields.update([len(feat.fields)])
        assert len(fields) == 1, fields
        return list(fields)[0]

    def each(self, func: Callable, *args, **kwargs) -> BedTool:
        """
        Modify each feature with a user-defined function.

        Applies user-defined function *func* to each feature.  *func* must
        accept an Interval as its first argument; *args and **kwargs will be
        passed to *func*.

        *func* must return an Interval object OR a value that evaluates to
        False, in which case the original feature will be removed from the
        output.  This way, an additional "filter" call is not necessary.

        >>> def truncate_feature(feature, limit=0):
        ...     feature.score = str(len(feature))
        ...     if len(feature) > limit:
        ...         feature.stop = feature.start + limit
        ...         feature.name = feature.name + '.short'
        ...     return feature

        >>> a = pybedtools.example_bedtool('a.bed')
        >>> b = a.each(truncate_feature, limit=100)
        >>> print(b) #doctest: +NORMALIZE_WHITESPACE
        chr1	1	100	feature1	99	+
        chr1	100	200	feature2	100	+
        chr1	150	250	feature3.short	350	-
        chr1	900	950	feature4	50	+
        <BLANKLINE>

        """

        def _generator():
            for f in self:
                result = func(f, *args, **kwargs)
                if result:
                    yield result

        return BedTool(_generator())

    def introns(self, gene: str = "gene", exon: str = "exon") -> BedTool:
        """
        Create intron features (requires specific input format).

        NOTE: this method assumes a simple file with non-overlapping exons. For
        more sophisticated features, consider the gffutils package instead.

        Given a BED12 or a GFF with exons, create a new `BedTool` with just
        introns.  The output is a bed6 file with the score column (5) being one
        of 'intron'/'utr5'/'utr3'
        """
        # iterate over all the features in the gene.
        s = self.sort()
        if self.file_type == "gff":
            exon_iter = BedTool((f for f in s if f[2] == exon)).saveas()
            gene_iter = BedTool((f for f in s if f[2] == gene)).saveas()

        elif self.file_type == "bed":
            if s.field_count() == 12:
                exon_iter = s.bed6().saveas()
                gene_iter = s.saveas()
            else:
                # TODO: bed6. groupby on name and find smallest start,
                # largest stop.
                exon_iter = s
                gene_iter = None
                raise NotImplementedError(
                    ".introns() only supported for bed12" "and GFF"
                )

        else:
            raise NotImplementedError(".introns() only supported for BED and GFF")

        with open(BedTool._tmp(), "w") as fh:
            # group on the name.
            exon_intervals = IntervalFile(exon_iter.fn)
            for g in gene_iter:
                # search finds all, but we just want the ones that completely
                # overlap this gene.
                exons = [
                    e
                    for e in exon_intervals.search(g, same_strand=True)
                    if e.start >= g.start and e.end <= g.end
                ]

                for i, exon_instance in enumerate(exons):
                    exon_instance: pybedtools.Interval
                    # 5' utr between gene start and first intron
                    if i == 0 and exon_instance.start > g.start:
                        utr = {"+": "utr5", "-": "utr3"}[g.strand]
                        print(
                            "%s\t%i\t%i\t%s\t%s\t%s"
                            % (g.chrom, g.start, exon_instance.start, g.name, utr, g.strand),
                            file=fh,
                        )
                    elif i == len(exons) - 1 and exon_instance.end < g.end:
                        utr = {"+": "utr3", "-": "utr5"}[g.strand]
                        print(
                            "%s\t%i\t%i\t%s\t%s\t%s"
                            % (g.chrom, exon_instance.end, g.end, g.name, utr, g.strand),
                            file=fh,
                        )
                    elif i != len(exons) - 1:
                        istart = exon_instance.end
                        iend = exons[i + 1].start
                        print(
                            "%s\t%i\t%i\t%s\tintron\t%s"
                            % (g.chrom, istart, iend, g.name, g.strand),
                            file=fh,
                        )
        return BedTool(fh.name)

    def features(self):
        """
        Returns an iterable of features
        """
        if hasattr(self, "next") or hasattr(self, "__next__"):
            return self
        return iter(self)

    FileType = Literal['bed', 'vcf', 'gff', 'bam', 'sam', 'empty']

    @property
    def file_type(self) -> Optional[FileType]:
        """
        Return the type of the current file.  One of ('bed','vcf','gff', 'bam',
        'sam', 'empty').

        >>> a = pybedtools.example_bedtool('a.bed')
        >>> print(a.file_type)
        bed
        """
        if not isinstance(self.fn, str):
            raise ValueError(
                "Checking file_type not supported for "
                "non-file BedTools. Use .saveas() to "
                "save as a temp file first."
            )
        if self._isbam:
            self._file_type = "bam"
        else:
            try:
                self._file_type = next(iter(self)).file_type
            except StopIteration:
                self._file_type = "empty"

        return self._file_type

    def cut(self, indexes: list[int], stream: bool = False) -> BedTool:
        """
        Analogous to unix `cut`.

        Similar to unix `cut` except indexes are 0-based, must be a list
        and the columns are returned in the order requested.

        This method returns a BedTool of results, which means that the indexes
        returned must be valid GFF/GTF/BED/SAM features.

        If you would like arbitrary columns -- say, just chrom and featuretype
        of a GFF, which would not comprise a valid feature -- then instead of
        this method, simply use indexes on each feature, e.g,

        >>> gff = pybedtools.example_bedtool('d.gff')
        >>> results = [(f[0], f[2]) for f in gff]

        In addition, `indexes` can contain keys of the GFF/GTF attributes, in
        which case the values are returned. e.g. 'gene_name' will return the
        corresponding name from a GTF, or 'start' will return the start
        attribute of a BED Interval.
        """
        if stream:
            return BedTool(([f[attr] for attr in indexes] for f in self))
        else:
            with open(self._tmp(), "w") as fh:
                for f in self:
                    print("\t".join(map(str, [f[attr] for attr in indexes])), file=fh)
            return BedTool(fh.name)

    @classmethod
    def _tmp(cls) -> str:
        """
        Makes a tempfile and registers it in the BedTool.TEMPFILES class
        variable.  Adds a "pybedtools." prefix and ".tmp" extension for easy
        deletion if you forget to call pybedtools.cleanup().
        """
        tmpfn = tempfile.NamedTemporaryFile(
            prefix=settings.tempfile_prefix,
            suffix=settings.tempfile_suffix,
            delete=False,
        )
        tmpfn = tmpfn.name
        cls.TEMPFILES.append(tmpfn)
        return tmpfn

    def __iter__(self):
        """
        Dispatches the right iterator depending on how this BedTool was
        created
        """
        if self._isbam:
            # Note: BAM class takes filename or stream, so self.fn is OK
            # here
            return BAM(self.fn)

        # Plain ol' filename
        if isinstance(self.fn, str):
            if not os.path.exists(self.fn):
                raise BedToolsFileError("{0} does not exist".format(self.fn))
            if isGZIP(self.fn):
                return IntervalIterator(gzip.open(self.fn, "rt"))
            else:
                return IntervalIterator(open(self.fn, "r"))
        # Any other kind of input (streaming string from stdout; iterable of
        # Intervals, iterable of (chrom, start, stop) tuples, etc are handled
        # appropriately by IntervalIterator.
        else:
            return IntervalIterator(self.fn)

    @property
    def intervals(self):
        if isinstance(self.fn, str):
            return IntervalFile(self.fn)
        else:
            raise ValueError("Please convert to a file-based BedTool using saveas")

    def __repr__(self):
        if isinstance(self.fn, str):
            if os.path.exists(self.fn) or self.remote:
                return "<BedTool(%s)>" % self.fn
            else:
                return "<BedTool(MISSING FILE: %s)>" % self.fn
        elif isinstance(self.fn, BedTool):
            return repr(self.fn)
        else:
            return "<BedTool(%s)>" % repr(self.fn)

    def __str__(self):
        """
        Returns the string representation of the whole `BedTool`
        """
        items = []
        for i in iter(self):
            i = str(i)
            if isinstance(i, bytes):
                i = i.decode("UTF-8")
            items.append(i)
        return "".join(items)

    def __len__(self):
        return self.count()

    def __eq__(self, other: object) -> bool:
        if isinstance(other, BedTool):
            if not isinstance(self.fn, str) or not isinstance(
                other.fn, str
            ):
                raise NotImplementedError(
                    "Testing equality only supported for"
                    " BedTools that point to files"
                )
        elif not isinstance(other, str):
            raise NotImplementedError(
                "Testing equality only supported for"
                " BedTools that point to files or str of content"
            )
        return str(self) == str(other)

    def __ne__(self, other:object):
        return not self.__eq__(other)

    def __getitem__(self, key: slice|int):
        if isinstance(key, slice):
            return islice(self, key.start, key.stop, key.step)
        elif isinstance(key, int):
            return list(islice(self, key, key + 1))[0]
        else:
            raise ValueError(
                "Only slices or integers allowed for indexing " "into a BedTool"
            )

    def __add__(self, other: BedTool) -> BedTool:
        try:
            result = self.intersect(other, u=True)
        except BEDToolsError as e:
            # BEDTools versions <2.20 would raise BEDToolsError
            if (self.file_type == "empty") or (other.file_type == "empty"):
                result = pybedtools.BedTool("", from_string=True)
            else:
                raise e
        return result

    def __sub__(self, other: BedTool) -> BedTool:
        result = None

        try:
            result = self.intersect(other, v=True)
        except BEDToolsError:
            # BEDTools versions <2.20 would raise BEDToolsError

            if (self.file_type == "empty") and (other.file_type == "empty"):
                result = pybedtools.BedTool("", from_string=True)
            elif other.file_type == "empty":
                result = self.saveas()
            elif self.file_type == "empty":
                result = pybedtools.BedTool("", from_string=True)
        if result is None:
            raise ValueError("Subtraction operation failed.")

        return result

    def head(self, n: int = 10, as_string: bool = False):
        """
        Prints the first *n* lines or returns them if as_string is True

        Note that this only opens the underlying file (gzipped or not), so it
        does not check to see if the file is a valid BED file.

        >>> a = pybedtools.example_bedtool('a.bed')
        >>> a.head(2) #doctest: +NORMALIZE_WHITESPACE
        chr1	1	100	feature1	0	+
        chr1	100	200	feature2	0	+
        <BLANKLINE>

        """
        if not isinstance(self.fn, str):
            raise NotImplementedError(
                "head() not supported for non file-based BedTools"
            )
        if as_string:
            return "".join(str(line) for line in self[:n])
        if self._isbam:
            raise NotImplementedError("head() not supported for BAM")
        else:
            if isGZIP(self.fn):
                openfunc = gzip.open
                openmode = "rt"
            else:
                openfunc = open
                openmode = "r"
            with openfunc(self.fn, openmode) as fin:
                for i, line in enumerate(fin):
                    if i == (n):
                        break
                    print(line, end=" ")

    def set_chromsizes(self, chromsizes:  str | dict):
        """
        Prepare BedTool for operations that require chromosome coords.

        Set the chromsizes for this genome. If *chromsizes* is a string, it
        will be considered a genome assembly name.  If that assembly name is
        not available in pybedtools.genome_registry, then it will be searched
        for on the UCSC Genome Browser.

        Example usage:

            >>> hg19 = pybedtools.chromsizes('hg19')
            >>> a = pybedtools.example_bedtool('a.bed')
            >>> a = a.set_chromsizes(hg19)
            >>> print(a.chromsizes['chr1'])
            (0, 249250621)

        """
        if isinstance(chromsizes, str):
            self.chromsizes = pybedtools.chromsizes(chromsizes)
        elif isinstance(chromsizes, dict):
            self.chromsizes = chromsizes
        else:
            raise ValueError(
                "Need to specify chromsizes either as a string"
                " (assembly name) or a dictionary"
            )
        return self

    def _collapse(
        self,
        iterable: Iterable,
        fn: Optional[str] = None,
        trackline: Optional[str] = None,
        in_compressed: bool = False,
        out_compressed: bool = False,
    ) -> str:
        """
        Collapses an iterable into file *fn* (or a new tempfile if *fn* is
        None).

        Returns the newly created filename.

        Parameters
        ----------

        iterable : iter
            Any iterable object whose items can be converted to an Interval.

        fn : str
            Output filename, if None then creates a temp file for output

        trackline : str
            If not None, string to be added to the top of the output. Newline
            will be added.

        in_compressed : bool
            Indicates whether the input is compressed

        out_compressed : bool
            Indicates whether the output should be compressed
        """
        if fn is None:
            fn = self._tmp()

        in_open_func = gzip.open if in_compressed else open
        out_open_func = gzip.open if out_compressed else open

        # special case: if BAM-format BedTool is provided, no trackline should
        # be supplied, and don't iterate -- copy the file wholesale
        if isinstance(iterable, BedTool) and iterable._isbam:
            if trackline:
                raise ValueError(
                    "trackline provided, but input is a BAM "
                    "file, which takes no track line"
                )
            with open(fn, "wb") as out_:
                out_.write(open(self.fn, "rb").read())
            return fn

        # If we're just working with filename-based BedTool objects, just copy
        # the files directly
        if isinstance(iterable, BedTool) and isinstance(iterable.fn, str):
            with out_open_func(fn, "wt") as out_:
                if sys.version_info > (3,0):
                    in_ = in_open_func(iterable.fn, "rt", errors="ignore")
                else:
                    in_ = in_open_func(iterable.fn, "rt")
                if trackline:
                    out_.write(trackline.strip() + "\n")
                out_.writelines(in_)
                in_.close()
        else:
            with out_open_func(fn, "wt") as out_:
                for i in iterable:
                    if isinstance(i, (list, tuple)):
                        i = create_interval_from_list(list(i))
                    out_.write(str(i))
        return fn

    def handle_kwargs(self, prog:str, arg_order: Optional[list[str]] = None, **kwargs):
        """
        Handle most cases of BEDTool program calls, but leave the specifics
        up to individual methods.

        *prog* is a BEDTools program name, e.g., 'intersectBed'.

        *arg_order* lists any arguments that are sensitive to order. Everything
        else will be reverse-sorted.

        *kwargs* are passed directly from the calling method (like
        self.intersect).

        This method figures out, given how this BedTool was constructed, what
        to send to BEDTools programs -- for example, an open file to stdin with
        the `-` argument, or a filename with the `-a` argument.
        """
        pybedtools.logger.debug(
            "BedTool.handle_kwargs() got these kwargs:\n%s", pprint.pformat(kwargs)
        )

        # If you pass in a list, how should it be converted to a BedTools arg?
        default_list_delimiter = " "
        list_delimiters = {
            "annotateBed": " ",
            "getOverlap": ",",
            "groupBy": ",",
            "multiIntersectBed": " ",
            "mergeBed": ",",
            "intersectBed": " ",
            "mapBed": ",",
        }
        stdin = None

        # -----------------------------------------------------------------
        # Decide how to send instream1 to BEDTools.  If there's no implicit
        # instream1 arg, then do nothing.
        #
        try:
            # e.g., 'a' for intersectBed
            if self._isbam:
                inarg1 = _bam_registry[prog]
            else:
                inarg1 = _implicit_registry[prog]

            # e.g., self.fn or 'a.bed' or an iterator...
            instream1 = kwargs[inarg1]

            # If it's a BedTool, then get underlying stream
            if isinstance(instream1, BedTool):
                instream1 = instream1.fn

            # Filename? No pipe, just provide the file
            if isinstance(instream1, str):
                kwargs[inarg1] = instream1
                stdin = None

            # Open file? Pipe it
            # elif isinstance(instream1, file):
            #     kwargs[inarg1] = 'stdin'
            #     stdin = instream1

            # A generator or iterator: pipe it as a generator of lines
            else:
                kwargs[inarg1] = "stdin"
                stdin = (str(i) for i in instream1)
        except KeyError:
            pass

        # -----------------------------------------------------------------
        # Decide how to send instream2 to BEDTools.
        try:
            # e.g., 'b' for intersectBed
            inarg2 = _other_registry[prog]

            # e.g., another BedTool
            instream2 = kwargs[inarg2]

            # Get stream if BedTool
            if isinstance(instream2, BedTool):
                instream2 = instream2.fn

            # Filename
            if isinstance(instream2, str):
                kwargs[inarg2] = instream2

            # If it's a list of strings, then we need to figure out if it's
            # a list of filenames or a list of intervals (see issue #156)
            #
            # Several options:
            #
            #   - assume intervals have tabs but filenames don't
            #   - assume that, upon being split on tabs, an interval is >=3 fields
            #   - try creating an interval out of the first thing, success means interval
            #
            # The last seems the most robust. It does allow filenames with
            # tabs; deciding whether or not such filenames are a good idea is
            # left to the user.
            #
            elif isinstance(instream2, (list, tuple)) and isinstance(
                instream2[0], str
            ):
                try:
                    _ = create_interval_from_list(instream2[0].split("\t"))
                    kwargs[inarg2] = self._collapse(instream2)
                except IndexError:
                    kwargs[inarg2] = instream2

            # Otherwise we need to collapse it in order to send to BEDTools
            # programs
            else:
                kwargs[inarg2] = self._collapse(instream2)

        except KeyError:
            pass

        # If stream not specified, then a tempfile will be created
        if kwargs.pop("stream", None):
            tmp = None
        else:
            output = kwargs.pop("output", None)
            if output:
                tmp = output
            else:
                tmp = self._tmp()

        additional_args = kwargs.pop("additional_args", None)

        # Parse the kwargs into BEDTools-ready args
        cmds = [prog]

        # arg_order mechanism added to fix #345
        if arg_order is None:
            arg_order = []

        for arg in arg_order:
            if arg in kwargs:
                val = kwargs.pop(arg)
                cmds.append("-" + arg)
                cmds.append(val)

        # The reverse-sort is a temp fix for issue #81
        for key, value in sorted(list(kwargs.items()), reverse=True):
            if isinstance(value, bool):
                if value:
                    cmds.append("-" + key)
                else:
                    continue
            elif isinstance(value, list) or isinstance(value, tuple):
                value = list(map(str, value))
                try:
                    delim = list_delimiters[prog]
                except KeyError:
                    delim = default_list_delimiter

                if delim == " ":
                    cmds.append("-" + key)
                    cmds.extend(value)

                # make comma-separated list if that's what's needed
                else:
                    cmds.append("-" + key)
                    cmds.append(delim.join(value))

            else:
                cmds.append("-" + key)
                cmds.append(str(value))

        if additional_args:
            cmds.append(additional_args)

        return cmds, tmp, stdin

    def check_genome(self, **kwargs):
        """
        Handles the different ways of specifying a genome in kwargs:

        g='genome.file' specifies a file directly
        genome='dm3' gets the file from genome registry
        self.chromsizes could be a dict.\
        """

        # If both g and genome are missing, assume self.chromsizes
        if ("g" not in kwargs) and ("genome" not in kwargs):
            if hasattr(self, "chromsizes"):
                kwargs["g"] = self.chromsizes
            else:
                raise ValueError(
                    'No genome specified. Use the "g" or '
                    '"genome" kwargs, or use the '
                    ".set_chromsizes() method"
                )

        # If both specified, rather than make an implicit decision, raise an
        # exception
        if "g" in kwargs and "genome" in kwargs:
            raise ValueError('Cannot specify both "g" and "genome"')

        # Something like genome='dm3' was specified
        if "g" not in kwargs and "genome" in kwargs:
            if isinstance(kwargs["genome"], dict):
                genome_dict = kwargs["genome"]
            else:
                genome_dict = pybedtools.chromsizes(kwargs["genome"])
            genome_file = pybedtools.chromsizes_to_file(genome_dict)
            kwargs["g"] = genome_file
            del kwargs["genome"]

        # By the time we get here, 'g' is specified.

        # If a dict was provided, convert to tempfile here
        if isinstance(kwargs["g"], dict):
            kwargs["g"] = pybedtools.chromsizes_to_file(kwargs["g"])

        if not os.path.exists(kwargs["g"]):
            msg = 'Genome file "%s" does not exist' % (kwargs["g"])
            raise FileNotFoundError(msg)

        return kwargs

    @_log_to_history
    def remove_invalid(self):
        """
        Remove invalid features that may break BEDTools programs.

        >>> a = pybedtools.BedTool("chr1 10 100\\nchr1 10 1",
        ... from_string=True)
        >>> print(a.remove_invalid()) #doctest: +NORMALIZE_WHITESPACE
        chr1	10	100
        <BLANKLINE>

        """
        tmp = self._tmp()

        # If it's a file-based BedTool -- which is likely, if we're trying to
        # remove invalid features -- then we need to parse it line by line.
        if isinstance(self.fn, str):
            i = IntervalIterator(open(self.fn, "r"))
        else:
            tmp = self.saveas()
            i = IntervalIterator(open(tmp.fn, "r"))

        def _generator():
            while True:
                try:
                    feature = next(i)
                    if feature.start <= feature.stop:
                        yield feature
                    else:
                        continue
                except pybedtools.MalformedBedLineError:
                    continue
                except OverflowError:
                    # This can happen if coords are negative
                    continue
                except IndexError:
                    continue
                except StopIteration:
                    break

        return BedTool(_generator())

    def all_hits(self, interval: Interval, same_strand: bool = False, overlap: float = 0.0):
        """
        Return all intervals that overlap `interval`.

        Calls the `all_hits` method of an IntervalFile to return all intervals
        in this current BedTool that overlap `interval`.

        Require that overlaps have the same strand with same_strand=True.

        Notes:
                If this current BedTool is generator-based, it will be
                converted into a file first.

                If this current BedTool refers to a BAM file, it will be
                converted to a BED file first using default arguments.  If you
                don't want this to happen, please convert to BED first before
                using this method.
        """
        if not isinstance(interval, Interval):
            raise ValueError("Need an Interval instance")
        fn = self.fn
        if not isinstance(fn, str):
            fn = self.saveas().fn
        if self._isbam:
            fn = self.bam_to_bed().fn
        interval_file = pybedtools.IntervalFile(fn)
        return interval_file.all_hits(interval, same_strand, overlap)

    def any_hits(self, interval: Interval, same_strand: bool = False, overlap: float=0.0):
        """
        Return whether or not any intervals overlap `interval`.

        Calls the `any_hits` method of an IntervalFile.  If there were any hits
        within `interval` in this BedTool, then return 1; otherwise 0.

        Require that overlaps have the same strand with same_strand=True.

        Notes:
                If this current BedTool is generator-based, it will be
                converted into a file first.

                If this current BedTool refers to a BAM file, it will be
                converted to a BED file first using default arguments.  If you
                don't want this to happen, please convert to BED first before
                using this method.
        """
        if not isinstance(interval, Interval):
            raise ValueError("Need an Interval instance")
        fn = self.fn
        if not isinstance(fn, str):
            fn = self.saveas().fn
        if self._isbam:
            fn = self.bam_to_bed().fn
        interval_file = pybedtools.IntervalFile(fn)
        return interval_file.any_hits(interval, same_strand, overlap)

    def count_hits(self, interval: Interval, same_strand: bool = False, overlap: float=0.0) -> int:
        """
        Return the number of intervals that overlap `interval`.

        Calls the `count_hits` method of an IntervalFile.  Returns the number
        of valid hits in this BedTool that overlap `interval`.

        Require that overlaps have the same strand with same_strand=True.

        Notes:
                If this current BedTool is generator-based, it will be
                converted into a file first.

                If this current BedTool refers to a BAM file, it will be
                converted to a BED file first using default arguments.  If you
                don't want this to happen, please convert to BED first before
                using this method.
        """
        if not isinstance(interval, Interval):
            raise ValueError("Need an Interval instance")
        fn = self.fn
        if not isinstance(fn, str):
            fn = self.saveas().fn
        if self._isbam:
            fn = self.bam_to_bed().fn
        interval_file = pybedtools.IntervalFile(fn)
        return interval_file.count_hits(interval, same_strand, overlap)

    @_log_to_history
    @_wraps(prog="bed12ToBed6", implicit="i", bam=None, other=None)
    def bed6(self, *args, **kwargs) -> BedTool: # type: ignore
        """
        Wraps `bedtools bed12tobed6`.
        """

    # Alias for backward compatibility
    bed12tobed6 = bed6

    @_log_to_history
    @_wraps(prog="bamToBed", implicit="i", other=None, nonbam="ALL", bam="i")
    def bam_to_bed(self, *args, **kwargs) -> BedTool: # type: ignore
        """
        Wraps `bedtools bamtobed`.
        """

    # Alias for backward compatibility
    bamtobed = bam_to_bed

    @_wraps(prog="bedToBam", implicit="i", uses_genome=True, force_bam=True)
    def _bed_to_bam(self, *args, **kwargs):
        """
        Wraps bedToBam and is called internally for BED/GFF/VCF files by
        self.to_bam (which needs to do something different for SAM files...)
        """

    @_log_to_history
    def to_bam(self, **kwargs):
        """
        Wraps `bedtools bedtobam`

        If self.fn is in BED/VCF/GFF format, call BEDTools' bedToBam.  If
        self.fn is in SAM format, then create a header out of the genome file
        and then convert using `samtools`.
        """
        if self.file_type == "bam":
            return self
        if self.file_type in ("bed", "gff", "vcf"):
            return self._bed_to_bam(**kwargs)

        # TODO: to maintain backwards compatibility we go from Interval to
        # AlignedSegment.
        if self.file_type == "sam":

            # Use pysam, but construct the header out of a provided genome
            # file.

            # construct a genome out of whatever kwargs were passed in
            kwargs = self.check_genome(**kwargs)

            # Build a header that we can use for the output BAM file.
            genome = dict(i.split() for i in open(kwargs["g"]))
            SQ = []
            ref_ids = {}
            text_header = ["@HD\tVN:1.0"]

            for i, (k, v) in enumerate(genome.items()):
                SQ.append(dict(SN=k, LN=int(v)))
                ref_ids[k] = i
                text_header.append("@SQ\tSN:{0}\tLN:{1}".format(k, v))

            # And the text-format header
            text_header = "\n".join(text_header) + "\n"

            # The strategy is to write an actual SAM file to disk, along with
            # a header, and then read that back in.
            #
            # Painfully inefficient, but this will change once all py2 tests
            # pass.
            sam_tmp = self._tmp()
            bam_tmp = self._tmp()
            with open(sam_tmp, "w") as fout:
                fout.write(text_header)
                for interval in self:
                    fout.write("\t".join(map(str, interval.fields)) + "\n")

            samfile = pysam.AlignmentFile(sam_tmp, "r")
            bamfile = pysam.AlignmentFile(bam_tmp, "wb", template=samfile)
            for alignment in samfile:
                bamfile.write(alignment)

            samfile.close()
            bamfile.close()
            new_bedtool = BedTool(bam_tmp)
            new_bedtool._isbam = True
            return new_bedtool

    # Alias for backward compatibility
    bedtobam = to_bam

    @_log_to_history
    @_wraps(prog="intersectBed", implicit="a", other="b", bam="abam",
            nonbam="bed", arg_order=["a", "abam"])
    def intersect(self, *args, **kwargs) -> BedTool:  # type: ignore
        """
        Wraps `bedtools intersect`.
        """

    @_log_to_history
    @_wraps(
        prog="fastaFromBed",
        implicit="bed",
        bam=None,
        other="fi",
        make_tempfile_for="fo",
        check_stderr=_check_sequence_stderr,
        add_to_bedtool={"fo": "seqfn"},
    )
    def sequence(self, *args, **kwargs) -> BedTool:  # type: ignore
        '''
        Wraps `bedtools getfasta`.

        *fi* is passed in by the user; *bed* is automatically passed in as the
        bedfile of this object; *fo* by default is a temp file.  Use
        save_seqs() to save as a file.

        The end result is that this BedTool will assign a value to the attribute , self.seqfn,
        that points to the new fasta file.

        Example usage:

        >>> a = pybedtools.BedTool("""
        ... chr1 1 10
        ... chr1 50 55""", from_string=True)
        >>> fasta = pybedtools.example_filename('test.fa')
        >>> a = a.sequence(fi=fasta)
        >>> print(open(a.seqfn).read())
        >chr1:1-10
        GATGAGTCT
        >chr1:50-55
        CCATC
        <BLANKLINE>

        '''

    # Alias for backwards compatibility
    getfasta = sequence

    @staticmethod
    def seq(loc, fasta) -> str:
        """
        Return just the sequence from a region string or a single location
        >>> fn = pybedtools.example_filename('test.fa')
        >>> BedTool.seq('chr1:2-10', fn)
        'GATGAGTCT'
        >>> BedTool.seq(('chr1', 1, 10), fn)
        'GATGAGTCT'
        """
        if isinstance(loc, str):
            chrom, start_end = loc.split(":")
            start, end = list(map(int, start_end.split("-")))
            start -= 1
        else:
            chrom, start, end = loc[0], loc[1], loc[2]

        loc = BedTool("%s\t%i\t%i" % (chrom, start, end), from_string=True)
        lseq = loc.sequence(fi=fasta)
        return "".join([l.rstrip() for l in open(lseq.seqfn, "r") if l[0] != ">"])

    @_log_to_history
    @_wraps(
        prog="nucBed", implicit="bed", other="fi", check_stderr=_check_sequence_stderr
    )
    def nucleotide_content(self) -> BedTool:  # type: ignore
        """
        Wraps `bedtools nuc`.

        Profiles nucleotide content.  The returned BED file contains extra
        information about the nucleotide content
        """

    # Alias for backwards compatibility
    nuc = nucleotide_content

    @_log_to_history
    @_wraps(prog="multiBamCov", implicit="bed")
    def multi_bam_coverage(self) -> BedTool:  # type: ignore
        """
        Wraps `bedtools multicov`.

        Pass a list of sorted and indexed BAM files as `bams`
        """

    # Alias for backwards compatibility
    multicov = multi_bam_coverage

    @_log_to_history
    @_wraps(prog="subtractBed", implicit="a", other="b", bam=None)
    def subtract(self, *args, **kwargs) -> BedTool: # type: ignore
        """
        Wraps `bedtools subtract`.

        Subtracts from another BED file and returns a new BedTool object.

        Example usage:

            >>> a = pybedtools.example_bedtool('a.bed')
            >>> b = pybedtools.example_bedtool('b.bed')

            Do a "stranded" subtraction:

            >>> c = a.subtract(b, s=True)

            Require 50% of features in `a` to overlap:

            >>> c = a.subtract(b, f=0.5)
        """
        if "a" not in kwargs:
            kwargs["a"] = self.fn

        if "b" not in kwargs:
            if len(args) > 0:
                kwargs["b"] = args[0]
            else:
                raise ValueError("Must specify a BED file to subtract, either as a positional argument or as the 'b' keyword argument.")

        cmds, tmp, stdin = self.handle_kwargs(prog="subtractBed", **kwargs)
        stream = call_bedtools(cmds, tmp, stdin=stdin)
        return BedTool(stream)

    @_log_to_history
    @_wraps(prog="slopBed", implicit="i", other=None, bam=None, uses_genome=True)
    def slop(self, *args, **kwargs) -> BedTool: # type: ignore
        """
        Wraps `bedtools slop`.
        """

    @_log_to_history
    @_wraps(prog="shiftBed", implicit="i", other=None, bam=None, uses_genome=True)
    def shift(self, *args, **kwargs) -> BedTool: # type: ignore
        """
        Wraps `bedtools shift`.

        Shift each feature by user-defined number of bases. Returns a new BedTool object.

        Example usage:

            >>> a = pybedtools.example_bedtool('a.bed')

            Shift every feature by 5bp:

            >>> b = a.shift(genome='hg19', s=5)
            >>> print(b) #doctest: +NORMALIZE_WHITESPACE
            chr1	6	105	feature1	0	+
            chr1	105	205	feature2	0	+
            chr1	155	505	feature3	0	-
            chr1	905	955	feature4	0	+
            <BLANKLINE>

            Shift features on the '+' strand by -1bp and on '-' strand by +3bp:

            >>> b = a.shift(genome='hg19', p=-1, m=3)
            >>> print(b) #doctest: +NORMALIZE_WHITESPACE
            chr1	0	99	feature1	0	+
            chr1	99	199	feature2	0	+
            chr1	153	503	feature3	0	-
            chr1	899	949	feature4	0	+
            <BLANKLINE>

            # Disabling, see https://github.com/arq5x/bedtools2/issues/807
            Shift features by a fraction of their length (0.50):

            #>>> b = a.shift(genome='hg19', pct=True, s=0.50)
            #>>> print(b) #doctest: +NORMALIZE_WHITESPACE
            #chr1	50	149	feature1	0	+
            #chr1	150	250	feature2	0	+
            #chr1	325	675	feature3	0	-
            #chr1	925	975	feature4	0	+
            #<BLANKLINE>

        """

    @_log_to_history
    @_wraps(prog="mergeBed", implicit="i", other=None, bam=None)
    def merge(self, *args, **kwargs) -> BedTool: # type: ignore
        """
        Wraps `bedtools merge`.

        Merge overlapping features together. Returns a new BedTool object.

        Example usage:

            >>> a = pybedtools.example_bedtool('a.bed')

            Merge:

            >>> c = a.merge()

            Allow merging of features 500 bp apart:

            >>> c = a.merge(d=500)

        """

    @_log_to_history
    @_wraps(prog="closestBed", implicit="a", other="b", bam=None)
    def closest(self, *args, **kwargs) -> BedTool: # type: ignore
        """
        Wraps `bedtools closest`.

        Return a new BedTool object containing closest features in *b*.  Note
        that the resulting file is no longer a valid BED format; use the
        special "_closest" methods to work with the resulting file.

        Example usage::

            a = BedTool('in.bed')

            # get the closest feature in 'other.bed' on the same strand
            b = a.closest('other.bed', s=True)

        """

    @_log_to_history
    @_wraps(prog="windowBed", implicit="a", other="b", bam="abam", nonbam="bed")
    def window(self, *args, **kwargs) -> BedTool: # type: ignore
        """
        Wraps `bedtools window`.

        Example usage::

            >>> a = pybedtools.example_bedtool('a.bed')
            >>> b = pybedtools.example_bedtool('b.bed')
            >>> print(a.window(b, w=1000)) #doctest: +NORMALIZE_WHITESPACE
            chr1	1	100	feature1	0	+	chr1	155	200	feature5	0	-
            chr1	1	100	feature1	0	+	chr1	800	901	feature6	0	+
            chr1	100	200	feature2	0	+	chr1	155	200	feature5	0	-
            chr1	100	200	feature2	0	+	chr1	800	901	feature6	0	+
            chr1	150	500	feature3	0	-	chr1	155	200	feature5	0	-
            chr1	150	500	feature3	0	-	chr1	800	901	feature6	0	+
            chr1	900	950	feature4	0	+	chr1	155	200	feature5	0	-
            chr1	900	950	feature4	0	+	chr1	800	901	feature6	0	+
            <BLANKLINE>
        """

    @_log_to_history
    @_wraps(prog="shuffleBed", implicit="i", other=None, bam=None, uses_genome=True)
    def shuffle(self, *args, **kwargs) -> BedTool: # type: ignore
        """
        Wraps `bedtools shuffle`.

        Example usage:

        >>> a = pybedtools.example_bedtool('a.bed')
        >>> seed = 1 # so this test always returns the same results
        >>> b = a.shuffle(genome='hg19', chrom=True, seed=seed)
        >>> print(b) #doctest: +NORMALIZE_WHITESPACE
        chr1	123081365	123081464	feature1	0	+
        chr1	243444570	243444670	feature2	0	+
        chr1	194620241	194620591	feature3	0	-
        chr1	172792873	172792923	feature4	0	+
        <BLANKLINE>
        """

    @_log_to_history
    @_wraps(prog="sortBed", implicit="i", uses_genome=True, genome_if=["g", "genome"])
    def sort(self, *args, **kwargs) -> BedTool: # type: ignore
        """
        Wraps `bedtools sort`.

        Note that chromosomes are sorted lexicographically, so chr12 will come
        before chr9.

        Example usage:

        >>> a = pybedtools.BedTool('''
        ... chr9 300 400
        ... chr1 100 200
        ... chr1 1 50
        ... chr12 1 100
        ... chr9 500 600
        ... ''', from_string=True)
        >>> print(a.sort()) #doctest: +NORMALIZE_WHITESPACE
        chr1	1	50
        chr1	100	200
        chr12	1	100
        chr9	300	400
        chr9	500	600
        <BLANKLINE>
        """

    @_log_to_history
    @_wraps(prog="annotateBed", implicit="i")
    def annotate(self, *args, **kwargs) -> BedTool: # type: ignore
        """
        Wraps  `bedtools annotate`.

        Annotate this BedTool with a list of other files.
        Example usage:

        >>> a = pybedtools.example_bedtool('a.bed')
        >>> b_fn = pybedtools.example_filename('b.bed')
        >>> print(a.annotate(files=b_fn)) #doctest: +NORMALIZE_WHITESPACE
        chr1	1	100	feature1	0	+	0.000000
        chr1	100	200	feature2	0	+	0.450000
        chr1	150	500	feature3	0	-	0.128571
        chr1	900	950	feature4	0	+	0.020000
        <BLANKLINE>
        """

    @_log_to_history
    @_wraps(prog="flankBed", implicit="i", uses_genome=True)
    def flank(self, *args, **kwargs) -> BedTool:
        """
        Wraps `bedtools flank`.

        Example usage:

        >>> a = pybedtools.example_bedtool('a.bed')
        >>> print(a.flank(genome='hg19', b=100)) #doctest: +NORMALIZE_WHITESPACE
        chr1	0	1	feature1	0	+
        chr1	100	200	feature1	0	+
        chr1	0	100	feature2	0	+
        chr1	200	300	feature2	0	+
        chr1	50	150	feature3	0	-
        chr1	500	600	feature3	0	-
        chr1	800	900	feature4	0	+
        chr1	950	1050	feature4	0	+
        <BLANKLINE>

        """
        kwargs = self.check_genome(**kwargs)

        if "i" not in kwargs:
            kwargs["i"] = self.fn

        cmds, tmp, stdin = self.handle_kwargs(prog="flankBed", **kwargs)
        stream = call_bedtools(cmds, tmp, stdin=stdin)
        return BedTool(stream)

    @_log_to_history
    @_wraps(
        prog="genomeCoverageBed",
        implicit="i",
        bam="ibam",
        genome_none_if=["ibam"],
        genome_ok_if=["ibam"],
        uses_genome=True,
        nonbam="ALL",
    )
    def genome_coverage(self, *args, **kwargs):
        """
        Wraps `bedtools genomecov`.

        Note that some invocations of `bedtools genomecov` do not result in
        a properly-formatted BED file. For example, the default behavior is to
        report a histogram of coverage. Iterating over the resulting,
        non-BED-format file will raise exceptions in pybedtools' parser.

        Consider using the `BedTool.to_dataframe` method to convert these
        non-BED files into a pandas DataFrame for further use.

        Example usage:

        BAM file input does not require a genome:

        >>> a = pybedtools.example_bedtool('x.bam')
        >>> b = a.genome_coverage(bg=True)
        >>> b.head(3) #doctest: +NORMALIZE_WHITESPACE
        chr2L	9329	9365	1
        chr2L	10212	10248	1
        chr2L	10255	10291	1

        Other input does require a genome:

        >>> a = pybedtools.example_bedtool('x.bed')
        >>> b = a.genome_coverage(bg=True, genome='dm3')
        >>> b.head(3) #doctest: +NORMALIZE_WHITESPACE
        chr2L	9329	9365	1
        chr2L	10212	10248	1
        chr2L	10255	10291	1

        Non-BED format results:
        >>> a = pybedtools.example_bedtool('x.bed')
        >>> b = a.genome_coverage(genome='dm3')
        >>> df = b.to_dataframe(names=['chrom', 'depth', 'n', 'chromsize', 'fraction'])
        """

    # Alias for backwards compatibility
    genomecov = genome_coverage

    @_log_to_history
    @_wraps(prog="coverageBed", implicit="a", other="b", bam="abam", nonbam="ALL")
    def coverage(self, *args, **kwargs) -> BedTool: # type: ignore
        """
        Wraps `bedtools coverage`.

        Note that starting in version 2.24.0, BEDTools swapped the semantics of
        the "a" and "b" files.

        Example usage:

        >>> a = pybedtools.example_bedtool('a.bed')
        >>> b = pybedtools.example_bedtool('b.bed')
        >>> c = b.coverage(a)
        >>> c.head(3) #doctest: +NORMALIZE_WHITESPACE
        chr1	155	200	feature5	0	-	2	45	45	1.0000000
        chr1	800	901	feature6	0	+	1	1	101	0.0099010
        """

    @_log_to_history
    @_wraps(
        prog="maskFastaFromBed",
        implicit="bed",
        other="fi",
        make_tempfile_for="fo",
        add_to_bedtool={"fo": "seqfn"},
        check_stderr=_check_sequence_stderr,
    )
    def mask_fasta(self, *args, **kwargs) -> BedTool: # type: ignore
        """
        Wraps `bedtools maskfasta`.

        Masks a fasta file at the positions in a BED file and saves result as
        'out' and stores the filename in seqfn.

        >>> a = pybedtools.BedTool('chr1 100 110', from_string=True)
        >>> fasta_fn = pybedtools.example_filename('test.fa')
        >>> a = a.mask_fasta(fi=fasta_fn, fo='masked.fa.example')
        >>> b = a.slop(b=2, genome='hg19')
        >>> b = b.sequence(fi=a.seqfn)
        >>> print(open(b.seqfn).read())
        >chr1:98-112
        TTNNNNNNNNNNAT
        <BLANKLINE>
        >>> os.unlink('masked.fa.example')
        >>> if os.path.exists('masked.fa.example.fai'):
        ...     os.unlink('masked.fa.example.fai')
        """

    # Alias for backwards compatibility
    maskfasta = mask_fasta

    @_log_to_history
    @_wraps(prog="complementBed", implicit="i", uses_genome=True)
    def complement(self, *args, **kwargs) -> BedTool: # type: ignore
        """
        Wraps `bedtools complement`.
        Example usage:

        >>> a = pybedtools.example_bedtool('a.bed')
        >>> a.complement(genome='hg19').head(5) #doctest: +NORMALIZE_WHITESPACE
        chr1	0	1
        chr1	500	900
        chr1	950	249250621
        chr2	0	243199373
        chr3	0	198022430
        """

    @_log_to_history
    @_wraps(prog="getOverlap", implicit="i")
    def overlap(self, *args, **kwargs) -> BedTool: # type: ignore
        """
        Wraps `bedtools overlap`.

        Example usage:

        >>> a = pybedtools.example_bedtool('a.bed')
        >>> b = pybedtools.example_bedtool('b.bed')
        >>> c = a.window(b, w=10).overlap(cols=[2,3,8,9])
        >>> print(c) #doctest: +NORMALIZE_WHITESPACE
        chr1	100	200	feature2	0	+	chr1	155	200	feature5	0	-	45
        chr1	150	500	feature3	0	-	chr1	155	200	feature5	0	-	45
        chr1	900	950	feature4	0	+	chr1	800	901	feature6	0	+	1
        <BLANKLINE>
        """

    # TODO: needs test files and doctests written
    @_log_to_history
    @_wraps(prog="pairToBed", implicit="a", other="b", bam="abam", nonbam="bedpe")
    def pair_to_bed(self, *args, **kwargs) -> BedTool: # type: ignore
        """
        Wraps `bedtools pairtobed`.
        """

    # Alias for backwards compatibility
    pairtobed = pair_to_bed

    @_log_to_history
    @_wraps(prog="pairToPair", implicit="a", other="b")
    def pair_to_pair(self) -> BedTool: # type: ignore
        """
        Wraps `bedtools pairtopair`.
        """

    # Alias for backwards compatibility
    pairtopair = pair_to_pair

    @_log_to_history
    @_wraps(prog="groupBy", implicit="i")
    def groupby(self, *args, **kwargs) -> BedTool:
        """
        Wraps `bedtools groupby`.

        Example usage:

        >>> a = pybedtools.example_bedtool('gdc.gff')
        >>> b = pybedtools.example_bedtool('gdc.bed')
        >>> c = a.intersect(b, c=True)
        >>> d = c.groupby(g=[1, 4, 5], c=10, o=['sum'])
        >>> print(d) #doctest: +NORMALIZE_WHITESPACE
        chr2L	41	70	0
        chr2L	71	130	2
        chr2L	131	170	4
        chr2L	171	200	0
        chr2L	201	220	1
        chr2L	41	130	2
        chr2L	171	220	1
        chr2L	41	220	7
        chr2L	161	230	6
        chr2L	41	220	7
        <BLANKLINE>

        """

    @_log_to_history
    @_wraps(prog="tagBam", implicit="i", bam="i")
    def tag_bam(self, *args, **kwargs) -> pysam.AlignmentFile: # type: ignore
        """
        Wraps `bedtools tag`.

        `files` and `labels` should lists of equal length.

        """
        if "i" in kwargs:
            if "labels" in kwargs:
                if len(kwargs["i"]) != len(kwargs["labels"]):
                    raise ValueError("files and labels must be lists of equal length")


    # Alias for backwards compatibility
    tag = tag_bam

    @_log_to_history
    @_wraps(prog="mapBed", implicit="a", other="b")
    def map(self, *args, **kwargs) -> BedTool: # type: ignore
        """
        Wraps  `bedtools map`; See also :meth:`BedTool.each`.
        """

    @_log_to_history
    @_wraps(prog="multiIntersectBed", uses_genome=True, genome_if=["empty"])
    def multi_intersect(self, *args, **kwargs) -> BedTool: # type: ignore
        """
        Wraps `bedtools multiintersect`.

        Provide a list of filenames as the "i" argument. e.g. if you already
        have BedTool objects then use their `.fn` attribute, like this::

            >>> x = pybedtools.BedTool()
            >>> a = pybedtools.example_bedtool('a.bed')
            >>> b = pybedtools.example_bedtool('b.bed')
            >>> result = x.multi_intersect(i=[a.fn, b.fn])
            >>> print(result)   #doctest: +NORMALIZE_WHITESPACE
            chr1	1	155	1	1	1	0
            chr1	155	200	2	1,2	1	1
            chr1	200	500	1	1	1	0
            chr1	800	900	1	2	0	1
            chr1	900	901	2	1,2	1	1
            chr1	901	950	1	1	1	0
            <BLANKLINE>

        """

    # Alias for backwards compatibility
    multiinter = multi_intersect

    @_log_to_history
    @_wraps(prog="randomBed", uses_genome=True)
    def random(self, *args, **kwargs) -> BedTool: # type: ignore
        """
        Wraps `bedtools random`.

        Since this method does not operate on an existing file, create
        a BedTool with no arguments and then call this method, e.g.,

        >>> x = BedTool()
        >>> y = x.random(l=100, n=10, genome='hg19')
        """

    @_log_to_history
    @_wraps("bedpeToBam", implicit="i", uses_genome=True, force_bam=True)
    def bedpe_to_bam(self, *args, **kwargs) -> pysam.AlignmentFile: # type: ignore
        """
        Wraps `bedtools bedpetobam`.
        """

    # Alias for backwards compatibility
    bedpetobam = bedpe_to_bam

    @_log_to_history
    @_wraps(prog="clusterBed", implicit="i")
    def cluster(self, *args, **kwargs) -> BedTool: # type: ignore
        """
        Wraps  `bedtools cluster`.
        """

    @_log_to_history
    @_wraps(prog="unionBedGraphs")
    def union_bedgraphs(self, *args, **kwargs) -> BedTool: # type: ignore
        """
        Wraps `bedtools unionbedg`.

        Warning: using the `header=True` kwarg will result in a file that is
        not in true BED format, which may break downstream analysis.
        """
        if "header" in kwargs:
            if kwargs["header"] is True:
                warn("Using header=True with unionbedg will result in a file that is not in true BED format, which may break downstream analysis.")


    # Alias for backwards compatibility
    unionbedg = union_bedgraphs

    @_log_to_history
    @_wraps(prog="windowMaker", uses_genome=True, genome_none_if=["b"], other="b", arg_order=["w"])
    def window_maker(self, *args, **kwargs) -> BedTool: # type: ignore
        """
        Wraps `bedtools makewindows`.
        """

    # Alias for backwards compatibility
    makewindows = window_maker

    @_log_to_history
    @_wraps(prog="expandCols", implicit="i")
    def expand(self, *args, **kwargs) -> BedTool: # type: ignore
        """
        Wraps `bedtools expand`
        """

    @_log_to_history
    @_wraps(prog="linksBed", implicit="i", add_to_bedtool={"stdout": "links_html"})
    def links(self, *args, **kwargs) -> BedTool: # type: ignore
        """
        Wraps `linksBed`.

        The resulting BedTool will assign a value to the attribute `links_html`.  This
        attribute is a temp filename containing the HTML links.
        """

    @_log_to_history
    @_wraps(prog="bedToIgv", implicit="i", add_to_bedtool={"stdout": "igv_script"})
    def igv(self, *args, **kwargs) -> BedTool: # type: ignore
        """
        Wraps `bedtools igv`.

        The resulting BedTool will assign a value to the attribute `igv_script`.  This
        attribute is a temp filename containing the IGV script.
        """

    @_log_to_history
    @_wraps(
        prog="bamToFastq",
        implicit="i",
        bam="i",
        make_tempfile_for="fq",
        add_to_bedtool={"fq": "fastq"},
    )
    def bam_to_fastq(self, *args, **kwargs) -> BedTool: # type: ignore
        """
        Wraps `bedtools bamtofastq`.

        The `fq` argument is required.

        The resulting BedTool will assign a value to the attribute `fastq`.
        """

    # Alias for backwards compatibility
    bamtofastq = bam_to_fastq

    @_wraps(
        prog="jaccard",
        implicit="a",
        other="b",
        does_not_return_bedtool=_jaccard_output_to_dict,
    )
    def jaccard(self, *args, **kwargs) -> dict[str, Any]: # type: ignore
        """
        Returns a dictionary with keys (intersection, union, jaccard).
        """

    @_wraps(
        prog="reldist",
        implicit="a",
        other="b",
        does_not_return_bedtool=_reldist_output_handler,
    )
    def reldist(self, *args, **kwargs) -> BedTool | dict[str, Any]: # type: ignore
        """
        If detail=False, then return a dictionary with keys (reldist, count,
        total, fraction), which is the summary of the bedtools reldist.

        Otherwise return a BedTool, with the relative distance for each
        interval in A in the last column.
        """
        if "detail" in kwargs:
            if kwargs["detail"] is False:
                warn("Using detail=False with reldist will return a dictionary with keys (reldist, count, total, fraction), which is the summary of the bedtools reldist."
                     "Not a BedTool object.")


    @_wraps(prog="sample", implicit="i", bam="i")
    def sample(self, *args, **kwargs) -> BedTool: # type: ignore
        """
        Wraps 'sample'.
        """

    @_wraps(
        prog="fisher",
        implicit="a",
        other="b",
        uses_genome=True,
        does_not_return_bedtool=FisherOutput,
    )
    def fisher(self, *args, **kwargs) -> FisherOutput: # type: ignore
        """
        Wraps 'fisher'. Returns an object representing the output.

        >>> a = pybedtools.example_bedtool('a.bed')
        >>> b = pybedtools.example_bedtool('b.bed')
        >>> f = a.fisher(b, genome='hg19')
        >>> print(f)  # doctest: +NORMALIZE_WHITESPACE
        # Number of query intervals: 4
        # Number of db intervals: 2
        # Number of overlaps: 3
        # Number of possible intervals (estimated): 13958448
        # phyper(3 - 1, 4, 13958448 - 4, 2, lower.tail=F)
        # Contingency Table Of Counts
        #_________________________________________
        #           |  in -b       | not in -b    |
        #     in -a | 3            | 1            |
        # not in -a | 0            | 13958444     |
        #_________________________________________
        # p-values for fisher's exact test
        left	right	two-tail	ratio
        1	8.8247e-21	8.8247e-21	inf
        <BLANKLINE>


        >>> f.table['not in -a']['in -b']
        0

        >>> f.table['not in -a']['not in -b']
        13958444

        >>> f.table['in -a']['in -b']
        3

        >>> f.table['in -a']['not in -b']
        1

        >>> f.two_tail
        8.8247e-21
        """

    @_wraps(prog="split", implicit="i", does_not_return_bedtool=SplitOutput)
    def splitbed(self, *args, **kwargs) -> SplitOutput:  # type: ignore
        """
        Wraps 'bedtools split'.

        BedTool objects have long had a `split` method which splits intervals
        according to a custom function. Now that BEDTools has a `split` tool,
        the method name conflicts. To maintain backwards compatibility, the
        method wrapping the BEDTools command is called `splitbed`.

        Since this tool does not return a single BED file, the method parses
        the output and returns a SplitOutput object, which includes an
        attribute, `bedtools`, that is a list of BedTool objects created from
        the split files.

        To keep the working directory clean, you may want to consider using
        `prefix=BedTool._tmp()` to get a temp file that will be deleted when
        Python exits cleanly.

        >>> a = pybedtools.example_bedtool('a.bed')
        >>> s = a.splitbed(n=2, p="split")
        >>> assert len(a) == 4, len(a)
        >>> assert len(s.bedtools) == 2
        >>> print(s.bedtools[0]) # doctest: +NORMALIZE_WHITESPACE
        chr1	150	500	feature3	0	-
        <BLANKLINE>
        >>> print(s.bedtools[1]) # doctest: +NORMALIZE_WHITESPACE
        chr1	100	200	feature2	0	+
        chr1	1	100	feature1	0	+
        chr1	900	950	feature4	0	+
        <BLANKLINE>
        """

    @_wraps(prog="spacing", implicit="i")
    def spacing(self, *args, **kwargs) -> BedTool:  # type: ignore
        """
        Wraps `bedtools spacing`

        >>> a = pybedtools.example_bedtool('a.bed')
        >>> print(a.spacing())  # doctest: +NORMALIZE_WHITESPACE
        chr1	1	100	feature1	0	+	.
        chr1	100	200	feature2	0	+	0
        chr1	150	500	feature3	0	-	-1
        chr1	900	950	feature4	0	+	400
        """

    def count(self) -> int:
        """
        Count the number features in this BedTool.

        Number of features in BED file. Does the same thing as len(self), which
        actually just calls this method.

        Only counts the actual features.  Ignores any track lines, browser
        lines, lines starting with a "#", or blank lines.

        Example usage:

        >>> a = pybedtools.example_bedtool('a.bed')
        >>> a.count()
        4
        """
        if hasattr(self, "next") or hasattr(self, "__next__"):
            return sum(1 for _ in self)
        return sum(1 for _ in iter(self))

    def print_sequence(self) -> str:
        """
        Print the sequence that was retrieved by BedTool.sequence.
        See usage example in :meth:`BedTool.sequence`.

        Returns:
            str: The sequence as a string.

        Raises:
            ValueError: If the sequence has not been generated using .sequence(fasta).
        """
        if self.seqfn is None:
            raise ValueError("Use .sequence(fasta) to get the sequence first")
        f = open(self.seqfn)
        s = f.read()
        f.close()
        return s

    def save_seqs(self, fn:str) -> BedTool:
        """
        Save sequences, after calling BedTool.sequence.

        In order to use this function, you need to have called
        the :meth:`BedTool.sequence()` method.

        A new BedTool object is returned which references the newly saved file.

        Example usage:

        >>> a = pybedtools.BedTool('''
        ... chr1 1 10
        ... chr1 50 55''', from_string=True)
        >>> fasta = pybedtools.example_filename('test.fa')
        >>> a = a.sequence(fi=fasta)
        >>> print(open(a.seqfn).read())
        >chr1:1-10
        GATGAGTCT
        >chr1:50-55
        CCATC
        <BLANKLINE>
        >>> b = a.save_seqs('example.fa')
        >>> assert open(b.fn).read() == open(a.fn).read()
        >>> if os.path.exists('example.fa'):
        ...     os.unlink('example.fa')
        """

        if self.seqfn is None:
            raise ValueError("Use .sequence(fasta) to get the sequence first")

        with open(fn, "w") as fout:
            if self.seqfn is None:
                raise ValueError("Use .sequence(fasta) to get the sequence first")
            with open(self.seqfn) as seqfile:
                fout.write(seqfile.read())

        new_bedtool = BedTool(self.fn)
        new_bedtool.seqfn = fn
        return new_bedtool

    def randomstats(
        self,
        other: BedTool,
        iterations: int,
        new: bool = False,
        genome_fn: Optional[str] = None,
        include_distribution: bool = False,
        **kwargs
    ) -> dict[str, Any]:
        """
        Dictionary of results from many randomly shuffled intersections.

        Sends args and kwargs to :meth:`BedTool.randomintersection` and
        compiles results into a dictionary with useful stats.  Requires
        numpy.

        If `include_distribution` is True, then the dictionary will include the
        full distribution; otherwise, the distribution is deleted and cleaned
        up to save on memory usage.

        This is one possible way of assigning significance to overlaps between
        two files. See, for example:

            Negre N, Brown CD, Shah PK, Kheradpour P, Morrison CA, et al. 2010
            A Comprehensive Map of Insulator Elements for the Drosophila
            Genome. PLoS Genet 6(1): e1000814. doi:10.1371/journal.pgen.1000814

        Example usage:

        Make chromsizes a very small genome for this example:

        >>> chromsizes = {'chr1':(1,1000)}
        >>> a = pybedtools.example_bedtool('a.bed').set_chromsizes(chromsizes)
        >>> b = pybedtools.example_bedtool('b.bed')
        >>> try:
        ...     results = a.randomstats(b, 100, debug=True)
        ... except ImportError:
        ...     pass

        *results* is a dictionary that you can inspect.

        (Note that the following examples are not run as part of the doctests
        to avoid forcing users to install NumPy just to pass tests)

        The actual overlap::

            print(results['actual'])
            3

        The median of all randomized overlaps::

            print(results['median randomized'])
            2.0

        The percentile of the actual overlap in the distribution of randomized
        overlaps, which can be used to get an empirical p-value::

            print(results['percentile'])
            90.0
        """
        if ("intersect_kwargs" not in kwargs) or (kwargs["intersect_kwargs"] is None):
            kwargs["intersect_kwargs"] = {"u": True}
        try:
            import numpy as np
        except ImportError:
            raise ImportError("Need to install NumPy for stats...")

        def percentileofscore(a, score):
            """
            copied from scipy.stats.percentileofscore, to avoid dependency on
            scipy.
            """
            a = np.array(a)
            n = len(a)

            if not (np.any(a == score)):
                a = np.append(a, score)
                a_len = np.array(list(range(len(a))))
            else:
                a_len = np.array(list(range(len(a)))) + 1.0

            a = np.sort(a)
            idx = tuple([a == score])
            pct = (np.mean(a_len[idx]) / n) * 100.0
            return pct

        if isinstance(other, str):
            other = BedTool(other)
        else:
            assert isinstance(
                other, BedTool
            ), "Either filename or another BedTool instance required"

        # Actual (unshuffled) counts.
        i_kwargs = kwargs["intersect_kwargs"]
        actual = len(self.intersect(other, **i_kwargs))

        # List of counts from randomly shuffled versions.
        # Length of counts == *iterations*.

        if not new:
            distribution = self.randomintersection(
                other, iterations=iterations, **kwargs
            )
        else:
            # use new mechanism
            if genome_fn is None:
                raise ValueError(
                    "`genome_fn` must be provided if using the "
                    "new _randomintersection mechanism"
                )
            distribution = self._randomintersection(
                other, iterations=iterations, genome_fn=genome_fn, **kwargs
            )

        distribution = np.array(list(distribution))

        # Median of distribution
        med_count = np.median(distribution)

        n = float(len(distribution))

        frac_above = sum(distribution > actual) / n
        frac_below = sum(distribution < actual) / n

        normalized = actual / med_count

        lower_thresh = 2.5
        upper_thresh = 97.5
        lower, upper = np.percentile(distribution, [lower_thresh, upper_thresh])

        actual_percentile = percentileofscore(distribution, actual)
        d = {
            "iterations": iterations,
            "actual": actual,
            "file_a": self.fn,
            "file_b": other.fn,
            self.fn: len(self),
            other.fn: len(other),
            "self": len(self),
            "other": len(other),
            "frac randomized above actual": frac_above,
            "frac randomized below actual": frac_below,
            "median randomized": med_count,
            "normalized": normalized,
            "percentile": actual_percentile,
            "lower_%sth" % lower_thresh: lower,
            "upper_%sth" % upper_thresh: upper,
        }
        if include_distribution:
            d["distribution"] = distribution
        else:
            del distribution
        return d

    def random_op(self, *args, **kwargs):
        """
        For backwards compatibility; see BedTool.parallel_apply instead.
        """
        return self.parallel_apply(*args, **kwargs)

    def parallel_apply(
        self,
        iterations: int,
        func: Callable,
        func_args: Iterable,
        func_kwargs: dict,
        processes: int=1,
        _orig_pool: Optional[Pool] = None  # type: ignore
    ):
        """
        Generalized method for applying a function in parallel.

        Typically used when having to do many random shufflings.

        `func_args` and `func_kwargs` will be passed to `func` each time in
        `iterations`, and these iterations will be split across `processes`
        processes.

        Notes on the function, `func`:

            * the function should manually remove any tempfiles created.  This
              is because the BedTool.TEMPFILES list of auto-created tempfiles
              does not share state across processes, so things will not get
              cleaned up automatically as they do in a single-process
              pybedtools session.

            * this includes deleting any "chromsizes" or genome files --
              generally it will be best to require a genome filename in
              `func_kwargs` if you'll be using any BedTool methods that accept
              the `g` kwarg.

            * the function should be a module-level function (rather than a
              class method) because class methods can't be pickled across
              process boundaries

            * the function can have any signature and have any return value

        `_orig_pool` can be a previously-created multiprocessing.Pool instance;
        otherwise, a new Pool will be created with `processes`
        """
        if processes == 1:
            for _ in range(iterations):
                yield func(*func_args, **func_kwargs)
            raise StopIteration

        if _orig_pool:
            p = _orig_pool
        else:
            p = Pool(processes)
        iterations_each = [iterations / processes] * processes
        iterations_each[-1] += iterations % processes

        # FYI some useful info on apply_async:
        # http://stackoverflow.com/questions/8533318/
        #      python-multiprocessing-pool-when-to-use-apply-apply-async-or-map
        #
        # Here, we don't care about the order, and don't want the subprocesses
        # to block.
        results = [
            p.apply_async(func, func_args, func_kwargs) for _ in range(iterations)
        ]
        for r in results:
            yield r.get()
        raise StopIteration

    def random_jaccard(
        self,
        other: BedTool,
        genome_fn: Optional[str] = None,
        iterations: Optional[int] = None,
        processes: int = 1,
        _orig_pool: Optional[Pool] = None,
        shuffle_kwargs: Optional[dict[str, Any]] = None,
        jaccard_kwargs: Optional[dict[str, Any]] = None,
    ) -> list:
        """
        Computes the naive Jaccard statistic (intersection divided by union).

        .. note::

            If you don't need the randomization functionality of this method,
            you can use the simpler BedTool.jaccard method instead.

        See Favorov et al. (2012) PLoS Comput Biol 8(5): e1002529 for more
        info on the Jaccard statistic for intersections.

        If `iterations` is None, then do not perform random shufflings.

        If `iterations` is an integer, perform `iterations` random shufflings,
        each time computing the Jaccard statistic to build an empirical
        distribution.  `genome_fn` will also be needed; optional `processes`
        will split the iterations across multiple CPUs.

        Returns a tuple of the observed Jaccard statistic and a list of the
        randomized statistics (which will be an empty list if `iterations` was
        None).
        """
        if shuffle_kwargs is None:
            shuffle_kwargs = {}
        if jaccard_kwargs is None:
            jaccard_kwargs = {}
        if not genome_fn:
            raise ValueError("Need a genome filename in order to perform randomization")
        return list(
            self.parallel_apply(
                iterations=iterations if iterations else 1,
                func=pybedtools.stats.random_jaccard,
                func_args=(self, other),
                func_kwargs=dict(
                    genome_fn=genome_fn,
                    shuffle_kwargs=shuffle_kwargs,
                    jaccard_kwargs=jaccard_kwargs,
                ),
                processes=processes,
                _orig_pool=_orig_pool,
            )
        )

    def _randomintersection(
        self,
        other: BedTool,
        iterations: int,
        genome_fn: str,
        intersect_kwargs: Optional[dict[str, Any]] = None,
        _orig_pool:Optional[Pool] = None,
        shuffle_kwargs: Optional[dict[str, Any]] = None,
        processes: int = 1,
    ):
        """
        Re-implementation of BedTool.randomintersection using the new
        `random_op` method
        """
        if shuffle_kwargs is None:
            shuffle_kwargs = {}
        if intersect_kwargs is None:
            intersect_kwargs = dict(u=True)
        if not genome_fn:
            raise ValueError("Need a genome filename in order to perform randomization")
        return list(
            self.parallel_apply(
                iterations=iterations,
                func=pybedtools.stats.random_intersection,
                func_args=(self, other),
                func_kwargs=dict(
                    genome_fn=genome_fn,
                    shuffle_kwargs=shuffle_kwargs,
                    intersect_kwargs=intersect_kwargs,
                ),
                processes=processes,
                _orig_pool=_orig_pool,
            )
        )

    def randomintersection_bp(
        self,
        other: BedTool,
        iterations: int,
        genome_fn: str,
        intersect_kwargs: Optional[dict[str, Any]] = None,
        shuffle_kwargs: Optional[dict[str, Any]] = None,
        processes: int = 1,
        _orig_pool:Optional[Pool] = None,
    ) -> list[int]:
        """
        Like randomintersection, but return the bp overlap instead of the
        number of intersecting intervals.
        """
        if shuffle_kwargs is None:
            shuffle_kwargs = {}
        if intersect_kwargs is None:
            intersect_kwargs = {}
        if not genome_fn:
            raise ValueError("Need a genome filename in order to perform randomization")
        return list(
            self.parallel_apply(
                iterations=iterations,
                func=pybedtools.stats.random_intersection_bp,
                func_args=(self, other),
                func_kwargs=dict(
                    genome_fn=genome_fn,
                    shuffle_kwargs=shuffle_kwargs,
                    intersect_kwargs=intersect_kwargs,
                ),
                processes=processes,
                _orig_pool=_orig_pool,
            )
        )

    def randomintersection(
        self,
        other: BedTool,
        iterations: int,
        intersect_kwargs: Optional[dict[str, Any]] = None,
        shuffle_kwargs: Optional[dict[str, Any]] = None,
        debug: bool = False,
        report_iterations: bool = False,
        processes: Optional[int] = None,
        _orig_processes: Optional[int] = None,
    ) -> Iterator[int]:
        """
        Perform `iterations` shufflings, each time intersecting with `other`.

        Returns a generator of integers where each integer is the number of
        intersections of a shuffled file with *other*. This distribution can
        be used in downstream analysis for things like empirical p-values.

        *intersect_kwargs* and *shuffle_kwargs* are passed to self.intersect()
        and self.shuffle() respectively.  By default for intersect, u=True is
        specified -- but s=True might be a useful option for strand-specific
        work.

        Useful kwargs for *shuffle_kwargs* are chrom, excl, or incl.  If you
        use the "seed" kwarg, that seed will be used *each* time shuffleBed is
        called -- so all your randomization results will be identical for each
        iteration.  To get around this and to allow for tests, debug=True will
        set the seed to the iteration number.  You may also break up the
        intersections across multiple processes with *processes* > 1.

        Example usage:

            >>> chromsizes = {'chr1':(0, 1000)}
            >>> a = pybedtools.example_bedtool('a.bed')
            >>> a = a.set_chromsizes(chromsizes)
            >>> b = pybedtools.example_bedtool('b.bed')
            >>> results = a.randomintersection(b, 10, debug=True)
            >>> print(list(results))
            [1, 0, 1, 2, 4, 2, 2, 1, 2, 4]

        """
        if processes is not None:
            p = Pool(processes)
            iterations_each = [iterations // processes] * processes
            iterations_each[-1] += iterations % processes
            results = [
                p.apply_async(
                    _call_randomintersect,
                    (self, other, it),
                    dict(
                        intersect_kwargs=intersect_kwargs,
                        shuffle_kwargs=shuffle_kwargs,
                        debug=debug,
                        report_iterations=report_iterations,
                        _orig_processes=processes,
                    ),
                )
                for it in iterations_each
            ]
            for r in results:
                for value in r.get():
                    yield value
            raise StopIteration

        if shuffle_kwargs is None:
            shuffle_kwargs = {}
        if intersect_kwargs is None:
            intersect_kwargs = {"u": True}

        if "u" not in intersect_kwargs:
            intersect_kwargs["u"] = True

        resort = intersect_kwargs.get("sorted", False)

        for i in range(iterations):
            if debug:
                shuffle_kwargs["seed"] = i
            if report_iterations:
                if _orig_processes > 1:
                    msg = "\rapprox (total across %s processes): %s" % (
                        _orig_processes,
                        i * _orig_processes,
                    )
                else:
                    msg = "\r%s" % i
                sys.stderr.write(msg)
                sys.stderr.flush()

            # Re-sort if sorted=True in kwargs
            if resort:
                tmp0 = self.shuffle(**shuffle_kwargs)
                tmp = tmp0.sort()
            else:
                tmp = self.shuffle(**shuffle_kwargs)

            tmp2 = tmp.intersect(other, stream=True, **intersect_kwargs)

            yield len(tmp2)

            # Close the open stdouts from subprocess.Popen calls.  Note: doing
            # this in self.__del__ doesn't fix the open file limit bug; it
            # needs to be done here.
            # if resort:
            #     tmp0.fn.close()
            # tmp.fn.close()
            tmp2.fn.close()
            del tmp
            del tmp2

    @_log_to_history
    def cat(self, *others: Iterable[str|Path|BedTool], **kwargs) -> BedTool:
        """
        Concatenate interval files together.

        Concatenates two BedTool objects (or an object and a file) and does an
        optional post-merge of the features.

        *postmerge=True* by default; use *postmerge=False* if you want to keep
        features separate.

        *force_truncate=False* by default; *force_truncate=True* to truncate
        all files to chrom, start, stop.

        When *force_truncate=False* and *postmerge=False*, the output will
        contain the smallest number of fields observed across all inputs. This
        maintains compatibility with BEDTools programs, which assume constant
        number of fields in all lines of a file.

        Other kwargs are sent to :meth:`BedTool.merge` (and assuming that
        *postmerge=True*).

        Example usage:

        >>> a = pybedtools.example_bedtool('a.bed')
        >>> b = pybedtools.example_bedtool('b.bed')
        >>> print(a.cat(b)) #doctest: +NORMALIZE_WHITESPACE
        chr1	1	500
        chr1	800	950
        <BLANKLINE>
        >>> print(a.cat(*[b,b],
        ...   postmerge=False)) #doctest: +NORMALIZE_WHITESPACE
        chr1	1	100	feature1	0	+
        chr1	100	200	feature2	0	+
        chr1	150	500	feature3	0	-
        chr1	900	950	feature4	0	+
        chr1	155	200	feature5	0	-
        chr1	800	901	feature6	0	+
        chr1	155	200	feature5	0	-
        chr1	800	901	feature6	0	+
        <BLANKLINE>
        """
        same_type = None
        same_field_num = None
        field_nums = set()

        assert len(others) > 0, "You must specify at least one other bedfile!"
        other_beds = []
        for other in others:
            if isinstance(other, (str, Path)):
                other = BedTool(other)
            else:
                assert isinstance(
                    other, BedTool
                ), "Either filename or another BedTool instance required"
            other_beds.append(other)

        # postmerge and force_truncate don't get passed on to merge
        postmerge = kwargs.pop("postmerge", True)
        force_truncate = kwargs.pop("force_truncate", False)
        stream_merge = kwargs.get("stream", False)
        if stream_merge and postmerge:
            raise ValueError(
                "The post-merge step in the `cat()` method "
                "performs a sort, which uses stream=True.  Using "
                "stream=True for the merge as well will result in a "
                "deadlock!"
            )

        # if filetypes and field counts are the same, don't truncate
        if not force_truncate:
            try:
                filetypes = set(
                    [self.file_type] + [i.file_type for i in other_beds]
                ).difference(["empty"])
                field_nums = (
                    set([self.field_count()] + [i.field_count() for i in other_beds])
                    .difference([None])
                    .difference([0])
                )
                same_field_num = len(field_nums) == 1
                same_type = len(set(filetypes)) == 1
            except ValueError:
                raise ValueError(
                    "Can't check filetype or field count -- "
                    "is one of the files you're merging a 'streaming' "
                    "BedTool?  If so, use .saveas() to save to file first"
                )

        tmp = self._tmp()

        if not force_truncate and same_type and same_field_num:
            with open(tmp, "w") as TMP:
                for f in self:
                    TMP.write(str(f))
                for other in other_beds:
                    for f in other:
                        TMP.write(str(f))

        # Types match, so we can use the min number of fields observed across
        # all inputs
        elif not force_truncate and same_type:
            minfields = min(field_nums)
            with open(tmp, "w") as TMP:
                for f in self:
                    TMP.write("\t".join(f.fields[:minfields]) + "\n")
                for other in other_beds:
                    for f in other:
                        TMP.write("\t".join(f.fields[:minfields]) + "\n")

        # Otherwise, use the zero-based chrom/start/stop to create a BED3,
        # which will work when catting a GFF and a BED together.
        else:
            with open(tmp, "w") as TMP:
                for f in self:
                    TMP.write("%s\t%i\t%i\n" % (f.chrom, f.start, f.end))
                for other in other_beds:
                    for f in other:
                        TMP.write("%s\t%i\t%i\n" % (f.chrom, f.start, f.end))

        c = BedTool(tmp)
        if postmerge:
            d = c.sort(stream=True).merge(**kwargs)

            # Explicitly delete -- needed when using multiprocessing
            os.unlink(tmp)
            return d
        else:
            return c

    @_log_to_history
    def saveas(self, fn: Optional[str] = None, trackline: Optional[str] = None, compressed: Optional[bool] = None) -> BedTool:
        """
        Make a copy of the BedTool.

        Optionally adds `trackline` to the beginning of the file.

        Optionally compresses output using gzip.

        if the filename extension is .gz, or compressed=True,
        the output is compressed using gzip

        Returns a new BedTool for the newly saved file.

        A newline is automatically added to the trackline if it does not
        already have one.

        Example usage:

        >>> a = pybedtools.example_bedtool('a.bed')
        >>> b = a.saveas('other.bed')
        >>> b.fn
        'other.bed'
        >>> print(b == a)
        True

        >>> b = a.saveas('other.bed', trackline="name='test run' color=0,55,0")
        >>> open(b.fn).readline()
        "name='test run' color=0,55,0\\n"
        >>> if os.path.exists('other.bed'):
        ...     os.unlink('other.bed')
        """
        if fn is None:
            fn = self._tmp()

        # Default to compressed if extension is .gz
        if compressed is None:
            __, extension = os.path.splitext(fn)
            if extension == ".gz":
                compressed = True
            else:
                compressed = False

        in_compressed = isinstance(self.fn, str) and isGZIP(self.fn)

        fn = self._collapse(
            self,
            fn=fn,
            trackline=trackline,
            in_compressed=in_compressed,
            out_compressed=compressed,
        )
        return BedTool(fn)

    @_log_to_history
    def moveto(self, fn: Optional[str]=None) -> BedTool:
        """
        Move to a new filename (can be much quicker than BedTool.saveas())

        Move BED file to new filename, `fn`.

        Returns a new BedTool for the new file.

        Example usage:

        >>> # make a copy so we don't mess up the example file
        >>> a = pybedtools.example_bedtool('a.bed').saveas()
        >>> a_contents = str(a)
        >>> b = a.moveto('other.bed')
        >>> b.fn
        'other.bed'
        >>> b == a_contents
        True
        """
        if not isinstance(self.fn, str):
            fn = self._collapse(self, fn=fn)
        else:
            shutil.move(self.fn, fn)
        return BedTool(fn)

    @_log_to_history
    def random_subset(self, n: Optional[int] = None, f: Optional[float] = None, seed: Optional[float|int]=None) -> BedTool:
        """
        Return a BedTool containing a random subset.

        NOTE: using `n` will be slower and use more memory than using `f`.

        Parameters
        ----------

        n : int
            Number of features to return. Only one of `n` or `f` can be provided.

        f : float, 0 <= f <= 1
            Fraction of features to return. Cannot be provided with `n`.

        seed : float or int
            Set random.seed

        Example
        -------

        >>> seed = 0  # only for test, otherwise use None

        `n` will always give the same number of returned features, but will be
        slower since it is creating an index and then shuffling it.

        >>> a = pybedtools.example_bedtool('a.bed')
        >>> b = a.random_subset(n=2)
        >>> len(b)
        2

        Using a fraction `f` will be faster but depending on seed will result
        in slightly different total numbers.

        >>> a = pybedtools.example_bedtool('x.bam')
        >>> len(a)
        45593
        >>> b = a.random_subset(f=0.4, seed=seed)
        >>> len(b)
        18316

        Check that we have approximately the right fraction
        >>> print('{0:.2f}'.format(len(b) / len(a)))
        0.40

        """
        if ((n is None) and (f is None)) or ((n is not None) and (f is not None)):
            raise ValueError("Exactly one of `n` or `f` must be provided")

        tmpfn = self._tmp()
        if seed is not None:
            random.seed(seed)

        if n:
            idxs = list(range(len(self)))
            random.shuffle(idxs)
            idxs = idxs[:n]
            with open(tmpfn, "w") as tmp:
                for i, feature in enumerate(self):
                    if i in idxs:
                        tmp.write(str(feature))

        elif f:
            with open(tmpfn, "w") as tmp:
                for i in self:
                    if random.random() <= f:
                        tmp.write(str(i))

        return BedTool(tmpfn)

    def total_coverage(self) -> int:
        """
        Return the total number of bases covered by this interval file.

        Does a self.merge() first to remove potentially multiple-counting
        bases.

        Example usage:

        >>> a = pybedtools.example_bedtool('a.bed')

        This does a merge() first, so this is what the total coverage is
        counting:

        >>> print(a.merge()) #doctest: +NORMALIZE_WHITESPACE
        chr1	1	500
        chr1	900	950
        <BLANKLINE>

        >>> print(a.total_coverage())
        549
        """
        b = self.merge()
        total_bp = 0
        for feature in b.features():
            total_bp += len(feature)
        return total_bp

    @_log_to_history
    def with_attrs(self, **kwargs) -> BedTool:
        """
        Helper method for adding attributes in the middle of a pipeline.

        Given arbitrary keyword arguments, turns the keys and values into
        attributes.  Useful for labeling BedTools at creation time.

        Example usage:

        >>> # add a "label" attribute to each BedTool
        >>> a = pybedtools.example_bedtool('a.bed')\
                                   .with_attrs(label='transcription factor 1')
        >>> b = pybedtools.example_bedtool('b.bed')\
                                   .with_attrs(label='transcription factor 2')
        >>> for i in [a, b]:
        ...     print('{0} features for {1}'.format(i.count(), i.label))
        4 features for transcription factor 1
        2 features for transcription factor 2

        """
        for key, value in list(kwargs.items()):
            setattr(self, key, value)
        return self

    def as_intervalfile(self) -> IntervalFile:
        """
        Returns an IntervalFile of this BedTool for low-level interface.
        """
        if not isinstance(self.fn, str):
            fn = self._collapse(self.fn)
        else:
            fn = self.fn
        return IntervalFile(fn)

    def liftover(self, chainfile: str, unmapped: Optional[str] = None, liftover_args: str = "") -> BedTool:
        """
        Returns a new BedTool of the liftedOver features, saving the unmapped
        ones as `unmapped`.  If `unmapped` is None, then discards the unmapped
        features.

        `liftover_args` is a string of additional args that is passed,
        verbatim, to liftOver.

        Needs `liftOver` from UCSC to be on the path and a `chainfile`
        downloaded from UCSC.
        """
        result = BedTool._tmp()
        if unmapped is None:
            unmapped = BedTool._tmp()
        cmds = ["liftOver", liftover_args, self.fn, chainfile, result, unmapped]
        os.system(" ".join(cmds))
        return BedTool(result)

    def absolute_distance(self, other: BedTool, closest_kwargs: Optional[dict[str, Any]]=None, use_midpoints: bool=False) -> Iterator[int]:
        """
        Returns an iterator of the *absolute* distances between features in
        self and other.

        If `use_midpoints` is True, then only use the midpoints of features
        (which will return values where features are overlapping).  Otherwise,
        when features overlap the value will always be zero.

        `closest_kwargs` are passed to self.closest(); either `d` or
        'D` are required in order to get back distance values (`d=True` is
        default)
        """
        from .featurefuncs import midpoint

        if closest_kwargs is None:
            closest_kwargs = {"d": True}

        if "D" not in closest_kwargs:
            closest_kwargs.update(dict(d=True))

        if use_midpoints:
            mid_self = self.each(midpoint).saveas()
            mid_other = other.each(midpoint).saveas()
            c = mid_self.closest(mid_other, stream=True, **closest_kwargs)
        else:
            c = self.closest(other, stream=True, **closest_kwargs)
        for i in c:
            yield int(i[-1])

    def relative_distance(self, other: BedTool, genome:Optional[dict|str] =None, g: Optional[str]=None) -> Iterator[float]:
        """
        Returns an iterator of relative distances between features in self and
        other.

        First computes the midpoints of self and other, then returns distances
        of each feature in `other` relative to the distance between `self`
        features.

        Requires either `genome` (dictionary of chromsizes or assembly name) or
        `g` (filename of chromsizes file).
        """
        g_dict = {}
        if (genome is None) and (g is None):
            raise ValueError("Need either `genome` or `g` arg for relative distance")
        elif genome and g:
            raise ValueError("Please specify only one of `genome` or `g`")
        elif genome:
            g_dict = dict(genome=genome)
        elif g:
            g_dict = dict(g=g)

        from .featurefuncs import midpoint

        # This gets the space between features in self.
        c = self.each(midpoint).complement(**g_dict)

        hits = c.intersect(other, wao=True, stream=True) # TODO: should this be other or mid_other?
        for i in hits:
            yield float(i[-1]) / len(i)

    def colormap_normalize(self,
                            vmin: Optional[float|int]=None,
                            vmax: Optional[float|int]=None,
                            percentile: bool=False,
                            log: bool=False) -> mcolors.LogNorm|mcolors.Normalize:
        """
        Returns a normalization instance for use by featurefuncs.add_color().

        Parameters
        ----------
        vmin, vmax : float, int, or None
            `vmin` and `vmax` set the colormap bounds; if None then
            these will be determined from the scores in the BED file.

        log : bool
            If True, put the scores on a log scale; of course be careful
            if you have negative scores

        percentile : bool
            If True, interpret vmin and vmax as a percentile in the range
            [0,100] rather than absolute values.
        """
        field_count = self.field_count()
        if (self.file_type != "bed") or (field_count < 5):
            raise ValueError("colorizing only works for BED files with score " "fields")
        try:
            import matplotlib.colors as mcolors
        except ImportError:
            raise ImportError("matplotlib.colors must be installed to use colormap_normalize")

        try:
            import numpy as np
        except ImportError:
            raise ImportError("numpy must be installed to use colormap_normalize")

        if log:
            norm = mcolors.LogNorm()
        else:
            norm = mcolors.Normalize()

        scores = np.array([i.score for i in self], dtype=float)
        scores = scores[np.isfinite(scores)]
        norm.autoscale(scores)

        if vmin is not None:
            if percentile:
                vmin = float(np.percentile(scores, vmin))
            norm.vmin = vmin
        if vmax is not None:
            if percentile:
                vmax = float(np.percentile(scores, vmax))
            norm.vmax = vmax

        return norm

    def at(self, inds: list[int]) -> BedTool:
        """
        Returns a new BedTool with only intervals at lines `inds`

        Parameters
        ----------
        inds : List[int]
            List of line numbers

        Returns
        -------
        BedTool
            New BedTool with only intervals at `inds`
        """
        length = len(inds)

        def _gen():
            k = 0
            for i, feature in enumerate(self):
                if i == inds[k]:
                    yield feature
                    k += 1
                    if k == length:
                        break

        return BedTool(_gen()).saveas()

    def to_dataframe(self, disable_auto_names: bool = False, *args, **kwargs) -> pd.DataFrame:
        """
        Create a pandas.DataFrame, passing args and kwargs to pandas.read_csv
        The separator kwarg `sep` is given a tab `\\t` as value by default.

        Parameters
        ----------
        disable_auto_names : bool
            By default, the created dataframe fills in column names
            automatically according to the detected filetype (e.g., "chrom",
            "start", "end" for a BED3 file). Set this argument to True to
            disable this behavior.
        """
        # Complain if BAM or if not a file
        if self._isbam:
            raise ValueError("BAM not supported for converting to DataFrame")
        if not isinstance(self.fn, str):
            raise ValueError("use .saveas() to make sure self.fn is a file")

        try:
            import pandas
        except ImportError:
            raise ImportError("pandas must be installed to convert to pandas.DataFrame")
        # Otherwise we're good:
        names = kwargs.get("names", None)
        if names is None and not disable_auto_names:
            try:
                _names = settings._column_names[self.file_type][: self.field_count()]
                if len(_names) < self.field_count():
                    warn(
                        "Default names for filetype %s are:\n%s\nbut file has "
                        "%s fields; you can supply custom names with the "
                        "`names` kwarg" % (self.file_type, _names, self.field_count())
                    )
                    _names = None
            except KeyError:
                _names = None
            kwargs["names"] = _names

        if os.path.isfile(self.fn) and os.path.getsize(self.fn) > 0:
            return pandas.read_csv(self.fn, *args, sep="\t", **kwargs) # type: ignore
        else:
            return pandas.DataFrame()

    def tail(self, lines:int = 10, as_string: bool = False) -> Optional[str]:
        """
        Like `head`, but prints last 10 lines of the file by default.

        To avoid consuming iterables, this only works with file-based, non-BAM
        BedTool objects.

        Use `as_string=True` to return a string.
        """
        if self._isbam:
            raise ValueError("tail() not yet implemented for BAM files")
        if not isinstance(self.fn, str):
            raise ValueError(
                "tail() not implemented for non-file-based "
                "BedTool objects.  Please use saveas() first."
            )
        bufsize = 8192
        offset = bufsize
        f = open(self.fn, "rb")

        # whence=2 arg means relative to end (i.e., go to the end)
        f.seek(0, 2)
        file_size = f.tell()
        data = []
        while True:
            if file_size < bufsize:
                offset = file_size
            f.seek(-offset, 2)
            chunk = f.read(offset)
            data.extend(chunk.splitlines(True))
            if len(data) >= lines or offset == file_size:
                break
            offset += bufsize

        result = "".join([i.decode() for i in data[-lines:]])
        if as_string:
            return result
        else:
            print(result)


class BAM(object):
    def __init__(self, stream):
        """
        Wraps pysam.Samfile so that it yields pybedtools.Interval objects when
        iterated over.

        The pysam.Samfile can be accessed via the .pysam_bamfile attribute.
        """
        self.stream = stream
        if not isinstance(self.stream, str):
            raise ValueError("Only files are supported, not streams")
        self.pysam_bamfile = pysam.Samfile(self.stream)

    def _aligned_segment_to_interval(self, r):

        if r.rname >= 0:
            rname = self.pysam_bamfile.get_reference_name(r.rname)
        else:
            rname = "*"

        if r.rnext >= 0:
            if r.rnext == r.rname:
                rnext = "="
            else:
                rnext = self.pysam_bamfile.get_reference_name(r.rnext)
        else:
            rnext = "*"

        # SAM spec says if unavailable should be set to 0. Pysam sets to -1.

        if r.pnext <= 0:
            pnext = "0"
        else:
            # +1 here because cbedtools.pyx expects SAM -- which is 1-based --
            # but pysam uses 0-based.
            pnext = str(r.pnext + 1)

        if r.cigarstring:
            cigarstring = r.cigarstring
        else:
            cigarstring = "*"

        # Rudimentary support.
        # TODO: remove when refactoring to new BAM iterating
        tags = []
        for k, v in r.tags:
            if isinstance(v, int):
                t = "i"
            elif isinstance(v, float):
                t = "f"
            else:
                t = "Z"
            tags.append("{0}:{1}:{2}".format(k, t, v))

        tags = "\t".join(tags)

        if r.seq:
            seq = r.seq
        else:
            seq = "*"

        if r.qual:
            qual = r.qual
        else:
            qual = "*"

        fields = [
            r.qname,
            str(r.flag),
            rname,
            # +1 here because cbedtools.pyx expects SAM -- which is 1-based --
            # but pysam uses 0-based.
            str(r.pos + 1),
            str(r.mapq),
            cigarstring,
            rnext,
            pnext,
            str(r.tlen),
            seq,
            qual,
        ]
        if tags:
            fields.append(tags)

        if None in fields:
            raise ValueError("Found 'None' in fields: %s" % fields)
        return create_interval_from_list(fields)

    def __iter__(self):
        return self

    # TODO: this is PAINFUL but it ensures that existing tests work.  Once all
    # tests work, the new behavior will be to yield pysam AlignedSegment
    # objects directly.
    def __next__(self):
        return self._aligned_segment_to_interval(next(self.pysam_bamfile))

    def next(self):
        return self.__next__()


class History(list):
    def __init__(self):
        """
        Represents one or many HistorySteps.  Mostly used for nicely formatting
        a series of HistorySteps.
        """
        list.__init__(self)


class HistoryStep(object):
    def __init__(self, method, args, kwargs, bedtool_instance, parent_tag, result_tag):
        """
        Class to represent one step in the history.

        Mostly used for its __repr__ method, to try and exactly replicate code
        that can be pasted to re-do history steps
        """
        try:
            self.method = method._name
        except AttributeError:
                self.method = method.__name__

        self.args = args
        self.kwargs = kwargs
        self.fn = bedtool_instance.fn
        self.parent_tag = parent_tag
        self.result_tag = result_tag

    def _clean_arg(self, arg: str|BedTool) -> str:
        """
        Wrap strings in quotes and convert bedtool instances to filenames.
        """
        if isinstance(arg, BedTool):
            arg = arg.fn
        if isinstance(arg, str):
            arg = '"%s"' % arg
        return arg

    def __repr__(self):
        # Still not sure whether to use pybedtools.bedtool() or bedtool()
        s = ""
        s += "<HistoryStep> "
        if os.path.exists(self.fn):
            s += 'BedTool("%(fn)s").%(method)s(%%s%%s)' % self.__dict__
        else:
            s += 'BedTool("MISSING FILE: %(fn)s")' % self.__dict__
            s += ".%(method)s(%%s%%s)" % self.__dict__

        # Format args and kwargs
        args_string = ",".join(map(self._clean_arg, self.args))
        kwargs_string = ",".join(
            ["%s=%s" % (i[0], self._clean_arg(i[1])) for i in list(self.kwargs.items())]
        )
        # stick a comma on the end if there's something here
        if len(args_string) > 0:
            args_string += ", "

        s = s % (args_string, kwargs_string)
        s += ", parent tag: %s" % self.parent_tag
        s += ", result tag: %s" % self.result_tag
        return s


def example_bedtool(fn):
    """
    Return a bedtool using a bed file from the pybedtools examples directory.
    Use :func:`list_example_files` to see a list of files that are included.
    """
    fn = os.path.join(filenames.data_dir(), fn)
    if not os.path.exists(fn):
        msg = "%s does not exist" % fn
        raise FileNotFoundError(msg)
    return BedTool(fn)


if __name__ == "__main__":
    import doctest

    doctest.testmod(optionflags=doctest.ELLIPSIS | doctest.NORMALIZE_WHITESPACE)