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

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序列的行，因此不需要搜索）

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

^{pr2}$

这不是我经常运行的东西，所以不介意它是否有点低效，但不想让它运行一整天。我想在几秒钟或几分钟内得到结果，而不是几个小时。在

Answer 1

tre library提供快速的approximate matching函数。{cd1>可以指定最大字符数为

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

还有一个regex module，它支持模糊搜索选项：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')