Dav*_*K11 14 python optimization for-loop list-comprehension
我写了一个函数来计算某些字符(的出现次数A,C,G并T在同一位置的多个字符串中),并保存在字典中出现的次数.
例如,使用这两个字符串'ACGG'和'CAGT',它应该返回:
{'A': [1, 1, 0, 0], 'C': [1, 1, 0, 0], 'G': [0, 0, 2, 1], 'T': [0, 0, 0, 1]}
我想将下面的代码转换为列表理解以优化速度.它使用两个嵌套的for循环,输入Motifs是一个包含A的C的G和T的字符串列表.
def CountWithPseudocounts(Motifs):
count = {}
k = len(Motifs[0])
t = len(Motifs)
for s in 'ACGT':
count[s] = [0] * k
for i in range(t):
for j in range(k):
symbol = Motifs[i][j]
count[symbol][j] += 1
return count
我已经尝试替换函数底部的嵌套for循环以获得此列表理解:
count = [ [ count[Motifs[i][j]][j] += 1 ] for i in range(0, t) ] for j in range(0, k)]
它不起作用,可能是因为我不允许在列表理解中进行+ = 1的值赋值.我该如何解决这个问题？
Kas*_*mvd 11
你可以使用zip():
In [10]: a = 'ACGG'
In [11]: b = 'CAGT'
In [12]: chars = ['A', 'C', 'G', 'T']
In [13]: [[(ch==i) + (ch==j) for i, j in zip(a, b)] for ch in chars]
Out[13]: [[1, 1, 0, 0], [1, 1, 0, 0], [0, 0, 2, 1], [0, 0, 0, 1]]
如果你想要一本字典,你可以使用字典理解:
In [25]: {ch:[(ch==i) + (ch==j) for i, j in zip(a, b)] for ch in chars}
Out[25]: {'T': [0, 0, 0, 1], 'G': [0, 0, 2, 1], 'C': [1, 1, 0, 0], 'A': [1, 1, 0, 0]}
或者,如果您希望结果与字符列表的顺序相同,则可以使用collections.OrderedDict:
In [26]: from collections import OrderedDict
In [27]: OrderedDict((ch, [(ch==i) + (ch==j) for i, j in zip(a, b)]) for ch in chars)
Out[28]: OrderedDict([('A', [1, 1, 0, 0]), ('C', [1, 1, 0, 0]), ('G', [0, 0, 2, 1]), ('T', [0, 0, 0, 1])])
如果你仍然需要更多的性能和/或你正在处理长字符串和更大的数据集,你可以使用Numpy通过矢量化方法来解决这个问题.
In [61]: pairs = np.array((list(a), list(b))).T
In [62]: chars
Out[62]:
array(['A', 'C', 'G', 'T'],
dtype='<U1')
In [63]: (chars[:,None,None] == pairs).sum(2)
Out[63]:
array([[1, 1, 0, 0],
[1, 1, 0, 0],
[0, 0, 2, 1],
[0, 0, 0, 1]])
你确实不能在列表理解中做任务(你可以 - 通过调用函数 - 执行副作用).列表理解需要一个表达式.此外,您想要分配count并同时更新旧的很奇怪count.
使用不太有效的字典理解和列表理解来实现此目的的方法是:
chars = 'ACGT'
a = 'ACGG'
b = 'CAGT'
sequences = list(zip(a,b))
counts = {char:[seq.count(char) for seq in sequences] for char in chars}
(归功于@Chris_Rands的seq.count(char)建议)
这会产生:
{'G': [0, 0, 2, 1], 'A': [1, 1, 0, 0], 'C': [1, 1, 0, 0], 'T': [0, 0, 0, 1]}
通过调用更多字符串,您可以轻松地概括解决方案以计算zip(..)更多字符串.
您还可以决定优化算法本身.这可能会更有效,因为那时你只需要遍历字符串一次就可以使用字典的查找,例如:
def CountWithPseudocounts(sequences):
k = len(sequences[0])
count = {char:[0]*k for char in 'ACGT'}
for sequence in sequences:
j = 0
for symbol in sequence:
count[symbol][j] += 1
j += 1
return count
编辑:
如果要为计数中的所有元素添加一个,可以使用:
counts = {char:[seq.count(char)+1 for seq in sequences] for char in chars}如果速度表现真的很重要,我建议@Kasramvs提供的numpy方法.
此外,即使对于现代计算机来说,计算字符也不友好,也许你可以在计数之前使用一些关于输入字符串的索引/散列的技巧.例如,因为每个字符串只有4个字符,每个字符只包含4个可能的字母,'A','C','G','T',因此它可以很容易地代表所有'ACGT'组合从"AAAA"到"TTTT",带有数字,哈希或神秘的代码.这里组合的体积应等于或小于4x4x4x4 = 256个不同的数字.
然后计算代码.例如,每次看到'AAAA'然后在python列表或numpy数组中将其计为0x0时,请参阅'AAAC'并计为0x1,反之亦然.之后,你将获得一个索引范围为0x0~0xFF(255)的binning数组以及相应的出现,对吧？现在记住,一个0x0代表A:{1,1,1,1}在你的情况下,或七个0x1代表A:{7,7,7,0}以及C:{0,0,0,7} ...仔细总结它们,这就是结果.
在这种情况下,这些技巧应该可以帮助提高速度性能.速度有两个因素:第一个是现在计算机处理数字而不是字符,而数字比字符更容易排序,计数,分组,分区或索引; 第二个来自自然界更高的缓存命中率,因为这些技巧中的内存跟踪减少了很多.
我希望这可能有所帮助.:)
那么,添加一些代码如下是一个明确的参考.
首先,进口:
import itertools
import numpy as np
并且encode_sequence()下面的函数应该足够快,前提t02是dict 不是太大,比如通常少于1M的键值对:
def encode_sequence(tsq):
t00 = itertools.product('ACGT', repeat=4)
t01 = np.array(list(t00)).view('|U4').ravel()
t02 = dict(zip(t01, np.arange(len(t01))))
t03 = [t02[i] for i in tsq]
return t03
并使用下面的代码片段生成一个张量map_decode,代表约``统计代码"东西" ......此外,还有称为增强矩阵变换这部分下方的数学技巧,称它转换ll0和'ACGT'以同样的方式来跨度map_decode供以后使用.
ll0 = np.array([[1, 0, 0, 0], [0, 1, 0, 0], [0, 0, 1, 0], [0, 0, 0, 1]])
map_decode = np.array(list(itertools.product(ll0, repeat=4)))
喂一个测试序列并翻译,
test_seq = ('ACGG', 'CAGT', 'CCGA', 'AAAT', 'TTGC', 'GCAT', 'ACGG', 'AAAT')
u01 = encode_sequence(test_seq)
计算出现次数; 请注意下面的块应该是速度增益的主要来源,因为计算机擅长处理数字u01,
p01, q01 = np.unique(u01, return_counts=True)
毕竟,生成输出...这是一个有点棘手这里,例如p01是排序的散列码的test_seq和q01确实是相应的计数,而map_decode作为我说什么,一个张量在映射的哈希码p01我们想要另一个向量例如,将0x0(或'AAAA')映射到A:[1,1,1,1],C:[0,0,0,0],G:[0,0,0,0]和T: [0,0,0,0].因此,映射map_decode[p01]由计数加权q01并准备为报告求和:
np.sum(map_decode[p01]*q01[:, None, None], axis=0).T
它说,
array([[4, 3, 3, 1],
[2, 4, 0, 1],
[1, 0, 5, 2],
[1, 1, 0, 4]])
这相当于A:{4,3,3,1},C:{2,4,0,1},G:{1,0,5,2}和T:{1,1,0,4 }.检查它是否符合答案.
这是笨重的实施; 它在主体中不包含显式循环.除此之外,encode_sequence()可以通过一些离线准备输入来取代提前提高性能.虽然我没有上述片段的速度测量,但我认为它们应该加速到一定程度.:)
好吧,让我们讨论如果有长字符串会发生什么.
我们以此序列为例,
test_seq0 = ((
'A'*40, 'A'*40, 'A'*40, 'C'*40, 'C'*40,
'C'*40, 'C'*40, 'G'*40, 'G'*40, 'T'*40
))*4
将test_seq0包含40个字符串和每串有40个字符.看起来像这样,
In: len(test_seq0), test_seq0
Out: (40,
('AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA',
'AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA',
'AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA',
'CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC',
'CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC',
... skip 30 lines ...
'CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC',
'CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC',
'GGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGG',
'GGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGG',
'TTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTT'))
相当有趣的观点,对吗？
然后我们必须重新设置encode_sequence()一个长字符串版本,
def encode_sequence_longstring(tsq_np):
t00 = itertools.product('ACGT', repeat=4)
t01 = np.array(list(t00)).view('|U4').ravel()
t02 = dict(zip(t01, np.arange(len(t01))))
t03 = np.empty_like(tsq_np, dtype=np.uint)
t03.ravel()[:] = [t02[i] for i in tsq_np.ravel()]
return t03
请注意,tsq_np这里不再是一个简单的字符串列表.后缀_np意味着它现在是一个numpy数组.
然后test_seq0以一种笨拙的方式划分原件,
In: v01 = np.asarray(test_seq0).view('|U4').reshape(-1, int(40/4))
In: v01
Out:
array([['AAAA', 'AAAA', 'AAAA', 'AAAA', 'AAAA', 'AAAA', 'AAAA', 'AAAA', 'AAAA', 'AAAA'],
['AAAA', 'AAAA', 'AAAA', 'AAAA', 'AAAA', 'AAAA', 'AAAA', 'AAAA', 'AAAA', 'AAAA'],
['AAAA', 'AAAA', 'AAAA', 'AAAA', 'AAAA', 'AAAA', 'AAAA', 'AAAA', 'AAAA', 'AAAA'],
['CCCC', 'CCCC', 'CCCC', 'CCCC', 'CCCC', 'CCCC', 'CCCC', 'CCCC', 'CCCC', 'CCCC'],
['CCCC', 'CCCC', 'CCCC', 'CCCC', 'CCCC', 'CCCC', 'CCCC', 'CCCC', 'CCCC', 'CCCC'],
... skip 30 lines ...
['CCCC', 'CCCC', 'CCCC', 'CCCC', 'CCCC', 'CCCC', 'CCCC', 'CCCC', 'CCCC', 'CCCC'],
['CCCC', 'CCCC', 'CCCC', 'CCCC', 'CCCC', 'CCCC', 'CCCC', 'CCCC', 'CCCC', 'CCCC'],
['GGGG', 'GGGG', 'GGGG', 'GGGG', 'GGGG', 'GGGG', 'GGGG', 'GGGG', 'GGGG', 'GGGG'],
['GGGG', 'GGGG', 'GGGG', 'GGGG', 'GGGG', 'GGGG', 'GGGG', 'GGGG', 'GGGG', 'GGGG'],
['TTTT', 'TTTT', 'TTTT', 'TTTT', 'TTTT', 'TTTT', 'TTTT', 'TTTT', 'TTTT', 'TTTT']],
dtype='<U4')
另一个有趣的观点v01.:)
并用于v01计算这样的哈希码u02.它涉及围绕这些变量和函数的一些numpy约定.只是习惯那些花哨的技巧; 他们值得,
In: u02 = encode_sequence_longstring(v01)
In: u02
Out:
array([[ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
[ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
[ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
[ 85, 85, 85, 85, 85, 85, 85, 85, 85, 85],
[ 85, 85, 85, 85, 85, 85, 85, 85, 85, 85],
... skip 30 lines ...
[ 85, 85, 85, 85, 85, 85, 85, 85, 85, 85],
[ 85, 85, 85, 85, 85, 85, 85, 85, 85, 85],
[170, 170, 170, 170, 170, 170, 170, 170, 170, 170],
[170, 170, 170, 170, 170, 170, 170, 170, 170, 170],
[255, 255, 255, 255, 255, 255, 255, 255, 255, 255]],
dtype=uint64)
通过观察,你可以告诉u02它实际上是一对一的映射v01.它只是按预期将每个'AAAA'映射到0x0.
从现在开始,地图u02包含您需要的所有信息test_seq0.u02在numpy的帮助下将其提取出来,
s01 = np.empty((4, 0))
for u03 in u02.T:
p02, q02 = np.unique(u03, return_counts=True)
s02 = np.sum(map_decode[p02]*q02[:, None, None], axis=0).T
s01 = np.hstack((s01, s02))
python的本机循环有一个拇指规则:你可以使用它,但在任何对性能敏感的区域之外使用它.但是,这需要经验来判断手头的情况.
现在s01的预期报告如下,
In: s01
Out:
array([[ 12., 12., 12., 12., 12., 12., 12., 12., 12., 12., 12.,
12., 12., 12., 12., 12., 12., 12., 12., 12., 12., 12.,
12., 12., 12., 12., 12., 12., 12., 12., 12., 12., 12.,
12., 12., 12., 12., 12., 12., 12.],
[ 16., 16., 16., 16., 16., 16., 16., 16., 16., 16., 16.,
16., 16., 16., 16., 16., 16., 16., 16., 16., 16., 16.,
16., 16., 16., 16., 16., 16., 16., 16., 16., 16., 16.,
16., 16., 16., 16., 16., 16., 16.],
[ 8., 8., 8., 8., 8., 8., 8., 8., 8., 8., 8.,
8., 8., 8., 8., 8., 8., 8., 8., 8., 8., 8.,
8., 8., 8., 8., 8., 8., 8., 8., 8., 8., 8.,
8., 8., 8., 8., 8., 8., 8.],
[ 4., 4., 4., 4., 4., 4., 4., 4., 4., 4., 4.,
4., 4., 4., 4., 4., 4., 4., 4., 4., 4., 4.,
4., 4., 4., 4., 4., 4., 4., 4., 4., 4., 4.,
4., 4., 4., 4., 4., 4., 4.]])
从上到下阅读4行,它们分别是'A','C','G','T'.
与此同时,我尝试过test_seq0像这样的40x10x4000 ,
test_seq0 = ((
'A'*40, 'A'*40, 'A'*40, 'C'*40, 'C'*40,
'C'*40, 'C'*40, 'G'*40, 'G'*40, 'T'*40
))*4000
报告说,
array([[ 12000., 12000., 12000., 12000., 12000., 12000., 12000.,
12000., 12000., 12000., 12000., 12000., 12000., 12000.,
12000., 12000., 12000., 12000., 12000., 12000., 12000.,
12000., 12000., 12000., 12000., 12000., 12000., 12000.,
12000., 12000., 12000., 12000., 12000., 12000., 12000.,
12000., 12000., 12000., 12000., 12000.],
[ 16000., 16000., 16000., 16000., 16000., 16000., 16000.,
16000., 16000., 16000., 16000., 16000., 16000., 16000.,
16000., 16000., 16000., 16000., 16000., 16000., 16000.,
16000., 16000., 16000., 16000., 16000., 16000., 16000.,
16000., 16000., 16000., 16000., 16000., 16000., 16000.,
16000., 16000., 16000., 16000., 16000.],
[ 8000., 8000., 8000., 8000., 8000., 8000., 8000.,
8000., 8000., 8000., 8000., 8000., 8000., 8000.,
8000., 8000., 8000., 8000., 8000., 8000., 8000.,
8000., 8000., 8000., 8000., 8000., 8000., 8000.,
8000., 8000., 8000., 8000., 8000., 8000., 8000.,
8000., 8000., 8000., 8000., 8000.],
[ 4000., 4000., 4000., 4000., 4000., 4000., 4000.,
4000., 4000., 4000., 4000., 4000., 4000., 4000.,
4000., 4000., 4000., 4000., 4000., 4000., 4000.,
4000., 4000., 4000., 4000., 4000., 4000., 4000.,
4000., 4000., 4000., 4000., 4000., 4000., 4000.,
4000., 4000., 4000., 4000., 4000.]])
它在我的MacBook Pro上完成不到1秒钟(按ENTER键完成),这不是怪物.:)