pybedtools

To simplify usage, pybedtools automatically handles common BEDTools input flags (-i, -a, and -b):

• Single Input (-i): For methods like .merge(), pybedtools assumes the input file is the one the current BedTool object points to. You do not need to pass the i argument explicitly.
• Two Inputs (-a and -b): For methods like .intersect(), the BedTool object itself is treated as -a, and the first non-keyword argument provided is treated as -b.
• Argument Types: The first non-keyword argument for -b can be either a filename (string) or another BedTool object.

Methods that use these defaults will include a note in their docstring: "For convenience, the file this BedTool object points to is passed as '-i'" or "This method accepts either a BedTool or a file name as the first unnamed argument".

# Default behavior: 'a' is used as input (-i)
result1 = a.merge(d=100, s=True)

# Explicit behavior (identical to above)
result2 = a.merge(i=a.fn, d=100, s=True)

# Two-input behavior: 'a' is -a, 'b' is -b
result3 = a.intersect(b)

# Using a filename or BedTool object for the second input
result5 = a.intersect(b=b.fn)
result6 = a.intersect(b=b)

To save disk space and memory, you can use BedTool objects as streaming generators. Instead of loading all data into memory or writing to disk, data is created one line at a time.

To enable streaming, pass stream=True as a keyword argument to operations like .intersect().

Warning: Streaming iterators are consumed when you read their features. For example, calling len() on a streaming BedTool will consume the iterator, making subsequent calls return 0.

Note on Rendering: Because many BEDTools programs require the second input (-b) to be a file, pybedtools may automatically "render" a streaming BedTool to a temporary file if it is passed as a secondary argument in an operation.

a = pybedtools.example_bedtool('a.bed')
b = pybedtools.example_bedtool('b.bed')
# Create a streaming BedTool
c = a.intersect(b, stream=True)

# len(c) consumes the iterator
print(len(c)) # returns length
print(len(c)) # returns 0

Some BEDTools commands (like groupby) can produce output files that do not follow standard interval formats (e.g., BED, GFF, GTF, BAM, SAM, VCF).

Limitations:

• You cannot use .saveas() on these results because pybedtools cannot wrap an unsupported file format into a new BedTool object.
• You cannot iterate over these results directly (e.g., [i for i in d]) because the lines cannot be converted into Interval objects.

Workarounds:

1. To save the file: Use the output='filename' keyword argument within the command that generates the non-interval data (e.g., c.groupby(..., output='counts.txt')). Note: This behavior depends on the underlying bedtools version.
2. To read the data: Use standard Python file handling to iterate over the lines of the file using the .fn attribute of the BedTool object to get the filename.

# To iterate over lines of a non-interval file:
for line in open(d.fn):
    featuretype, count = line.strip().split()

Different file formats use different coordinate systems:

• BED files: 0-based, stop position is not included.
• GFF, GTF, and VCF files: 1-based, stop position is included.

pybedtools normalizes these so you can treat all Interval objects identically:

1. Interval.start is always 0-based: Even if the source is a 1-based GFF file, feature.start will return the 0-based coordinate.
2. len(interval) is always correct: It always returns stop - start.
3. Interval.fields preserves original format: The fields list (and string indexing) always contains the original strings from the file. For a GFF file, feature[3] will return the 1-based string coordinate, while feature.start returns the 0-based integer.

Summary Table:

# GFF (1-based) vs BED (0-based) comparison
# GFF line: chr1, start 51, stop 300
gff = pybedtools.create_interval_from_list(["chr1", "fake", "mRNA", "51", "300", ".", "+", ".", "ID=mRNA1;"])

# BED line representing same region: chr1, start 50, stop 300
bed = pybedtools.create_interval_from_list(["chr1", "50", "300", "mRNA1", ".", "+"])

# Both have the same 0-based start attribute
assert gff.start == bed.start  # 50

# The GFF string field remains 1-based
assert gff[3] == '51'

# Length is consistent
assert len(gff) == len(bed)  # 250

When calling methods on a BedTool object (such as .intersect(), .subtract(), or .closest()), pybedtools automatically assumes that the file associated with the current BedTool instance (stored in the .fn attribute) is the first input file (the -a argument in bedtools command line terms).

For methods that require a second file, the first unnamed argument passed to the method is automatically treated as the second input file (the -b argument). If you pass another BedTool object as that argument, pybedtools will automatically use its .fn attribute.

>>> import pybedtools
>>> exons = pybedtools.example_bedtool('a.bed')
>>> snps = pybedtools.example_bedtool('b.bed')
>>> # The following are all equivalent:
>>> x1 = exons.intersect(snps)
>>> x2 = exons.intersect(a=exons.fn, b=snps.fn)
>>> x3 = exons.intersect(b=snps.fn)
>>> x4 = exons.intersect(snps, a=exons.fn)
>>> x1 == x2 == x3 == x4
True

When transforming genomic data, choose between each, filter, and map based on your desired output:

1. BedTool.each(func): Use this when you want to modify the data within the intervals (e.g., updating the score field) and want the result to remain a BedTool object for further processing (chaining).
2. BedTool.filter(func): Use this when you want to subset the data. The function must return True or False.
3. map(func, bedtool_obj): Use the standard Python map if you want to apply a function to each interval but do not need the result to be a BedTool object (e.g., extracting a list of lengths).

# Returns a BedTool object (chainable)
a.intersect(b).each(normalize_count)

# Returns a BedTool object (subsetting)
a.filter(lambda x: float(x[4]) < 1e-5)

# Returns an iterator of the function results (not a BedTool)
feature_lengths = map(len, a)