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
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
|
# Copyright 2015 by Gert Hulselmans. All rights reserved.
# This file is part of the Biopython distribution and governed by your
# choice of the "Biopython License Agreement" or the "BSD 3-Clause License".
# Please see the LICENSE file that should have been included as part of this
# package.
"""Parse various position frequency matrix format files."""
import re
from Bio import motifs
class Record(list):
"""Class to store the information in a position frequency matrix table.
The record inherits from a list containing the individual motifs.
"""
def __str__(self):
return "\n".join(str(motif) for motif in self)
def read(handle, pfm_format):
"""Read motif(s) from a file in various position frequency matrix formats.
Return the record of PFM(s).
Call the appropriate routine based on the format passed.
"""
# Supporting underscores here for backward compatibility
pfm_format = pfm_format.lower().replace("_", "-")
if pfm_format == "pfm-four-columns":
record = _read_pfm_four_columns(handle)
return record
elif pfm_format == "pfm-four-rows":
record = _read_pfm_four_rows(handle)
return record
else:
raise ValueError("Unknown Position Frequency matrix format '%s'" % pfm_format)
def _read_pfm_four_columns(handle):
"""Read motifs in Cluster Buster position frequency matrix format from a file handle.
Cluster Buster motif format: http://zlab.bu.edu/cluster-buster/help/cis-format.html
#cisbp
Pos A C G T
1 0.00961538461538462 0.00961538461538462 0.00961538461538462 0.971153846153846
2 0.00961538461538462 0.00961538461538462 0.00961538461538462 0.971153846153846
3 0.971153846153846 0.00961538461538462 0.00961538461538462 0.00961538461538462
4 0.00961538461538462 0.00961538461538462 0.00961538461538462 0.971153846153846
5 0.00961538461538462 0.971153846153846 0.00961538461538462 0.00961538461538462
6 0.971153846153846 0.00961538461538462 0.00961538461538462 0.00961538461538462
7 0.00961538461538462 0.971153846153846 0.00961538461538462 0.00961538461538462
8 0.00961538461538462 0.00961538461538462 0.00961538461538462 0.971153846153846
#c2h2 zfs
Gene ENSG00000197372
Pos A C G T
1 0.341303 0.132427 0.117054 0.409215
2 0.283785 0.077066 0.364552 0.274597
3 0.491055 0.078208 0.310520 0.120217
4 0.492621 0.076117 0.131007 0.300256
5 0.250645 0.361464 0.176504 0.211387
6 0.276694 0.498070 0.197793 0.027444
7 0.056317 0.014631 0.926202 0.002850
8 0.004470 0.007769 0.983797 0.003964
9 0.936213 0.058787 0.002387 0.002613
10 0.004352 0.004030 0.002418 0.989200
11 0.013277 0.008165 0.001991 0.976567
12 0.968132 0.002263 0.002868 0.026737
13 0.397623 0.052017 0.350783 0.199577
14 0.000000 0.000000 1.000000 0.000000
15 1.000000 0.000000 0.000000 0.000000
16 0.000000 0.000000 1.000000 0.000000
17 0.000000 0.000000 1.000000 0.000000
18 1.000000 0.000000 0.000000 0.000000
19 0.000000 1.000000 0.000000 0.000000
20 1.000000 0.000000 0.000000 0.000000
#c2h2 zfs
Gene FBgn0000210
Motif M1734_0.90
Pos A C G T
1 0.25 0.0833333 0.0833333 0.583333
2 0.75 0.166667 0.0833333 0
3 0.833333 0 0 0.166667
4 1 0 0 0
5 0 0.833333 0.0833333 0.0833333
6 0.333333 0 0 0.666667
7 0.833333 0 0 0.166667
8 0.5 0 0.333333 0.166667
9 0.5 0.0833333 0.166667 0.25
10 0.333333 0.25 0.166667 0.25
11 0.166667 0.25 0.416667 0.166667
# flyfactorsurvey (cluster buster)
>AbdA_Cell_FBgn0000014
1 3 0 14
0 0 0 18
16 0 0 2
18 0 0 0
1 0 0 17
0 0 6 12
15 1 2 0
# homer
>ATGACTCATC AP-1(bZIP)/ThioMac-PU.1-ChIP-Seq(GSE21512)/Homer 6.049537 -1.782996e+03 0 9805.3,5781.0,3085.1,2715.0,0.00e+00
0.419 0.275 0.277 0.028
0.001 0.001 0.001 0.997
0.010 0.002 0.965 0.023
0.984 0.003 0.001 0.012
0.062 0.579 0.305 0.054
0.026 0.001 0.001 0.972
0.043 0.943 0.001 0.012
0.980 0.005 0.001 0.014
0.050 0.172 0.307 0.471
0.149 0.444 0.211 0.195
# hocomoco
> AHR_si
40.51343240527031 18.259112547756697 56.41253757072521 38.77363485291994
10.877470982533044 11.870876719950774 34.66312982331297 96.54723985087516
21.7165707818416 43.883079837598544 20.706746561638717 67.6523201955933
2.5465132509466635 1.3171620263517245 145.8637051322628 4.231336967110781
0.0 150.35847450464382 1.4927836298652875 2.1074592421627525
3.441039751299748 0.7902972158110341 149.37613720253387 0.3512432070271259
0.0 3.441039751299748 0.7024864140542533 149.81519121131782
0.0 0.0 153.95871737667187 0.0
43.07922333291745 66.87558226865211 16.159862546986584 27.844049228115868
# neph
UW.Motif.0001 atgactca
0.772949 0.089579 0.098612 0.038860
0.026652 0.004653 0.025056 0.943639
0.017663 0.023344 0.918728 0.040264
0.919596 0.025414 0.029759 0.025231
0.060312 0.772259 0.104968 0.062462
0.037406 0.020643 0.006667 0.935284
0.047316 0.899024 0.026928 0.026732
0.948639 0.019497 0.005737 0.026128
# tiffin
T A G C
30 0 28 40
0 0 0 99
0 55 14 29
0 99 0 0
20 78 0 0
0 52 7 39
19 46 11 22
0 60 38 0
0 33 0 66
73 0 25 0
99 0 0 0
"""
record = Record()
motif_name = None
motif_nbr = 0
motif_nbr_added = 0
default_nucleotide_order = ["A", "C", "G", "T"]
nucleotide_order = default_nucleotide_order
nucleotide_counts = {"A": [], "C": [], "G": [], "T": []}
for line in handle:
line = line.strip()
if line:
columns = line.split()
nbr_columns = len(columns)
if line.startswith("#"):
# Skip comment lines.
continue
elif line.startswith(">"):
# Parse ">AbdA_Cell_FBgn0000014" and "> AHR_si" like lines and put the part after ">" as motif name.
if motif_nbr != 0 and motif_nbr_added != motif_nbr:
# Add the previous motif to the record.
motif = motifs.Motif(alphabet="GATC", counts=nucleotide_counts)
motif.name = motif_name
record.append(motif)
motif_nbr_added = motif_nbr
# Reinitialize variables for the new motif.
motif_name = line[1:].strip()
nucleotide_order = default_nucleotide_order
elif columns[0] == "Gene":
# Parse "Gene ENSG00000197372" like lines and put the gene name as motif name.
if motif_nbr != 0 and motif_nbr_added != motif_nbr:
# Add the previous motif to the record.
motif = motifs.Motif(alphabet="GATC", counts=nucleotide_counts)
motif.name = motif_name
record.append(motif)
motif_nbr_added = motif_nbr
# Reinitialize variables for the new motif.
motif_name = columns[1]
nucleotide_order = default_nucleotide_order
elif columns[0] == "Motif":
# Parse "Motif M1734_0.90" like lines.
if motif_nbr != 0 and motif_nbr_added != motif_nbr:
# Add the previous motif to the record.
motif = motifs.Motif(alphabet="GATC", counts=nucleotide_counts)
motif.name = motif_name
record.append(motif)
motif_nbr_added = motif_nbr
# Reinitialize variables for the new motif.
motif_name = columns[1]
nucleotide_order = default_nucleotide_order
elif columns[0] == "Pos":
# Parse "Pos A C G T" like lines and change nucleotide order if necessary.
if nbr_columns == 5:
# If the previous line was not a "Gene ENSG00000197372" like line, a new motif starts here.
if motif_nbr != 0 and motif_nbr_added != motif_nbr:
# Add the previous motif to the record.
motif = motifs.Motif(alphabet="GATC", counts=nucleotide_counts)
motif.name = motif_name
record.append(motif)
motif_nbr_added = motif_nbr
nucleotide_order = default_nucleotide_order
if set(columns[1:]) == set(default_nucleotide_order):
nucleotide_order = columns[1:]
elif columns[0] in default_nucleotide_order:
# Parse "A C G T" like lines and change nucleotide order if necessary.
if nbr_columns == 4:
nucleotide_order = default_nucleotide_order
if set(columns) == set(default_nucleotide_order):
nucleotide_order = columns
else:
# Parse matrix columns lines and use the correct nucleotide order.
if nbr_columns == 4:
matrix_columns = columns
elif nbr_columns == 5:
matrix_columns = columns[1:]
else:
continue
if motif_nbr == motif_nbr_added:
# A new motif matrix starts here, so reinitialize variables for the new motif.
nucleotide_counts = {"A": [], "C": [], "G": [], "T": []}
motif_nbr += 1
[
nucleotide_counts[nucleotide].append(float(nucleotide_count))
for nucleotide, nucleotide_count in zip(
nucleotide_order, matrix_columns
)
]
else:
# Empty lines can be separators between motifs.
if motif_nbr != 0 and motif_nbr_added != motif_nbr:
# Add the previous motif to the record.
motif = motifs.Motif(alphabet="GATC", counts=nucleotide_counts)
motif.name = motif_name
record.append(motif)
motif_nbr_added = motif_nbr
# Reinitialize variables for the new motif.
motif_name = None
nucleotide_order = default_nucleotide_order
# nucleotide_counts = {'A': [], 'C': [], 'G': [], 'T': []}
if motif_nbr != 0 and motif_nbr_added != motif_nbr:
motif = motifs.Motif(alphabet="GATC", counts=nucleotide_counts)
motif.name = motif_name
record.append(motif)
return record
def _read_pfm_four_rows(handle):
"""Read motifs in position frequency matrix format from a file handle.
Cluster Buster motif format: http://zlab.bu.edu/cluster-buster/help/cis-format.html
#hdpi
A 0 5 6 5 1 0
C 1 1 0 0 0 4
G 5 0 0 0 3 0
T 0 0 0 1 2 2
# yetfasco
A 0.5 0.0 0.0 0.25 0.25 0.25 0.25 0.25 0.25 0.25 0.25 0.25 0.5 0.0 0.0833333334583333
T 0.0 0.0 0.0 0.25 0.25 0.25 0.25 0.25 0.25 0.25 0.25 0.25 0.0 0.0 0.0833333334583333
G 0.0 1.0 0.0 0.25 0.25 0.25 0.25 0.25 0.25 0.25 0.25 0.25 0.0 1.0 0.249999999875
C 0.5 0.0 1.0 0.25 0.25 0.25 0.25 0.25 0.25 0.25 0.25 0.25 0.5 0.0 0.583333333208333
#flyfactorsurvey ZFP finger
A | 92 106 231 135 0 1 780 28 0 700 739 94 60 127 130
C | 138 82 129 81 774 1 3 1 0 6 17 49 193 122 148
G | 270 398 54 164 7 659 1 750 755 65 1 41 202 234 205
T | 290 204 375 411 9 127 6 11 36 20 31 605 335 307 308
# scertf pcm
A | 9 1 1 97 1 94
T | 80 1 97 1 1 2
C | 9 97 1 1 1 2
G | 2 1 1 1 97 2
# scertf pfm
A | 0.090 0.010 0.010 0.970 0.010 0.940
C | 0.090 0.970 0.010 0.010 0.010 0.020
G | 0.020 0.010 0.010 0.010 0.970 0.020
T | 0.800 0.010 0.970 0.010 0.010 0.020
#idmmpmm
> abd-A
0.218451749734889 0.0230646871686108 0.656680805938494 0.898197242841994 0.040694591728526 0.132953340402969 0.74907211028632 0.628313891834571
0.0896076352067868 0.317338282078473 0.321580063626723 0.0461293743372216 0.0502386002120891 0.040694591728526 0.0284994697773065 0.0339342523860021
0.455991516436904 0.0691940615058324 0.0108695652173913 0.0217391304347826 0.0284994697773065 0.0284994697773065 0.016304347826087 0.160127253446448
0.235949098621421 0.590402969247084 0.0108695652173913 0.0339342523860021 0.880567338282079 0.797852598091198 0.206124072110286 0.17762460233298
# JASPAR
>MA0001.1 AGL3
A [ 0 3 79 40 66 48 65 11 65 0 ]
C [94 75 4 3 1 2 5 2 3 3 ]
G [ 1 0 3 4 1 0 5 3 28 88 ]
T [ 2 19 11 50 29 47 22 81 1 6 ]
or::
>MA0001.1 AGL3
0 3 79 40 66 48 65 11 65 0
94 75 4 3 1 2 5 2 3 3
1 0 3 4 1 0 5 3 28 88
2 19 11 50 29 47 22 81 1 6
"""
record = Record()
name_pattern = re.compile(r"^>\s*(.+)\s*")
row_pattern_with_nucleotide_letter = re.compile(
r"\s*([ACGT])\s*[[]*[|]*\s*([0-9.\s]+)\s*[]]*\s*"
)
row_pattern_without_nucleotide_letter = re.compile(r"\s*([0-9.\s]+)\s*")
motif_name = None
nucleotide_counts = {}
row_count = 0
nucleotides = ["A", "C", "G", "T"]
for line in handle:
line = line.strip()
name_match = name_pattern.match(line)
row_match_with_nucleotide_letter = row_pattern_with_nucleotide_letter.match(
line
)
row_match_without_nucleotide_letter = row_pattern_without_nucleotide_letter.match(
line
)
if name_match:
motif_name = name_match.group(1)
elif row_match_with_nucleotide_letter:
(nucleotide, counts_str) = row_match_with_nucleotide_letter.group(1, 2)
current_nucleotide_counts = counts_str.split()
nucleotide_counts[nucleotide] = [
float(current_nucleotide_count)
for current_nucleotide_count in current_nucleotide_counts
]
row_count += 1
if row_count == 4:
motif = motifs.Motif(alphabet="GATC", counts=nucleotide_counts)
if motif_name:
motif.name = motif_name
record.append(motif)
motif_name = None
nucleotide_counts = {}
row_count = 0
elif row_match_without_nucleotide_letter:
current_nucleotide_counts = row_match_without_nucleotide_letter.group(
1
).split()
nucleotide_counts[nucleotides[row_count]] = [
float(current_nucleotide_count)
for current_nucleotide_count in current_nucleotide_counts
]
row_count += 1
if row_count == 4:
motif = motifs.Motif(alphabet="GATC", counts=nucleotide_counts)
if motif_name:
motif.name = motif_name
record.append(motif)
motif_name = None
nucleotide_counts = {}
row_count = 0
return record
def write(motifs):
"""Return the representation of motifs in Cluster Buster position frequency matrix format."""
lines = []
for m in motifs:
line = f">{m.name}\n"
lines.append(line)
for ACGT_counts in zip(
m.counts["A"], m.counts["C"], m.counts["G"], m.counts["T"]
):
lines.append(
"{0:0.0f}\t{1:0.0f}\t{2:0.0f}\t{3:0.0f}\n".format(*ACGT_counts)
)
# Finished; glue the lines together.
text = "".join(lines)
return text
|
