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

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