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

Question

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

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

TAATACGACTCACTATAGGGTTAACTTTAAGAGGGGATATACATATGAGTCTTTTGGGTAAGAAGCCTTTTTGTCTGCTTTATGGTCCTATCTGCGGCAGGGCCAGCGGCAGCTAGGACGGGGGGCGGATAAGATCGGAAGAGCACTCGTCTGAACTCCAGTCACTAGAGGCAATCTCGT

+

AAABBFAABBBFGGGGFGGGGGAG5GHHHCH54BEEEEA5GHDHHHH5BAE5DF5GGCEB33AF3313GHHHE255D55D55D53@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 提供快速的approximate matching 函数。您可以使用maxerr 指定最大不匹配字符数，如下例所示：

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')