您的示例输出表明每个子序列只需要一个索引,但您的描述暗示您可能需要所有索引。
这是一个函数,它将构建一个字典,其中包含存在的子序列(给定长度)的所有索引:
from collections import defaultdict
def sequenceDict(seq,count):
result = defaultdict(list)
for i,subseq in enumerate(zip(*(seq[i:] for i in range(count)))):
result["".join(subseq)].append(i)
return result
r = sequenceDict('GACTCGTAGCTAGCTG',3)
print(r)
# {'GAC': [0], 'ACT': [1], 'CTC': [2], 'TCG': [3], 'CGT': [4], 'GTA': [5], 'TAG': [6, 10], 'AGC': [7, 11], 'GCT': [8, 12], 'CTA': [9], 'CTG': [13]})
如果你真的只想要每个 3 字母子序列的第一个索引,那么使用单个字典理解可以更快地获得字典:
from itertools import product
{ ss:sequence.index(ss) for p in product(*["ACGT"]*3)for ss in ["".join(p)] if ss in sequence}
我对 2.5 亿个字母的随机序列进行了性能测试,可以在几微秒内获得单索引字典。获取所有索引需要一分钟多一点(使用上面的函数):
import time
size = 250_000_000
print("loading sequence...",size)
start = time.time()
import random
sequence = "".join(random.choice("ACGT") for _ in range(size))
print("sequence ready",time.time()-start)
start = time.time()
from itertools import product
seqDict = { ss:sequence.index(ss) for p in product(*["ACGT"]*3)for ss in ["".join(p)] if ss in sequence}
print("\n1st index",time.time()-start)
start = time.time()
r = sequenceDict(sequence,3)
print("\nall indexes",time.time()-start)
输出:
loading sequence... 250000000
sequence ready 193.82172107696533
1st index 0.000141143798828125
all indexes 71.74848103523254
鉴于加载序列的时间比构建索引的时间长得多,您可能会放弃存储该索引字典并每次都从源数据中重建它(您似乎仍然需要为您的进程加载)
您也可以只存储一个计数字典并根据需要提取索引:
此函数获取每个子序列的出现次数:
from collections import Counter
def countSubSeqs(seq,size):
return Counter("".join(s) for s in zip(*(seq[i:] for i in range(size))))
它的运行时间与 sequenceDict 函数大致相同,但生成的字典要小得多。
要获取特定子序列(包括重叠位置)的索引,您可以使用:
subSeq = "ACT"
indexes = [ i for i in range(len(sequence)) if sequence[i:i+3]==subSeq ]
如果您不需要立即为所有子序列创建所有索引,则可以相应地构建代码并仅在需要时获取索引(并可能将它们存储在字典中以供查询和重用)