作为我自己对基因型网络研究的一部分，我有一个代码块，我试图从字符串列表中构建一个单差网络。过程如下:
- 遍历字符串的所有成对组合。
- 如果字符串相差一个位置，则在它们之间在网络中画一条边。
- 如果它们没有一个位置的差异，那么通过。
我现在的代码块是这样的:
from itertools import combinations
import Levenshtein as lev # a package that wraps a C-implemented levenshtein distance
import networkx as nx
strings = [...list of strings...]
G = nx.Graph()
for (s1, s2) in combinations(strings, 2):
if s1 not in G.nodes():
G.add_node(s1)
if s2 not in G.nodes():
G.add_node(s2)
if lev.distance(s1, s2) == 1:
G.add_edge(s1, s2)
显然没有办法提高图构建过程的计算复杂度——它总是O(n**2)。至少，就我有限的知识而言，这就是我的想法——也许我错了?
也就是说，考虑到我需要进行的比较次数的正常规模(~大约)。2000-5000)，如果我能得到几个数量级的加速，总的来说，那么现实世界的计算时间仍然是可以接受的——使用当前的Python实现，构建一个图需要~天。有了正确的导入(下面没有说明)，我尝试了下面的一个Cython实现，但不知道如何使它更快:
cpdef cython_genotype_network(sequences):
G = nx.Graph()
cdef:
unicode s1
unicode s2
for (s1, s2) in combinations(sequences, 2):
if lev.distance(s1, s2) == 1:
G.add_edge(s1, s2)
return G
特别是，对于s1和s2, Cython期望bytes，而不是str。该代码块抛出一个错误。
所以…我来回答我的两个问题:
- Q1: Cython的实现会有帮助吗?我如何修复字节与str错误?Q2:是否有可能使用
numpy来解决这个问题?从numpy矩阵转换为NetworkX很容易;然而，我似乎无法弄清楚如何在n × n空矩阵中广播Levenshtein距离函数，其中每一行和列对应于字符串。UPDATE 1:如何生成样本数据
生成字符串:
from random import choice
def create_random_nucleotides_python(num_nucleotides):
"""
Creates random nucleotides of length num_nucleotides.
"""
sequence = ''
letters = ['A', 'T', 'G', 'C']
for i in range(num_nucleotides):
sequence = sequence + choice(letters)
return sequence
def mutate_random_position(string):
"""
Mutates one position in the nucleotide sequence at random.
"""
positions = [i for i in range(len(string))]
pos_to_mut = choice(positions)
letters = ['A', 'T', 'G', 'C']
new_string = ''
for i, letter in enumerate(string):
if i == pos_to_mut:
new_string = new_string + choice(letters)
else:
new_string = new_string + letter
return new_string
# Create 100 Python srings by mutating a first sequence, then randomly choosing stuff to mutate a single position.
base_sequence = create_random_nucleotides_python(1000)
sequences = [base_sequence]
for i in range(99):
sequence = choice(sequences)
mutseq = mutate_random_position(sequence)
sequences.append(mutseq)
关于复杂性:
你正在考虑每一对字符串。你不需要那个。您可以考虑每个字符串的所有1-距离字符串:
# I would start by populating the whole graph:
for s1 in strings:
if s1 not in G.nodes():
G.add_node(s1)
# Then compute the leven-1:
for s1 in strings:
for s2 in variations(s1):
if s2 in G.nodes():
G.add_edge(s1, s2)
现在你所需要的是一个比O(n)短的variations(string):
返回距离为1的所有变量。(only 1 edit|delete|insert)
def variations(string):
for i in range(len(string)):
# delete
yield string[:i] + string[i+1:]
# edit
for l in 'ATGC':
yield string[:i] + l + string[i+1:]
# insert
for l in 'ATGC':
yield string[:i] + l + string[i:]
# insert at the end
for l in 'ATGC':
yield string + l
现在，复杂度是O(m^2)(因为字符串concat)，其中m是最长序列的大小。如果它是已知的，它就是一个常量，现在所有的O(1)
如果序列都是相同的大小，你可以只计算编辑。
否则，您可以将序列从大到小排序，并且只计算编辑和删除。
或者，您可以按大小对字符串进行排序，而不是将所有字符串与其他字符串进行比较，而是比较那些大小相差为<= 1的字符串。
关于cython代码如何更快的一些想法
for (s1, s2) in combinations(sequences, 2):
if lev.distance(s1, s2) == 1:
G.add_edge(s1, s2)
我认为itertools.combinations是编译的，所以它本身是快速的。你调用它一次，所以for ...部分可能是最快的。尽管如此，在cython中自己循环sequences并不难。
lev.distance似乎使用了编译后的代码。是否有可能直接导入和调用该代码?查看lev源代码。
我认为lev.distance最终调用c函数签名:
distance_py(PyObject *self, PyObject *args)
或更多
lev_edit_distance(size_t len1, const lev_byte *string1,
size_t len2, const lev_byte *string2,
int xcost)
G.add_edge做什么?有没有一种更简单的数据结构可以收集这些边。也许是一个元组列表?
networkx是否提供了一种批量添加节点和边的方法?你总是必须使用add_node和add_edge吗?或者你可以给它一个节点列表，一个节点对列表(元组)?
看起来graph是一个Python字典，或者字典的字典。
networkx允许您添加一堆边，通过元组列表:
>>> G.add_edges_from([(1,2),(1,3)])
在python代码中，这部分可能是不必要的:
if s1 not in G.nodes():
G.add_node(s1)
if s2 not in G.nodes():
G.add_node(s2)
只做:
G.add_nodes_from(strings)
节点http://networkx.github.io/documentation/latest/tutorial/tutorial.html