annotate tools/ilmn_pacbio/cov_model.py @ 2:c2a356708570

Uploaded
author xuebing
date Fri, 09 Mar 2012 19:45:42 -0500
parents 9071e359b9a3
children
Ignore whitespace changes - Everywhere: Within whitespace: At end of lines:
rev   line source
0
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
1 #!/usr/bin/env python
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
2 from optparse import OptionParser, SUPPRESS_HELP
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
3 import os, random, quake
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
4
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
5 ############################################################
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
6 # cov_model.py
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
7 #
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
8 # Given a file of kmer counts, reports the cutoff to use
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
9 # to separate trusted/untrusted kmers.
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
10 ############################################################
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
11
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
12 ############################################################
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
13 # main
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
14 ############################################################
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
15 def main():
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
16 usage = 'usage: %prog [options] <counts file>'
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
17 parser = OptionParser(usage)
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
18 parser.add_option('--int', dest='count_kmers', action='store_true', default=False, help='Kmers were counted as integers w/o the use of quality values [default: %default]')
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
19 parser.add_option('--ratio', dest='ratio', type='int', default=200, help='Likelihood ratio to set trusted/untrusted cutoff [default: %default]')
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
20 parser.add_option('--no_sample', dest='no_sample', action='store_true', default=False, help='Do not sample kmer coverages into kmers.txt because its already done [default: %default]')
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
21 # help='Model kmer coverage as a function of GC content of kmers [default: %default]'
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
22 parser.add_option('--gc', dest='model_gc', action='store_true', default=False, help=SUPPRESS_HELP)
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
23 (options, args) = parser.parse_args()
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
24
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
25 if len(args) != 1:
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
26 parser.error('Must provide kmers counts file')
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
27 else:
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
28 ctsf = args[0]
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
29
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
30 if options.count_kmers:
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
31 model_cutoff(ctsf, options.ratio)
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
32 print 'Cutoff: %s' % open('cutoff.txt').readline().rstrip()
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
33
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
34 else:
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
35 if options.model_gc:
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
36 model_q_gc_cutoffs(ctsf, 25000, options.ratio)
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
37 else:
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
38 model_q_cutoff(ctsf, 50000, options.ratio, options.no_sample)
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
39 print 'Cutoff: %s' % open('cutoff.txt').readline().rstrip()
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
40
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
41
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
42 ############################################################
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
43 # model_cutoff
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
44 #
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
45 # Make a histogram of kmers to give to R to learn the cutoff
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
46 ############################################################
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
47 def model_cutoff(ctsf, ratio):
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
48 # make kmer histogram
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
49 cov_max = 0
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
50 for line in open(ctsf):
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
51 cov = int(line.split()[1])
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
52 if cov > cov_max:
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
53 cov_max = cov
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
54
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
55 kmer_hist = [0]*cov_max
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
56 for line in open(ctsf):
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
57 cov = int(line.split()[1])
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
58 kmer_hist[cov-1] += 1
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
59
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
60 cov_out = open('kmers.hist', 'w')
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
61 for cov in range(0,cov_max):
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
62 if kmer_hist[cov]:
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
63 print >> cov_out, '%d\t%d' % (cov+1,kmer_hist[cov])
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
64 cov_out.close()
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
65
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
66 os.system('R --slave --args %d < %s/cov_model.r 2> r.log' % (ratio,quake.quake_dir))
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
67
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
68
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
69 ############################################################
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
70 # model_q_cutoff
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
71 #
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
72 # Sample kmers to give to R to learn the cutoff
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
73 # 'div100' is necessary when the number of kmers is too
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
74 # large for random.sample, so we only consider every 100th
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
75 # kmer.
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
76 ############################################################
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
77 def model_q_cutoff(ctsf, sample, ratio, no_sample=False):
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
78 if not no_sample:
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
79 # count number of kmer coverages
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
80 num_covs = 0
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
81 for line in open(ctsf):
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
82 num_covs += 1
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
83
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
84 # choose random kmer coverages
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
85 div100 = False
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
86 if sample >= num_covs:
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
87 rand_covs = range(num_covs)
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
88 else:
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
89 if num_covs > 1000000000:
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
90 div100 = True
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
91 rand_covs = random.sample(xrange(num_covs/100), sample)
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
92 else:
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
93 rand_covs = random.sample(xrange(num_covs), sample)
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
94 rand_covs.sort()
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
95
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
96 # print to file
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
97 out = open('kmers.txt', 'w')
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
98 kmer_i = 0
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
99 rand_i = 0
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
100 for line in open(ctsf):
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
101 if div100:
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
102 if kmer_i % 100 == 0 and kmer_i/100 == rand_covs[rand_i]:
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
103 print >> out, line.split()[1]
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
104 rand_i += 1
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
105 if rand_i >= sample:
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
106 break
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
107 else:
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
108 if kmer_i == rand_covs[rand_i]:
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
109 print >> out, line.split()[1]
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
110 rand_i += 1
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
111 if rand_i >= sample:
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
112 break
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
113 kmer_i += 1
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
114 out.close()
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
115
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
116 os.system('R --slave --args %d < %s/cov_model_qmer.r 2> r.log' % (ratio,quake.quake_dir))
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
117
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
118
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
119 ############################################################
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
120 # model_q_gc_cutoffs
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
121 #
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
122 # Sample kmers to give to R to learn the cutoff for each
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
123 # GC value
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
124 ############################################################
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
125 def model_q_gc_cutoffs(ctsf, sample, ratio):
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
126 # count number of kmer coverages at each at
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
127 k = len(open(ctsf).readline().split()[0])
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
128 num_covs_at = [0]*(k+1)
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
129 for line in open(ctsf):
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
130 kmer = line.split()[0]
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
131 num_covs_at[count_at(kmer)] += 1
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
132
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
133 # for each AT bin
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
134 at_cutoffs = []
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
135 for at in range(1,k):
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
136 # sample covs
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
137 if sample >= num_covs_at[at]:
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
138 rand_covs = range(num_covs_at[at])
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
139 else:
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
140 rand_covs = random.sample(xrange(num_covs_at[at]), sample)
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
141 rand_covs.sort()
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
142
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
143 # print to file
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
144 out = open('kmers.txt', 'w')
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
145 kmer_i = 0
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
146 rand_i = 0
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
147 for line in open(ctsf):
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
148 (kmer,cov) = line.split()
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
149 if count_at(kmer) == at:
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
150 if kmer_i == rand_covs[rand_i]:
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
151 print >> out, cov
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
152 rand_i += 1
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
153 if rand_i >= sample:
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
154 break
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
155 kmer_i += 1
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
156 out.close()
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
157
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
158 os.system('R --slave --args %d < %s/cov_model_qmer.r 2> r%d.log' % (ratio,quake.quake_dir,at))
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
159
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
160 at_cutoffs.append( open('cutoff.txt').readline().rstrip() )
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
161 if at in [1,k-1]: # setting extremes to next closests
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
162 at_cutoffs.append( open('cutoff.txt').readline().rstrip() )
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
163
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
164 os.system('mv kmers.txt kmers.at%d.txt' % at)
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
165 os.system('mv cutoff.txt cutoff.at%d.txt' % at)
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
166
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
167 out = open('cutoffs.gc.txt','w')
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
168 print >> out, '\n'.join(at_cutoffs)
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
169 out.close()
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
170
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
171
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
172 ############################################################
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
173 # model_q_gc_cutoffs_bigmem
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
174 #
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
175 # Sample kmers to give to R to learn the cutoff for each
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
176 # GC value
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
177 ############################################################
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
178 def model_q_gc_cutoffs_bigmem(ctsf, sample, ratio):
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
179 # input coverages
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
180 k = 0
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
181 for line in open(ctsf):
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
182 (kmer,cov) = line.split()
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
183 if k == 0:
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
184 k = len(kmer)
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
185 at_covs = ['']*(k+1)
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
186 else:
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
187 at = count_at(kmer)
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
188 if at_covs[at]:
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
189 at_covs[at].append(cov)
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
190 else:
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
191 at_covs[at] = [cov]
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
192
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
193 for at in range(1,k):
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
194 print '%d %d' % (at,len(at_covs[at]))
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
195
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
196 # for each AT bin
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
197 at_cutoffs = []
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
198 for at in range(1,k):
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
199 # sample covs
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
200 if sample >= len(at_covs[at]):
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
201 rand_covs = at_covs[at]
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
202 else:
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
203 rand_covs = random.sample(at_covs[at], sample)
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
204
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
205 # print to file
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
206 out = open('kmers.txt', 'w')
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
207 for rc in rand_covs:
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
208 print >> out, rc
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
209 out.close()
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
210
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
211 os.system('R --slave --args %d < %s/cov_model_qmer.r 2> r%d.log' % (ratio,quake.quake_dir,at))
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
212
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
213 at_cutoffs.append( open('cutoff.txt').readline().rstrip() )
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
214 if at in [1,k-1]: # setting extremes to next closests
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
215 at_cutoffs.append( open('cutoff.txt').readline().rstrip() )
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
216
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
217 os.system('mv kmers.txt kmers.at%d.txt' % at)
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
218 os.system('mv cutoff.txt cutoff.at%d.txt' % at)
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
219
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
220 out = open('cutoffs.gc.txt','w')
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
221 print >> out, '\n'.join(at_cutoffs)
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
222 out.close()
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
223
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
224
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
225 ############################################################
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
226 # count_at
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
227 #
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
228 # Count A's and T's in the given sequence
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
229 ############################################################
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
230 def count_at(seq):
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
231 return len([nt for nt in seq if nt in ['A','T']])
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
232
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
233
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
234 ############################################################
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
235 # __main__
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
236 ############################################################
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
237 if __name__ == '__main__':
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
238 main()