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

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

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

行番号 
1#!/usr/bin/python2.6
2
3import sys
4
5import bx.align.maf
6import bx.bitset
7from bx.bitset_builders import *
8from itertools import *
9from optparse import OptionParser
10from rpy import *
11
12def main():
13
14        # Parse the command line
15        parser = OptionParser(usage = "usage: %prog [options] maf_file snp_file neutral_file window_size step_size")
16       
17        parser.add_option("-o", "--outfile", help = "Specify file for output")
18       
19        parser.add_option("-s", "--species", type = "string", default = "panTro2")
20       
21        parser.add_option("-b", "--build", type = "string", default = "hg18")
22       
23        (options, args) = parser.parse_args()
24
25        if len(args) != 5:
26                parser.error("Incorrect number of arguments")
27        else:
28                maf_filename = args[0]
29                snp_filename = args[1]
30                neutral_filename = args[2]
31                window_size = int(args[3])
32                step_size = int(args[4])
33               
34        if options.outfile != None:
35                out_file = open(options.outfile, 'w')
36               
37        #Generate snp and neutral bitsets
38        AR_snp_bitsets = binned_bitsets_from_file(open(snp_filename))
39        neutral_bitsets = binned_bitsets_from_file(open(neutral_filename))
40
41        # Generate divergence bitset from maf file
42        AR_div_bitsets = dict()
43        chr_lens = dict()
44        reader = bx.align.maf.Reader( open (maf_filename) )
45       
46        for block in reader:
47                comp1 = block.get_component_by_src_start( options.build )
48                comp2 = block.get_component_by_src_start( options.species )
49               
50                if comp1 is None or comp2 is None:
51                        continue
52                       
53                # Chromosome, start, and stop of reference species alignment
54                chr = comp1.src.split( '.' )[1]
55                start = comp1.start
56                end = comp1.end
57               
58                # Get or create bitset for this chromosome
59                if chr in AR_div_bitsets:
60                        bitset = AR_div_bitsets[chr]
61                else:
62                        bitset = AR_div_bitsets[chr] = bx.bitset.BinnedBitSet()
63                        chr_lens[chr] = comp1.get_src_size()
64                       
65                # Iterate over text and set diverged bit
66                pos = start
67                for ch1, ch2 in izip( comp1.text.upper(), comp2.text.upper() ):
68                        if ch1 == '-': continue
69                        if ch2 == '-':
70                                pos += 1
71                                continue
72                       
73                        if ch1 != ch2 and not AR_snp_bitsets[chr][pos]:
74                                bitset.set( pos )
75                        pos += 1
76                       
77        # Debugging Code
78#       for chr in AR_div_bitsets:
79#               for pos in range(0, AR_div_bitsets[chr].size):
80#                       if AR_div_bitsets[pos]:
81#                               print >> sys.stderr, chr, pos, pos+1
82       
83        # Copy div and snp bitsets
84        nonAR_snp_bitsets = dict()
85        for chr in AR_snp_bitsets:
86                nonAR_snp_bitsets[chr] = bx.bitset.BinnedBitSet()
87                nonAR_snp_bitsets[chr].ior(AR_snp_bitsets[chr])
88               
89        nonAR_div_bitsets = dict()
90        for chr in AR_div_bitsets:
91                nonAR_div_bitsets[chr] = bx.bitset.BinnedBitSet()
92                nonAR_div_bitsets[chr].ior(AR_div_bitsets[chr])
93       
94        # Generates AR snps by intersecting with neutral intervals
95        for chr in AR_snp_bitsets:
96                AR_snp_bitsets[chr].iand(neutral_bitsets[chr])
97       
98        # Generates AR divs by intersecting with neutral intervals     
99        for chr in AR_div_bitsets:
100                AR_div_bitsets[chr].iand(neutral_bitsets[chr])
101
102        # Inverts the neutral intervals so now represents nonAR
103        for chr in neutral_bitsets:
104                neutral_bitsets[chr].invert()
105       
106        # Generates nonAR snps by intersecting with masked neutral intervals   
107        for chr in nonAR_snp_bitsets:
108                nonAR_snp_bitsets[chr].iand(neutral_bitsets[chr])
109        # Generates nonAR divs by intersecting with masked neutral intervals           
110        for chr in nonAR_div_bitsets:
111                nonAR_div_bitsets[chr].iand(neutral_bitsets[chr])
112
113        for chr in AR_div_bitsets:
114                for window in range(0, chr_lens[chr] - window_size, step_size):
115#                       neutral_size = neutral_bitsets[chr].count_range(window, window_size)
116#                       if neutral_size < 9200: continue
117                        AR_snp = AR_snp_bitsets[chr].count_range(window, window_size)
118                        AR_div = AR_div_bitsets[chr].count_range(window, window_size)
119                        nonAR_snp = nonAR_snp_bitsets[chr].count_range(window, window_size)
120                        nonAR_div = nonAR_div_bitsets[chr].count_range(window, window_size)
121                       
122                        if nonAR_snp >= 6 and nonAR_div >= 6 and AR_snp >= 6 and AR_div >= 6:
123                                MK_pval = MK_chi_pvalue(nonAR_snp, nonAR_div, AR_snp, AR_div)
124                        else:
125                                MK_pval = MK_fisher_pvalue(nonAR_snp, nonAR_div, AR_snp, AR_div)
126                               
127                        if options.outfile != None:
128                                out_file.write("%s\t%d\t%d\t%d\t%d\t%d\t%d\t%1.15f\n" % (chr, window, window+window_size, nonAR_snp, nonAR_div, AR_snp, AR_div, MK_pval))
129                        else:
130                                print "%s\t%d\t%d\t%d\t%d\t%d\t%d\t%1.15f" % (chr, window, window+window_size, nonAR_snp, nonAR_div, AR_snp, AR_div, MK_pval)
131       
132        if options.outfile != None:
133                out_file.close()
134                       
135def MK_fisher_pvalue(win_snp, win_div, AR_snp, AR_div):
136
137        if win_snp == 0 and win_div == 0 and AR_snp == 0 and AR_div == 0:
138                return 1.0
139       
140        fisher_result = r.fisher_test(r.matrix(r.c([win_snp, win_div, AR_snp, AR_div]), nr = 2))
141       
142        return fisher_result['p.value']
143
144def MK_chi_pvalue(win_snp, win_div, AR_snp, AR_div):
145       
146        chi_result = r.chisq_test(r.matrix(r.c([win_snp, win_div, AR_snp, AR_div]), nr = 2))
147       
148        return chi_result['p.value']
149
150main()
Note: リポジトリブラウザについてのヘルプは TracBrowser を参照してください。