root/galaxy-central/eggs/bx_python-0.5.0_dev_f74aec067563-py2.6-macosx-10.6-universal-ucs2.egg/EGG-INFO/scripts/interval_count_intersections.py

リビジョン 3, 1.2 KB (コミッタ: kohda, 14 年 前)

Install Unix tools  http://hannonlab.cshl.edu/galaxy_unix_tools/galaxy.html

行番号 
1#!/usr/bin/python2.6
2
3"""
4Read two lists of intervals (with chromosomes) and count the number of entries
5in the second set that intersect any entry in the first set.
6
7TODO: This could use bitsets rather than the intervals package, would it be
8      faster?
9
10usage: %prog bed1 bed2 > out
11"""
12
13from __future__ import division
14
15import psyco_full
16
17from bx import intervals
18from bx import misc
19import string
20import sys
21
22def main():
23
24    intersecters = {}
25
26    # Read ranges
27
28    for chr, start, end in read_intervals( misc.open_compressed( sys.argv[1] ) ):
29        if not intersecters.has_key( chr ): intersecters[ chr ] = intervals.Intersecter()
30        intersecters[ chr ].add_interval( intervals.Interval( start, end ) )
31
32    # Count intersection
33
34    total = 0
35
36    for chr, start, end in read_intervals( misc.open_compressed( sys.argv[2] ) ):
37        if intersecters.has_key( chr ):
38            intersection = intersecters[ chr ].find( start, end )
39            if intersection:
40                #print chr, intersection
41                total += 1
42
43    print total
44
45def read_intervals( input ):
46    for line in input:
47        fields = line.split()
48        yield fields[0], int( fields[1] ), int( fields[2] )
49
50main()
Note: リポジトリブラウザについてのヘルプは TracBrowser を参照してください。