1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
|
#!/usr/bin/env python
import sys
import os
from multiprocessing import Pool
from subprocess import *
import numpy as np
import pdb
fg_mers = {}
bg_mers = {}
if(len(sys.argv) == 5):
selectivity_fn = sys.argv[1]
fg_fasta_fn = sys.argv[2]
bg_fasta_fn = sys.argv[3]
output_file = sys.argv[4]
else:
print "please specify your inputs"
print "ex: select_mers.py fg_counts_file fg_fasta_file bg_counts_file bg_fasta_file output_file"
exit()
# empty class to fill up mer information with
class Mer:
pass
# import our variables
min_mer_range = int(os.environ.get("min_mer_range", 6));
max_mer_range = int(os.environ.get("max_mer_range", 10));
min_mer_count = int(os.environ.get("min_mer_count", 0));
max_select = int(os.environ.get("max_select", 15));
max_mer_distance = int(os.environ.get("max_mer_distance", 5000));
def populate_locations(input_fn, mers, mer):
''' Run the strstreamone command, and parse in the integers that are output
by the command, and add it to mers[mer].pts
'''
cmd = 'strstreamone ' + mer + " < " + input_fn
strstream = Popen(cmd, stdout=PIPE, shell=True)
for line in strstream.stdout:
mers[mer].pts.append(int(line))
def score_mers(selected):
from itertools import combinations
import time
scores = []
p = Pool()
fh = open(output_file, 'w');
fh.write("scores:\n");
for select_n in range(1, max_select+1):
print "scoring size ", select_n,
t = time.time()
scores_it = p.imap_unordered(score, combinations(selected, select_n))
for score_res in scores_it:
fh.write(str(score_res) + "\n");
print "size ", select_n, "took:", t - time.time()
return scores
def score(combination):
# input is a string of mers like
# ['ACCAA', 'ACCCGA', 'ACGTATA']
ret = [combination]
for mer in combination:
for other_mer in combination:
if not mer == other_mer:
if mer in other_mer:
ret.append("duplicates")
return ret
fg_pts = []
fg_dist = []
bg_pts = []
bg_dist = []
for mer in combination:
fg_pts = fg_pts + fg_mers[mer].pts
bg_pts = bg_pts + bg_mers[mer].pts
fg_pts.sort()
bg_pts.sort()
# remove any exact duplicates
# fg_pts = list(set(fg_pts))
# bg_pts = list(set(bg_pts))
# distances
min_mer_distance = max(len(i) for i in combination)
fg_dist = np.array([abs(fg_pts[i] - fg_pts[i-1]) for i in range(1, len(fg_pts))])
bg_dist = np.array([abs(bg_pts[i] - bg_pts[i-1]) for i in range(1, len(bg_pts))])
# return without calculating scores if any objects are higher than our max distance
if any(dist > max_mer_distance for dist in fg_dist):
ret.append("max")
ret.append(max(fg_dist))
return ret
# return without calculating scores if any mers are closer than the length of our longest mer in the combination
if any(dist < min_mer_distance for dist in fg_dist):
ret.append("min")
ret.append(min(fg_dist))
return ret
nb_primers = len(combination)
fg_mean_dist = np.mean(fg_dist)
fg_variance_dist = np.var(fg_dist)
bg_mean_dist = np.mean(bg_dist)
bg_variance_dist = np.var(bg_dist)
# this is our equation
score = (nb_primers * fg_mean_dist * fg_variance_dist) / ((bg_mean_dist * bg_variance_dist) + .000001)
ret.append(score)
ret.append(fg_mean_dist)
ret.append(fg_variance_dist)
ret.append(bg_mean_dist)
ret.append(bg_variance_dist)
return ret
def pop_fg(mer):
''' helper for map function '''
populate_locations(fg_fasta_fn, fg_mers, mer)
def pop_bg(mer):
''' helper for map function '''
populate_locations(bg_fasta_fn, bg_mers, mer)
def main():
import time
selected = []
selectivity_fh = open(selectivity_fn, "r")
# get our genome length
fg_genome_length = os.path.getsize(fg_fasta_fn)
bg_genome_length = os.path.getsize(bg_fasta_fn)
for row in selectivity_fh:
(mer, fg_count, bg_count, selectivity) = row.split()
fg_mers[mer] = Mer()
fg_mers[mer].pts = []
fg_mers[mer].count = fg_count
bg_mers[mer] = Mer()
bg_mers[mer].pts = []
bg_mers[mer].count = bg_count
selected.append([mer, selectivity])
# exhaustive = False
#
# if exhaustive:
# selected = fg_mers.keys()
# else:
# selected = select_mers(fg_mers, bg_mers, max_select)
selected = selected[-100:]
selected_mers = [row[0] for row in selected]
pdb.set_trace()
# print "searching through combinations of"
# print selected
print "Populating foreground locations"
map(pop_fg, selected_mers)
map(pop_bg, selected_mers)
scores = score_mers(selected_mers)
print "fg_genome_length", fg_genome_length
print "bg_genome_length", bg_genome_length
print "output_file:", output_file
if __name__ == "__main__":
sys.exit(main())
|
