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

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

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

行番号 
1#!/usr/bin/python2.6
2
3"""
4Complement the regions of a bed file. Requires a file that maps source names
5to sizes. This should be in the simple LEN file format (each line contains
6a source name followed by a size, separated by whitespace).
7
8usage: %prog bed_file chrom_length_file
9"""
10
11import sys
12
13from bx.bitset import *
14from bx.bitset_builders import *
15
16from bx.cookbook import doc_optparse
17
18def read_len( f ):
19    """Read a 'LEN' file and return a mapping from chromosome to length"""
20    mapping = dict()
21    for line in f:
22        fields = line.split()
23        mapping[ fields[0] ] = int( fields[1] )
24    return mapping
25
26options, args = doc_optparse.parse( __doc__ )
27try:
28    in_fname, len_fname = args
29except:
30    doc_optparse.exit()
31
32bitsets = binned_bitsets_from_file( open( in_fname ) )
33
34lens = read_len( open( len_fname ) )
35
36for chrom in lens:
37    if chrom in bitsets:
38        bits = bitsets[chrom]
39        bits.invert()
40        len = lens[chrom]
41        end = 0
42        while 1:
43            start = bits.next_set( end )
44            if start == bits.size: break
45            end = bits.next_clear( start )
46            if end > len: end = len
47            print "%s\t%d\t%d" % ( chrom, start, end )
48            if end == len: break
49    else:
50        print "%s\t%d\t%d" % ( chrom, 0, lens[chrom] )
Note: リポジトリブラウザについてのヘルプは TracBrowser を参照してください。