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

リビジョン 3, 2.1 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 (start, stop) and a maf. Produces a new maf
5containing the blocks from the original that overlapped the intervals.
6
7NOTE: See maf_extract_ranges_indexed.py which works better / faster for many
8      use cases.
9
10NOTE: chromosome/src information in the MAF is ignored by this variant.
11
12usage: %prog interval_file refindex [options] < maf_file
13   -m, --mincols=10: Minimum length (columns) required for alignment to be output
14"""
15
16import psyco_full
17
18from bx.cookbook import doc_optparse
19
20import bx.align.maf
21from bx import intervals
22import sys
23
24
25def __main__():
26
27    # Parse Command Line
28
29    options, args = doc_optparse.parse( __doc__ )
30
31    try:
32        range_filename = args[ 0 ]
33        refindex = int( args[ 1 ] )
34        if options.mincols: mincols = int( options.mincols )
35        else: mincols = 10
36    except:
37        doc_optparse.exit()
38
39    # Load Intervals
40
41    intersecter = intervals.Intersecter()
42    for line in file( range_filename ):
43        fields = line.split()
44        intersecter.add_interval( intervals.Interval( int( fields[0] ), int( fields[1] ) ) )
45
46    # Start MAF on stdout
47
48    out = bx.align.maf.Writer( sys.stdout )
49
50    # Iterate over input MAF
51
52    for maf in bx.align.maf.Reader( sys.stdin ):
53        ref = maf.components[ refindex ]
54        # Find overlap with reference component
55        intersections = intersecter.find( ref.get_forward_strand_start(), ref.get_forward_strand_end() )
56        # Keep output maf ordered
57        intersections.sort()
58        # Write each intersecting block
59        for interval in intersections:
60            start = max( interval.start, ref.get_forward_strand_start() )
61            end = min( interval.end, ref.get_forward_strand_end() )
62            sliced = maf.slice_by_component( refindex, start, end )
63            good = True
64            for c in sliced.components:
65                if c.size < 1:
66                    good = False
67            if good and sliced.text_size > mincols: out.write( sliced )
68         
69    # Close output MAF
70
71    out.close()
72
73if __name__ == "__main__": __main__()
Note: リポジトリブラウザについてのヘルプは TracBrowser を参照してください。