这更多的是要找到最快的方法来做到这一点。我有一个 file1,它在单独的行中包含大约一百万个字符串(长度 6-40)。我想在另一个包含大约 80,000 个字符串的文件 2 中搜索它们中的每一个并计算出现次数(如果在一个字符串中多次找到小字符串,则该字符串的出现次数仍然为 1)。如果有人有兴趣比较性能,有下载 file1 和 file2 的链接。dropbox.com/sh/oj62918p83h8kus/sY2WejWmhu?m
我现在正在做的是为文件 2 构建一个字典,使用字符串 ID 作为键和字符串作为值。(因为 file2 中的字符串有重复值,只有字符串 ID 是唯一的)我的代码是
for line in file1:
substring=line[:-1].split("\t")
for ID in dictionary.keys():
bigstring=dictionary[ID]
IDlist=[]
if bigstring.find(substring)!=-1:
IDlist.append(ID)
output.write("%s\t%s\n" % (substring,str(len(IDlist))))
我的代码需要几个小时才能完成。任何人都可以建议一种更快的方法吗?file1 和 file2 都只有 50M 左右,我的电脑有 8G 内存,你可以使用尽可能多的内存来让它更快。任何可以在一小时内完成的方法都是可以接受的:)
在这里,在我尝试了下面这些评论的一些建议之后,请查看性能比较,首先是代码,然后是运行时。
Mark Amery 和其他人提出的一些改进import sys
from Bio import SeqIO
#first I load strings in file2 to a dictionary called var_seq,
var_seq={}
handle=SeqIO.parse(file2,'fasta')
for record in handle:
var_seq[record.id]=str(record.seq)
print len(var_seq) #Here print out 76827, which is the right number. loading file2 to var_seq doesn't take long, about 1 second, you shall not focus here to improve performance
output=open(outputfilename,'w')
icount=0
input1=open(file1,'r')
for line in input1:
icount+=1
row=line[:-1].split("\t")
ensp=row[0] #ensp is just peptides iD
peptide=row[1] #peptides is the substrings i want to search in file2
num=0
for ID,bigstring in var_seq.iteritems():
if peptide in bigstring:
num+=1
newline="%s\t%s\t%s\n" % (ensp,peptide,str(num))
output.write(newline)
if icount%1000==0:
break
input1.close()
handle.close()
output.close()
完成需要1m4s。与我的旧版本相比,提高了 20 多岁
#######熵建议的NEXT METHODfrom collections import defaultdict
var_seq=defaultdict(int)
handle=SeqIO.parse(file2,'fasta')
for record in handle:
var_seq[str(record.seq)]+=1
print len(var_seq) # here print out 59502, duplicates are removed, but occurances of duplicates are stored as value
handle.close()
output=open(outputfilename,'w')
icount=0
with open(file1) as fd:
for line in fd:
icount+=1
row=line[:-1].split("\t")
ensp=row[0]
peptide=row[1]
num=0
for varseq,num_occurrences in var_seq.items():
if peptide in varseq:
num+=num_occurrences
newline="%s\t%s\t%s\n" % (ensp,peptide,str(num))
output.write(newline)
if icount%1000==0:
break
output.close()
这需要 1 分 10 秒,因为它避免搜索重复项,所以没有像预期的那样快,不明白为什么。
Mark Amery 建议的 Haystack 和 Needle 方法,结果是最快的,这种方法的问题是所有子串的计数结果都是 0,我还不明白。这是我实现他的方法的代码。
class Node(object):
def __init__(self):
self.words = set()
self.links = {}
base = Node()
def search_haystack_tree(needle):
current_node = base
for char in needle:
try:
current_node = current_node.links[char]
except KeyError:
return 0
return len(current_node.words)
input1=open(file1,'r')
needles={}
for line in input1:
row=line[:-1].split("\t")
needles[row[1]]=row[0]
print len(needles)
handle=SeqIO.parse(file2,'fasta')
haystacks={}
for record in handle:
haystacks[record.id]=str(record.seq)
print len(haystacks)
for haystack_id, haystack in haystacks.iteritems(): #should be the same as enumerate(list)
for i in xrange(len(haystack)):
current_node = base
for char in haystack[i:]:
current_node = current_node.links.setdefault(char, Node())
current_node.words.add(haystack_id)
icount=0
output=open(outputfilename,'w')
for needle in needles:
icount+=1
count = search_haystack_tree(needle)
newline="%s\t%s\t%s\n" % (needles[needle],needle,str(count))
output.write(newline)
if icount%1000==0:
break
input1.close()
handle.close()
output.close()
只需0m11s即可完成,比其他方法快得多。但是,我不知道将所有计数结果都设为0是我的错误,还是Mark的方法存在缺陷。