| 1 | #!/usr/bin/python2.6 |
|---|
| 2 | |
|---|
| 3 | """ |
|---|
| 4 | Match up intersecting intervals from two files. This performs a "full join", |
|---|
| 5 | any pair of intervals with any basewise overlap will be printed side-by-side. |
|---|
| 6 | |
|---|
| 7 | usage: %prog bed1 bed2 |
|---|
| 8 | """ |
|---|
| 9 | |
|---|
| 10 | from __future__ import division |
|---|
| 11 | |
|---|
| 12 | import psyco_full |
|---|
| 13 | |
|---|
| 14 | import string |
|---|
| 15 | import sys |
|---|
| 16 | |
|---|
| 17 | import bx.intervals.io |
|---|
| 18 | import bx.intervals.intersection |
|---|
| 19 | |
|---|
| 20 | def main(): |
|---|
| 21 | |
|---|
| 22 | intersecters = {} |
|---|
| 23 | |
|---|
| 24 | # Read second set into intersecter |
|---|
| 25 | for interval in bx.intervals.io.GenomicIntervalReader( open( sys.argv[2] ) ): |
|---|
| 26 | if not intersecters.has_key( interval.chrom ): |
|---|
| 27 | intersecters[ interval.chrom ] = bx.intervals.Intersecter() |
|---|
| 28 | intersecters[ interval.chrom ].add_interval( interval ) |
|---|
| 29 | |
|---|
| 30 | # Join with first set |
|---|
| 31 | for interval in bx.intervals.io.GenomicIntervalReader( open( sys.argv[1] ) ): |
|---|
| 32 | if intersecters.has_key( interval.chrom ): |
|---|
| 33 | intersection = intersecters[ interval.chrom ].find( interval.start, interval.end ) |
|---|
| 34 | for interval2 in intersection: |
|---|
| 35 | print "\t".join( [ str( interval ), str( interval2 ) ] ) |
|---|
| 36 | |
|---|
| 37 | |
|---|
| 38 | if __name__ == "__main__": |
|---|
| 39 | main() |
|---|