.. include:: includeme.rst .. _`BedTools as iterators`: Using BedTool objects as iterators/generators ============================================= Typically, :mod:`BedTool` objects are used somewhat like handles to individual files on disk that contain BED lines. To save disk space, :mod:`BedTool` objects also have the ability to "stream", much like piping in Unix. That is, the data are created only one line at a time in memory, instead of either creating a list of all data in memory or writing all data to disk. .. note:: You'll need to be careful when using :mod:`BedTool` objects as generators, since any operation that reads all the features of a :mod:`BedTool` will consume the iterable. To get a streaming BedTool, use the `stream=True` kwarg. This :class:`BedTool` will act a little differently from a standard, file-based :class:`BedTool`. .. doctest:: >>> a = pybedtools.example_bedtool('a.bed') >>> b = pybedtools.example_bedtool('b.bed') >>> c = a.intersect(b, stream=True) >>> # checking the length consumes the iterator >>> len(c) 3 >>> # nothing left, so checking length again returns 0 >>> len(c) 0 In some cases, a stream may be "rendered" to a temp file. This is because BEDTools programs can only accept one input file as `stdin`. This is typically the first input (`-i` or `-a`), while the other input (`-b`) must be a file. Consider this example, where the second intersection needs to convert the streaming BedTool to a file before sending to `intersectBed`: .. doctest:: >>> a = pybedtools.example_bedtool('a.bed') >>> b = pybedtools.example_bedtool('b.bed') >>> # first we set up a streaming BedTool: >>> c = a.intersect(b, stream=True) >>> # But supplying a streaming BedTool as the first unnamed argument >>> # means it is being passed as -b to intersectBed, and so must be a file. >>> # In this case, `c` is rendered to a tempfile before being passed. >>> d = a.intersect(c, stream=True) .. warning:: Chaining two streaming BedTool objects together? You'll need to be careful, because sometimes this will result in deadlocks. You'll see zero CPU usage as pybedtools tries to write to the stdin of a downstream BedTool object. For example, for two BedTool objects pointing to files that have >5000 features, the following usually blocks:: len(a.intersect(b, stream=True).intersect(a, stream=True)) The solution is to save to a tempfile first, or use non-streaming BedTools. All of the following will work fine:: >>> # only use file-based >>> len(a.intersect(b).intersect(a)) >>> # using the second streaming BedTool is fine >>> len(a.intersect(b).intersect(a, stream=True)) >>> # if you have a streaming BedTool, "render" it to a tempfile with >>> # saveas() >>> len(a.intersect(b, stream=True).saveas().intersect(a)) There's a nice explanation of blocking along with figures at http://www.pixelbeat.org/programming/stdio_buffering/. Most solutions to this blocking problem on stackoverflow suggest using threads, but in my test cases this tends to make interactive IPython sessions act strangely. Another option is to try pexpect, but I have been unable to get this to work and it requires an additional dependency. Contributions to help solve this would be most appreciated! Example references: * http://stackoverflow.com/questions/1595492/blocks-send-input-to-python-subprocess-pipeline * http://stackoverflow.com/questions/375427/non-blocking-read-on-a-subprocess-pipe-in-python * http://www.python.org/dev/peps/pep-3145/ * http://stackoverflow.com/questions/3076542/how-can-i-read-all-availably-data-from-subprocess-popen-stdout-non-blocking/3078292#3078292 * http://stackoverflow.com/questions/3140189/subprocess-popen-stdout-reading-stdout-in-real-time-again (and references therein) Creating a :class:`BedTool` from an iterable -------------------------------------------- You can create a :class:`BedTool` on the fly from a generator or iterator -- in fact, this is what the :meth:`BedTool.filter` method does for you: .. doctest:: :options: +NORMALIZE_WHITESPACE >>> a = pybedtools.example_bedtool('a.bed') >>> print(a) chr1 1 100 feature1 0 + chr1 100 200 feature2 0 + chr1 150 500 feature3 0 - chr1 900 950 feature4 0 + >>> b = pybedtools.BedTool(f for f in a if f.start > 200) >>> # this is the same as using filter: >>> c = a.filter(lambda x: x.start > 200) We need to "render" these BedTools to string before we can check equality -- consuming them both -- since they are both iterables for which `==` is not defined: .. doctest:: :options: +ELLIPSIS >>> b == c Traceback (most recent call last): ... NotImplementedError: Testing equality only supported for BedTools that point to a file >>> str(b) == str(c) True Indexing a :class:`BedTool` --------------------------- In some cases it may be useful to index into a :class:`BedTool` object. We can use standard list slice syntax, and get an iterable of :class:`Interval` objects as a result. This iterable can in turn be used to create a new :class:`BedTool` instance: .. doctest:: :options: +NORMALIZE_WHITESPACE +ELLIPSIS >>> a = pybedtools.example_bedtool('a.bed') >>> a[2:4] >>> for i in a[2:4]: ... print(i) chr1 150 500 feature3 0 - chr1 900 950 feature4 0 + >>> b = pybedtools.example_bedtool('b.bed') >>> print(pybedtools.BedTool(a[:3]).intersect(b)) chr1 155 200 feature2 0 + chr1 155 200 feature3 0 -