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

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

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

行番号 
1#!/usr/bin/python2.6
2
3"""
4Reads a list of intervals and an axt. Produces a new axt containing the
5portions of the original that overlapped the intervals
6
7usage: %prog interval_file refindex [options] < axt_file
8   -m, --mincols=10: Minimum length (columns) required for alignment to be output
9"""
10
11import psyco_full
12
13from bx.cookbook import doc_optparse
14
15import bx.align.axt
16from bx import intervals
17import sys
18
19
20def __main__():
21
22        # Parse Command Line
23
24        options, args = doc_optparse.parse( __doc__ )
25
26        try:
27                range_filename = args[ 0 ]
28                refindex = int( args[ 1 ] )
29                if options.mincols: mincols = int( options.mincols )
30                else: mincols = 10
31        except:
32                doc_optparse.exit()
33
34        # Load Intervals
35
36        intersecter = intervals.Intersecter()
37        for line in file( range_filename ):
38                fields = line.split()
39                intersecter.add_interval( intervals.Interval( int( fields[0] ), int( fields[1] ) ) )
40
41        # Start axt on stdout
42
43        out = bx.align.axt.Writer( sys.stdout )
44
45        # Iterate over input axt
46
47        for axt in bx.align.axt.Reader( sys.stdin ):
48                ref_component = axt.components[ refindex ]
49                # Find overlap with reference component
50                intersections = intersecter.find( ref_component.start, ref_component.end )
51                # Keep output axt ordered
52                intersections.sort()
53                # Write each intersecting block
54                for interval in intersections:
55                        start = max( interval.start, ref_component.start )
56                        end = min( interval.end, ref_component.end )
57                        sliced = axt.slice_by_component( refindex, start, end )
58                        good = True
59                        for c in sliced.components:
60                                if c.size < 1:
61                                        good = False
62                        if good and sliced.text_size > mincols: out.write( sliced )
63                 
64        # Close output axt
65
66        out.close()
67
68if __name__ == "__main__": __main__()
Note: リポジトリブラウザについてのヘルプは TracBrowser を参照してください。