1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
|
#!/usr/bin/python3
"""
Find regions of first bed file that overlap regions in a second bed file. This
program performs a base-by-base intersection, so only runs of bases that are
covered in both of the inputs will be output.
usage: %prog bed_file_1 bed_file_2
"""
from bx.bitset_builders import binned_bitsets_from_file
from bx.cookbook import doc_optparse
options, args = doc_optparse.parse(__doc__)
try:
in_fname, in2_fname = args
except ValueError:
doc_optparse.exit()
bits1 = binned_bitsets_from_file(open(in_fname))
bits2 = binned_bitsets_from_file(open(in2_fname))
bitsets = {}
for key in bits1:
if key in bits2:
bits1[key].iand(bits2[key])
bitsets[key] = bits1[key]
for chrom in bitsets:
bits = bitsets[chrom]
end = 0
while True:
start = bits.next_set(end)
if start == bits.size:
break
end = bits.next_clear(start)
print("%s\t%d\t%d" % (chrom, start, end))
|
