| 1 | #!/usr/bin/python2.6 |
|---|
| 2 | |
|---|
| 3 | """ |
|---|
| 4 | Print number of bases covered by all intervals in a bed file (bases covered by |
|---|
| 5 | more than one interval are counted only once). Multiple bed files can be |
|---|
| 6 | provided on the command line or to stdin. |
|---|
| 7 | |
|---|
| 8 | usage: %prog bed files ... |
|---|
| 9 | """ |
|---|
| 10 | |
|---|
| 11 | import psyco_full |
|---|
| 12 | import sys |
|---|
| 13 | from bx.bitset import BinnedBitSet |
|---|
| 14 | from bx.bitset_builders import * |
|---|
| 15 | from itertools import * |
|---|
| 16 | |
|---|
| 17 | bed_filenames = sys.argv[1:] |
|---|
| 18 | if bed_filenames: |
|---|
| 19 | input = chain( * imap( open, bed_filenames ) ) |
|---|
| 20 | else: |
|---|
| 21 | input = sys.stdin |
|---|
| 22 | |
|---|
| 23 | bitsets = binned_bitsets_from_file( input ) |
|---|
| 24 | |
|---|
| 25 | total = 0 |
|---|
| 26 | for chrom in bitsets: |
|---|
| 27 | total += bitsets[chrom].count_range( 0, bitsets[chrom].size ) |
|---|
| 28 | |
|---|
| 29 | print total |
|---|