intersecting multiple files produces different default output that bedtools

[keoughkath]

Hello,

thanks for the awesome and super useful package! I've noticed that intersecting multiple files produces different output than the command line bedtools, which seems like a potential problem.

From the bedtools tutorial, here's the command line code:

bedtools intersect -a exons.bed -b cpg.bed gwas.bed hesc.chromHmm.bed -sorted | head

and here's the pybedtools version:

exons.intersect(cpg+gwas+hesc, sorted=True).head()

which produces a different output. However, adding loj=True,

exons.intersect(cpg+gwas+hesc, sorted=True, loj=True).head()

renders the output the same as the command-line tool, and what's expected from the tutorial. Am I doing something incorrectly?

Thanks!
Kathleen

[wreckenrodd]

Good question!

[vroide]

Hello keoughkath,

It seems, that loj=True gives you the correct output by accident, because you are only looking at the head of the file.

bedtools intersect -a exons.bed -b cpg.bed gwas.bed hesc.chromHmm.bed -sorted |
 grep chr | wc -l
541210

from pybedtools import BedTool
exons, cpg, gwas, hesc = (BedTool('exons.bed'),
                          BedTool('cpg.bed'),
                          BedTool('gwas.bed'),
                          BedTool('hesc.chromHmm.bed'))

loj = exons.intersect(cpg + gwas + hesc, sorted=True, loj=True)
len(loj)
# 436

You get the correct result like this, by passing a list (or other iterable):

correct = exons.intersect([cpg, gwas, hesc])
len(correct)
# 541210

This is the same behaviour as on the command line, where you also pass a list
of files to -b.

As far as i understand, the '+' operator is short for BedTool.intersect() and
thus you may get additional, unwanted intersections.

Maybe another Issue or just me not having read the pybedtools docs properly:

What puzzles me though, is that this does not give the same result:

plus = exons + [cpg, gwas, hesc]
len(plus)
# 436556

Is this expected?

[daler]

Sorry for the late reply here . . .

The overloaded + operator does an intersect with the -u argument, so these should give the same results:

x1 = exons + cpg
x2 = exons.intersect(cpg, u=True)

so I think that last example would be the equivalent of:

x1 = exons + [cpg, gwas, hesc]
x2 = exons.intersect([cpg, gwas, hesc], u=True)

[keoughkath]

Ah, that makes sense. Thanks!