如何有效地遍历一个大字符串并记录每个3个字母的子字符串的索引？

Question

我有一个非常大的序列，我以字符串的形式读入（250000000个字母）。字母是G，A，C，T

例如：

'GACTCGTAGCTAGCTG'

我想创建一些方法来存储序列中每个3个字母子字符串的索引，以便稍后在另一个函数中使用

例如：

{'GAC': 1, 'ACT': 2 'CTC':3, 'TCG':4 ...}

对于我目前的方法，我的问题是我没有找到一种有效的方法来存储序列中每个3个字母子字符串的索引。一旦我知道每个子串的索引，我将根据给定的概率随机选择其中的一些，并将它们更改为另一个已知的子串

我尝试过使用for循环迭代，在滑动窗口中一次使用一个3个字母的子字符串，并将位置保存为字典值，使用3个字母的子字符串作为键。另外，当我保存这个字典文件时，我一直在使用pd.to_csv，但它似乎效率很低，并且创建了一个8 GB的文件。但是，对于一个非常大的字符串（250000000个字母），这需要太长的时间

Answer 1

可以将长字符串拆分为3个字母的子字符串：

string='GACTCGTAGCTAGCT'
substrings=[string[3*x:3*x+3] for x in range(int(len(string)/3))]

substrings将是：

['GAC', 'TCG', 'TAG', 'CTA', 'GCT']

将这些标记添加到另一个列表：

indices=[x for x in range(int(len(string)/3))]

这只会产生：

[0, 1, 2, 3, 4]

子字符串列表中的第n个元素将对应于索引列表中的第n个元素

有关如何将文件放入字符串变量，请参见： How to read a text file into a string variable and strip newlines?

Answer 2

您的示例输出表明，每个子序列只需要一个索引，但您的描述暗示您可能需要所有索引

下面是一个函数，该函数将为存在的子序列（给定长度）构建一个包含所有索引的字典：

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 ] 

如果您不需要立即为所有子序列创建所有索引，那么您可以相应地构造代码，只在需要时获取索引（并可能将它们存储在字典中以便查询和重用）