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
|
from unittest import TestCase
from numpy.testing import assert_allclose, assert_equal
from cogent3.parse.blast import (
FastacmdTaxonomyParser,
GenericBlastParser9,
LastProteinIds9,
PsiBlastFinder,
PsiBlastParser9,
PsiBlastQueryFinder,
PsiBlastTableParser,
QMEBlast9,
QMEPsiBlast9,
TableToValues,
fastacmd_taxonomy_splitter,
is_blast_junk,
is_blat_junk,
iter_finder,
iteration_set_finder,
make_label,
query_finder,
)
class BlastTests(TestCase):
"""Tests of top-level functions"""
def setUp(self):
"""Define some standard data"""
self.rec = """# BLASTP 2.2.10 [Oct-19-2004]
# Iteration: 1
# Query: ece:Z4181
# Database: db/everything.faa
# Fields: Query id, Subject id, % identity, alignment length, mismatches, gap openings, q. start, q. end, s. start, s. end, e-value, bit score
ece:Z4181 ece:Z4181 100.00 110 0 0 1 110 1 110 3e-47 187
ece:Z4181 ecs:ECs3717 100.00 110 0 0 1 110 1 110 3e-47 187
ece:Z4181 cvi:CV2421 41.67 72 42 0 39 110 29 100 2e-06 52.8
# BLASTP 2.2.10 [Oct-19-2004]
# Iteration: 2
# Query: ece:Z4181
# Database: db/everything.faa
# Fields: Query id, Subject id, % identity, alignment length, mismatches, gap openings, q. start, q. end, s. start, s. end, e-value, bit score
ece:Z4181 ece:Z4181 100.00 110 0 0 1 110 1 110 3e-54 211
ece:Z4181 ecs:ECs3717 100.00 110 0 0 1 110 1 110 3e-54 211
ece:Z4181 cvi:CV2421 41.67 72 42 0 39 110 29 100 2e-08 59.0
ece:Z4181 sfl:CP0138 33.98 103 57 2 8 110 6 97 6e-06 50.5
ece:Z4181 spt:SPA2730 37.50 72 45 0 39 110 30 101 1e-05 49.8
ece:Z4181 sec:SC2804 37.50 72 45 0 39 110 30 101 1e-05 49.8
ece:Z4181 stm:STM2872 37.50 72 45 0 39 110 30 101 1e-05 49.8""".split(
"\n"
)
self.rec2 = """# BLASTP 2.2.10 [Oct-19-2004]
# Iteration: 1
# Query: ece:Z4181
# Database: db/everything.faa
# Fields: Query id, Subject id, % identity, alignment length, mismatches, gap openings, q. start, q. end, s. start, s. end, e-value, bit score
ece:Z4181 ece:Z4181 100.00 110 0 0 1 110 1 110 3e-47 187
ece:Z4181 ecs:ECs3717 100.00 110 0 0 1 110 1 110 3e-47 187
ece:Z4181 cvi:CV2421 41.67 72 42 0 39 110 29 100 2e-06 52.8
# BLASTP 2.2.10 [Oct-19-2004]
# Iteration: 2
# Query: ece:Z4181
# Database: db/everything.faa
# Fields: Query id, Subject id, % identity, alignment length, mismatches, gap openings, q. start, q. end, s. start, s. end, e-value, bit score
ece:Z4181 ece:Z4181 100.00 110 0 0 1 110 1 110 3e-54 211
ece:Z4181 ecs:ECs3717 100.00 110 0 0 1 110 1 110 3e-54 211
ece:Z4181 cvi:CV2421 41.67 72 42 0 39 110 29 100 2e-08 59.0
ece:Z4181 sfl:CP0138 33.98 103 57 2 8 110 6 97 6e-06 50.5
ece:Z4181 spt:SPA2730 37.50 72 45 0 39 110 30 101 1e-05 49.8
ece:Z4181 sec:SC2804 37.50 72 45 0 39 110 30 101 1e-05 49.8
ece:Z4181 stm:STM2872 37.50 72 45 0 39 110 30 101 1e-05 49.8
# BLASTP 2.2.10 [Oct-19-2004]
# Iteration: 1
# Query: ece:Z4182
# Database: db/everything.faa
# Fields: Query id, Subject id, % identity, alignment length, mismatches, gap openings, q. start, q. end, s. start, s. end, e-value, bit score
ece:Z4182 ece:Z4182 100.00 110 0 0 1 110 1 110 3e-47 187
ece:Z4182 ecs:ECs3718 100.00 110 0 0 1 110 1 110 3e-47 187
ece:Z4182 cvi:CV2422 41.67 72 42 0 39 110 29 100 2e-06 52.8""".split(
"\n"
)
self.rec3 = """# BLASTP 2.2.10 [Oct-19-2004]
# Iteration: 1
# Query: ece:Z4181
# Database: db/everything.faa
# Fields: Query id, Subject id, % identity, alignment length, mismatches, gap openings, q. start, q. end, s. start, s. end, e-value, bit score
ece:Z4181 ece:Z4181 100.00 110 0 0 1 110 1 110 3e-47 187
ece:Z4181 ecs:ECs3717 100.00 110 0 0 1 110 1 110 3e-47 187
ece:Z4181 spt:SPA2730 37.50 72 45 0 39 110 30 101 1e-05 49.8
# BLASTP 2.2.10 [Oct-19-2004]
# Iteration: 2
# Query: ece:Z4181
# Database: db/everything.faa
# Fields: Query id, Subject id, % identity, alignment length, mismatches, gap openings, q. start, q. end, s. start, s. end, e-value, bit score
ece:Z4181 ecs:ECs3717 100.00 110 0 0 1 110 1 110 3e-54 211
ece:Z4181 cvi:CV2421 41.67 72 42 0 39 110 29 100 2e-08 59.0
# BLASTP 2.2.10 [Oct-19-2004]
# Iteration: 1
# Query: ece:Z4182
# Database: db/everything.faa
# Fields: Query id, Subject id, % identity, alignment length, mismatches, gap openings, q. start, q. end, s. start, s. end, e-value, bit score
ece:Z4182 ece:Z4182 100.00 110 0 0 1 110 1 110 3e-47 187
ece:Z4182 cvi:CV2422 41.67 72 42 0 39 110 29 100 2e-06 52.8""".split(
"\n"
)
def test_iter_finder(self):
"""iter_finder should split on lines starting with '# Iteration:'"""
lines = "abc\n# Iteration: 3\ndef".splitlines()
self.assertEqual(list(map(iter_finder, lines)), [False, True, False])
def test_query_finder(self):
"""query_finder should split on lines starting with '# Query:'"""
lines = "abc\n# Query: dfdsffsd\ndef".split("\n")
self.assertEqual(list(map(query_finder, lines)), [False, True, False])
def test_iteration_set_finder(self):
"""iter_finder should split on lines starting with '# Iteration:'"""
lines = "abc\n# Iteration: 3\ndef\n# Iteration: 1".split("\n")
self.assertEqual(
list(map(iteration_set_finder, lines)), [False, False, False, True]
)
def test_is_junk(self):
"""is_junk should reject an assortment of invalid lines"""
# Note: testing two functions that call it instead of function itself
lines = "abc\n# BLAST blah blah\n \n# BLAT blah\n123".split("\n")
self.assertEqual(
list(map(is_blast_junk, lines)), [False, True, True, False, False]
)
self.assertEqual(
list(map(is_blat_junk, lines)), [False, False, True, True, False]
)
def test_make_label(self):
"""make_label should turn comment lines into (key, val) pairs"""
a = "this test will fail: no # at start"
b = "#this test will fail because no colon"
c = "# Iteration: 1"
d = "# Query: ece:Z4147 ygdP; putative invasion protein [EC:3.6.1.-]"
e = "#Iteration: 1" # no space after the hash
self.assertRaises(ValueError, make_label, a)
self.assertRaises(ValueError, make_label, b)
# Note that we _do_ map the data type of known values value, so the
# value of the iteration will be 1, not '1'
self.assertEqual(make_label(c), ("ITERATION", 1))
self.assertEqual(
make_label(d),
("QUERY", "ece:Z4147 ygdP; putative invasion protein [EC:3.6.1.-]"),
)
self.assertEqual(make_label(e), ("ITERATION", 1))
def test_TableToValues(self):
"""TableToValues should convert itself into the correct type."""
constructors = {"a": int, "b": float, "c": str}
table = [["c", "b", "a", "d"], ["1.5", "3.5", "2", "2.5"], ["1", "2", "3", "4"]]
self.assertEqual(
TableToValues(table, constructors),
([["1.5", 3.5, 2, "2.5"], ["1", 2.0, 3, "4"]], ["c", "b", "a", "d"]),
)
# check that it works with supplied header
self.assertEqual(
TableToValues(table[1:], constructors, list("cbad")),
([["1.5", 3.5, 2, "2.5"], ["1", 2.0, 3, "4"]], ["c", "b", "a", "d"]),
)
def test_PsiBlastTableParser(self):
"""PsiBlastTableParser should wrap values in table."""
fields = [
v.strip()
for v in "Query id, Subject id, % identity, alignment length, mismatches, gap openings, q. start, q. end, s. start, s. end, e-value, bit score".split(
","
)
]
table = [
v.split()
for v in """ece:Z4147 ece:Z4147 100.00 176 0 0 1 176 1 176 2e-89 328
ece:Z4147 ecs:ECs3687 100.00 176 0 0 1 176 1 176 2e-89 328
ece:Z4147 ecc:c3425 100.00 176 0 0 1 176 1 176 2e-89 328
ece:Z4147 sfl:SF2840 100.00 176 0 0 1 176 1 176 2e-89 328""".splitlines()
]
headed_table = [fields] + table
new_table, new_fields = PsiBlastTableParser(headed_table)
self.assertEqual(new_fields, fields)
self.assertEqual(len(new_table), 4)
self.assertEqual(
new_table[1],
["ece:Z4147", "ecs:ECs3687", 100.0, 176, 0, 0, 1, 176, 1, 176, 2e-89, 328],
)
def test_GenericBlastParser9(self):
"""GenericBlastParser9 should read blast's tabular format (#9)."""
rec = self.rec
p = GenericBlastParser9(rec, PsiBlastFinder)
result = list(p)
self.assertEqual(len(result), 2)
first, second = result
self.assertEqual(
first[0],
{
"ITERATION": 1,
"QUERY": "ece:Z4181",
"DATABASE": "db/everything.faa",
"FIELDS": "Query id, Subject id, % identity, alignment length, mismatches, gap openings, q. start, q. end, s. start, s. end, e-value, bit score",
},
)
self.assertEqual(len(first[1]), 3)
self.assertEqual(second[0]["ITERATION"], 2)
self.assertEqual(len(second[1]), 7)
self.assertEqual(
second[1][-1],
"ece:Z4181 stm:STM2872 37.50 72 45 0 39 110 30 101 1e-05 49.8".split(),
)
def test_PsiBlastParser9(self):
"""PsiBlastParser9 should provide convenient results for format #9."""
result = PsiBlastParser9(self.rec2)
self.assertEqual(len(result), 2)
assert "ece:Z4181" in result
assert "ece:Z4182" in result
first = result["ece:Z4181"]
second = result["ece:Z4182"]
self.assertEqual(len(first), 2)
self.assertEqual(len(second), 1)
iter_1 = first[0]
iter_2 = first[1]
self.assertEqual(len(iter_1), 3)
self.assertEqual(len(iter_2), 7)
iter_1_2 = second[0]
self.assertEqual(len(iter_1_2), 3)
self.assertEqual(len(result["ece:Z4181"][1][3]), 12)
self.assertEqual(result["ece:Z4181"][1][3]["ALIGNMENT LENGTH"], 103)
def test_LastProteinIds9(self):
"""LastProteinIds9 should give last protein ids in iter"""
result = LastProteinIds9(self.rec)
self.assertEqual(
result,
[
"ece:Z4181",
"ecs:ECs3717",
"cvi:CV2421",
"sfl:CP0138",
"spt:SPA2730",
"sec:SC2804",
"stm:STM2872",
],
)
# should also work if threshold set
result = LastProteinIds9(self.rec, False, threshold=8e-6)
self.assertEqual(
result, ["ece:Z4181", "ecs:ECs3717", "cvi:CV2421", "sfl:CP0138"]
)
# should work on multiple records
result = list(map(LastProteinIds9, PsiBlastQueryFinder(self.rec2)))
self.assertEqual(len(result), 2)
self.assertEqual(
result[0],
[
"ece:Z4181",
"ecs:ECs3717",
"cvi:CV2421",
"sfl:CP0138",
"spt:SPA2730",
"sec:SC2804",
"stm:STM2872",
],
)
self.assertEqual(result[1], ["ece:Z4182", "ecs:ECs3718", "cvi:CV2422"])
def test_QMEBlast9(self):
"""QMEBlast9 should return expected lines from all iterations"""
expect = list(
zip(
*[
("ece:Z4181", "ece:Z4181", 3e-47),
("ece:Z4181", "ecs:ECs3717", 3e-47),
("ece:Z4181", "spt:SPA2730", 1e-5),
("ece:Z4181", "ecs:ECs3717", 3e-54), # WARNING: allows duplicates
("ece:Z4181", "cvi:CV2421", 2e-8),
("ece:Z4182", "ece:Z4182", 3e-47),
("ece:Z4182", "cvi:CV2422", 2e-6),
],
)
)
got = list(zip(*QMEBlast9(self.rec3)))
assert_equal(got[:-1], expect[:-1])
assert_allclose(got[-1], expect[-1])
def test_QMEPsiBlast9(self):
"""QMEPsiBlast9 should only return items from last iterations"""
expect = list(
zip(
*[
("ece:Z4181", "ecs:ECs3717", 3e-54),
("ece:Z4181", "cvi:CV2421", 2e-8),
("ece:Z4182", "ece:Z4182", 3e-47),
("ece:Z4182", "cvi:CV2422", 2e-6),
]
)
)
got = list(zip(*QMEPsiBlast9(self.rec3)))
assert_equal(got[:-1], expect[:-1])
assert_allclose(got[-1], expect[-1])
def test_fastacmd_taxonomy_splitter(self):
"""fastacmd_taxonomy_splitter should split records into groups"""
text = """NCBI sequence id: gi|3021565|emb|AJ223314.1|PSAJ3314
NCBI taxonomy id: 3349
Common name: Scots pine
Scientific name: Pinus sylvestris
NCBI sequence id: gi|37777029|dbj|AB108787.1|
NCBI taxonomy id: 228610
Common name: cf. Acremonium sp. KR21-2
Scientific name: cf. Acremonium sp. KR21-2
""".splitlines()
recs = list(fastacmd_taxonomy_splitter(text))
self.assertEqual(len(recs), 2)
self.assertEqual(recs[0], text[:5]) # includes trailing blank
def test_FastaCmdTaxonomyParser(self):
"""FastaCmdTaxonomyParser should parse taxonomy record to dict"""
text = """NCBI sequence id: gi|3021565|emb|AJ223314.1|PSAJ3314
NCBI taxonomy id: 3349
Common name: Scots pine
Scientific name: Pinus sylvestris
NCBI sequence id: gi|37777029|dbj|AB108787.1|
NCBI taxonomy id: 228610
Common name: cf. Acremonium sp. KR21-2
Scientific name: cf. Acremonium sp. KR21-2
""".splitlines()
recs = list(FastacmdTaxonomyParser(text))
self.assertEqual(len(recs), 2)
for r in recs:
self.assertEqual(
sorted(r.keys()), ["common_name", "scientific_name", "seq_id", "tax_id"]
)
r0, r1 = recs
self.assertEqual(r0["tax_id"], "3349")
self.assertEqual(r0["common_name"], "Scots pine")
self.assertEqual(r0["scientific_name"], "Pinus sylvestris")
self.assertEqual(r0["seq_id"], "gi|3021565|emb|AJ223314.1|PSAJ3314")
self.assertEqual(r1["tax_id"], "228610")
|
