查找允许某些不匹配的子字符串的快速方法

Question

我正在寻求帮助，以制定一种有效的方式来处理一些高通量 DNA 测序数据。数据在5个文件中，每个文件有几十万个序列，其中每个序列的格式如下：

@M01102:307:000000000-BCYH3:1:1102:19202:1786 1:N:0:TAGAGGCA+CTCTCTCT

TAATACGACTCACTATAGGGTTAACTTTAAGAGGGAGATATACATATGAGTCTTTTGGGTAAGAAGCCTTTTTGTCTGCTTTATGGTCCTATCTGCGGCAGGGCCAGCGGCAGCTAGGACGGGGGGCGGATAAGATCGGAAGAGCACTCGTCTGAACTCCAGTCACTAGAGGCAATCTCGT

+

AAABBFAABBBFGGGGFGGGGGAG5GHHHCH54BEEEEA5GGHDHHHH5BAE5DF5GGCEB33AF3313GHHHE255D55D55D53@5@B5DBD5@E/@//>/1??/?/E@///FDF0B?CC??CAAA;--./;/BBE?;AFFA./;/;.;AEA//BFFFF/BB/////;/..:.9999.;

我现在正在做的是迭代这些行，检查第一个和最后一个字母是否是 DNA 序列允许的字符（A/C/G/T 或 N），然后进行模糊搜索两个引物序列位于我感兴趣的编码序列片段的两侧。最后一步是出错的部分...

当我搜索完全匹配时，我会在合理的时间范围内获得可用数据。但是，我知道我错过了很多由于引物序列中的单个错配而被跳过的数据。发生这种情况是因为读取质量会随着长度的增加而降低，因此会出现更多不可读的碱基 ('N')。这些在我的分析中不是问题，而是简单的直接字符串搜索方法的问题——从 DNA 的角度来看，应该允许 N 与任何东西匹配，但从字符串搜索的角度来看则不是（我不太关心关于插入或删除）。出于这个原因，我正在尝试实施某种模糊或更符合生物学信息的搜索方法，但尚未找到一种有效的方法。

我现在拥有的东西确实适用于测试数据集，但速度太慢，无法用于全尺寸真实数据集。代码的相关片段是：

from Bio import pairwise2 
Sequence = 'NNNNNTAATACGACTCACTATAGGGTTAACTTTAAGAGGGAGATATACATATGAGTCTTTTGGGTAAGAAGCCTTTTTGTCTGCTTTATGGTCCTATCTGCGGCAGGGCCAGCGGCAGCTAGGACGGGGGGCGGATAAGATCGGAAGAGCACTCGTCTGAACTCCAGTCACTAGAGGCAATCTCGT'
fwdprimer = 'TAATACGACTCACTATAGGGTTAACTTTAAGAAGGAGATATACATATG'
revprimer = 'TAGGACGGGGGGCGGAAA'
if Sequence.endswith(('N','A','G','T','G')) and Sequence.startswith(('N','A','G','T','G')): 
    fwdalign = pairwise2.align.localxs(Sequence,fwdprimer,-1,-1, one_alignment_only=1)
    revalign = pairwise2.align.localxs(Sequence,revprimer,-1,-1, one_alignment_only=1)
    if fwdalign[0][2]>45 and revalign[0][2]>15:
        startIndex = fwdalign[0][3]+45
        endIndex = revalign[0][3]+3
        Sequence = Sequence[startIndex:endIndex]
        print Sequence

（显然在这个例子中不需要第一个条件，但有助于过滤掉其他 3/4 没有 DNA 序列的行，因此不需要搜索）

此方法使用 biopython 中的成对比对方法，该方法旨在查找允许错配的 DNA 序列比对。这部分做得很好，但是因为它需要用两个引物对每个序列进行序列比对，所以花费的时间太长而不实用。我需要它做的就是找到匹配序列，允许一两个不匹配。是否有另一种方法可以满足我的目标但在计算上更可行？为了进行比较，以前版本中的以下代码对我的完整数据集运行得非常快：

 if ('TAATACGACTCACTATAGGGTTAACTTTAAGAAGGAGATATACATATG' in Line) and ('TAGGACGGGGGGCGGAAA' in Line):
    startIndex = Line.find('TAATACGACTCACTATAGGGTTAACTTTAAGAAGGAGATATACATATG')+45
    endIndex = Line.find('TAGGACGGGGGGCGGAAA')+3
    Line = Line[startIndex:endIndex]
    print Line

这不是我运行经常做的事情，所以如果它有点低效请不要介意，但不想让它运行离开一整天。我想在几秒或几分钟内得到结果，而不是几小时。

Answer 1

tre library provides fast approximate matching函数。您可以使用 maxerr 指定最大不匹配字符数，如下例所示：

https://github.com/laurikari/tre/blob/master/python/example.py

还有regex module, which supports fuzzy searching options: https://pypi.org/project/regex/#additional-features

此外，您还可以使用简单的正则表达式来允许替代字符，如：

# Allow any character to be N
pattern = re.compile('[TN][AN][AN][TN]')

if pattern.match('TANN'):
    print('found')

查找允许某些不匹配的子字符串的快速方法

Fast way to find a substring with some mismatches allowed

python

bioinformatics

dna-sequence

biopython