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
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
|
#!/usr/bin/python3
import sys
import bx.pwm.position_weight_matrix as pwmx
from bx.align import maf as align_maf
def isnan(x):
# return ieeespecial.isnan(x)
if x == x:
return False
return True
NaN = float("nan")
def main():
pwm_file = sys.argv[1]
splist = sys.argv[2]
if len(sys.argv) == 4:
inmaf = open(sys.argv[3])
else:
inmaf = sys.stdin
# read alignment species
species = []
for sp in splist.split(","):
species.append(sp)
# read weight matrices
pwm = {}
for wm in pwmx.Reader(open(pwm_file), format="basic"):
pwm[wm.id] = wm
fbunch = {}
for scoremax, index, headers in MafScorer(pwm, species, inmaf):
for k, matrix in scoremax.items():
fname = k + ".mx"
if fname not in fbunch:
fbunch[fname] = open(fname, "w")
print("Writing", fname, file=sys.stderr)
for i in range(len(matrix)):
for j in range(len(matrix[i])):
print(f"{matrix[i][j]:.2f}", end=" ", file=fbunch[fname])
print(file=fbunch[fname])
for file in fbunch.values():
file.close()
def MafScorer(pwm, species, inmaf):
index = 0
scoremax, width = None, None
for maf in align_maf.Reader(inmaf):
# try:
if True:
val = MafBlockScorer(pwm, species, maf)
for scoremax, width, headers in val:
yield scoremax, index, headers
try:
pass
except Exception:
print("Failed on:", file=sys.stderr)
syserr = align_maf.Writer(sys.stderr)
syserr.write(maf)
if width:
print(width, file=sys.stderr)
if scoremax:
print(len(scoremax), file=sys.stderr)
syserr.close()
sys.exit(1)
index += width
yield scoremax, index, headers
def MafMotifSelect(mafblock, pwm, motif=None, threshold=0):
if motif is not None and len(motif) != len(pwm):
raise Exception("pwm and motif must be the same length")
# generic alignment
alignlist = [c.text for c in mafblock.components]
align = pwmx.Align(alignlist)
nrows, ncols = align.dims
# required sequence length
minSeqLen = len(motif)
# record the text sizes from the alignment rows
for start in range(ncols - minSeqLen):
if align.rows[0][start] == "-":
continue
subseq = ""
pwm_score_vec = []
motif_score_vec = []
max_cols = 0
for ir in range(nrows):
expanded = align.rows[ir].count("-", start, minSeqLen)
subtext = align.rows[ir][start : minSeqLen + expanded]
max_cols = max(len(subtext), max_cols)
subseq = subtext.replace("-", "")
revseq = pwmx.reverse_complement(subseq)
# pwm score
nill, f_score = pwm.score_seq(subseq)[0]
r_score, nill = pwm.score_seq(revseq)[0]
pwm_score_vec.append(max(f_score, r_score))
# consensus score
if motif is not None:
for_score = int(pwmx.match_consensus(subseq, motif))
rev_score = int(pwmx.match_consensus(revseq, motif))
motif_score_vec.append(max(for_score, rev_score))
# check threshold
try:
assert not isnan(max(pwm_score_vec))
assert not isnan(max(motif_score_vec))
except AssertionError:
print(pwm_score_vec, motif_score_vec, file=sys.stderr)
print(len(subseq), len(pwm), file=sys.stderr)
if max(pwm_score_vec) < threshold:
continue
if max(motif_score_vec) < threshold:
continue
# chop block
col_start = start
col_end = max_cols + 1
motifmaf = mafblock.slice(col_start, col_end)
yield motifmaf, pwm_score_vec, motif_score_vec
"""
for ir in range(nrows):
# scan alignment row for motif subsequences
for start in range(ncols):
if align.rows[ir][start] == '-': continue
elif align.rows[ir][start] == 'n': continue
elif align.rows[ir][start] == 'N': continue
# gather enough subseq for motif
for ic in range(start,ncols):
char = align.rows[ir][ic].upper()
if char == '-' or char == 'N': continue
else: subseq += char
if len(subseq) == minSeqLen:
revseq = pwmx.reverse_complement( subseq )
align_match_lens.append( ic )
# pwm score
nill,f_score = pwm.score_seq( subseq )[0]
r_score, nill = pwm.score_seq( revseq )[0]
pwm_score_vec.append( max(f_score, r_score) )
# consensus score
if motif is not None:
for_score = int( pwmx.match_consensus(subseq,motif) )
rev_score = int( pwmx.match_consensus(revseq,motif) )
motif_score_vec.append( max(for_score, rev_score) )
#check threshold
try:
assert not isnan(max(pwm_score_vec) )
assert not isnan(max(motif_score_vec) )
except:
print >>sys.stderr, pwm_score_vec, motif_score_vec
print >>sys.stderr, len(subseq), len(pwm)
if max(pwm_score_vec) < threshold: continue
if max(motif_score_vec) < threshold: continue
# chop block
col_start = start
col_end = max( align_match_lens ) + 1
motifmaf = mafblock.slice( col_start, col_end )
print subseq,revseq,ic
print align_match_lens
yield motifmaf, pwm_score_vec, motif_score_vec
"""
def MafBlockScorer(pwm, species, maf):
width = len(maf.components[0].text)
headers = [(c.src, c.start, c.end) for c in maf.components]
# expand block rows to full
mafBlockSpecies = [specName.src.split(".")[0] for specName in maf.components]
alignlist = []
for sp in species:
try:
i = mafBlockSpecies.index(sp)
alignlist.append(maf.components[i].text)
except ValueError:
alignlist.append([NaN for n in range(width)])
alignrows = pwmx.Align(alignlist)
scoremax = {}
# record gap positions
filter = pwmx.score_align_gaps(alignrows)
# score pwm models
for model in pwm.keys():
scoremax[model] = pwm[model].score_align(alignrows, filter)
yield scoremax, width, headers
def MafMotifScorer(species, maf, motifs):
width = len(maf.components[0].text)
headers = [(c.src, c.start, c.end) for c in maf.components]
# expand block rows to full
mafBlockSpecies = [specName.src.split(".")[0] for specName in maf.components]
alignlist = []
for sp in species:
try:
i = mafBlockSpecies.index(sp)
alignlist.append(maf.components[i].text)
except ValueError:
alignlist.append([NaN for n in range(width)])
alignrows = pwmx.Align(alignlist, headers)
# record gap positions
filter = pwmx.score_align_gaps(alignrows)
# score motif
if isinstance(motifs, list):
scoremax = {}
for string in motifs:
scoremax[string] = pwmx.score_align_motif(alignrows, string, filter)
else:
scoremax = pwmx.score_align_motif(alignrows, motifs, filter)
yield scoremax, width, headers
if __name__ == "__main__":
main()
|
