搜索字符串,允许在字符串的任何位置发生一次不匹配

Question

我正在使用长度为25的DNA序列(参见下面的例子).我有一个230,000的清单,需要寻找整个基因组中的每个序列(弓形虫寄生虫).我不确定基因组有多大,但比230,000个序列长得多.

我需要查找每个25个字符的序列,例如,(AGCCTCCCATGATTGAACAGATCAT).

基因组被格式化为连续的字符串,即(CATGGGAGGCTTGCGGAGCCTGAGGGCGGAGCCTGAGGTGGGAGGCTTGCGGAGTGCGGAGCCTGAGCCTGAGGGCGGAGCCTGAGGTGGGAGGCTT ....)

我不关心它被发现的地点和次数,只关注它是否存在.
我认为这很简单 -

str.find(AGCCTCCCATGATTGAACAGATCAT)

但我还要找到在任何位置定义为错误(不匹配)的近距离匹配,但只有一个位置,并记录序列中的位置.我不知道怎么做到这一点.我唯一能想到的是使用通配符并在每个位置使用通配符执行搜索.即,搜索25次.

例如,
AGCCTCCCATGATTGAACAGATCAT
AGCCTCCCATGATAGAACAGATCAT

与位置13处的不匹配密切匹配.

速度不是一个大问题,因为我只做了3次,但如果它很快就会很好.

有些程序可以执行此操作 - 查找匹配项和部分匹配项 - 但我正在寻找一种使用这些应用程序无法发现的部分匹配项.

这是perl的类似帖子,虽然它们只是比较序列而不是搜索连续的字符串:

相关文章

Answer 1

在你继续阅读之前,你看过biopython吗？

您似乎希望找到具有一个替换错误的近似匹配,以及零插入/删除错误,即汉明距离为1.

如果你有汉明距离匹配功能(参见例如Ignacio提供的链接),你可以像这样使用它来搜索第一场比赛:

any(Hamming_distance(genome[x:x+25], sequence) == 1 for x in xrange(len(genome)))

但这会相当缓慢,因为(1)汉明距离函数会在失败后第二次替换错误(2)之后继续磨削,它会将光标前进一次,而不是根据它看到的东西(如Boyer-)向前跳过摩尔搜索确实).

你可以用这样的函数克服(1):

def Hamming_check_0_or_1(genome, posn, sequence):
    errors = 0
    for i in xrange(25):
        if genome[posn+i] != sequence[i]:
            errors += 1
            if errors >= 2:
                return errors
    return errors 

注意:这是故意不是Pythonic,它是Cic,因为你需要使用C(可能通过Cython)来获得合理的速度.

Navarro和Raffinot(google"Navarro Raffinot nrgrep")已经完成了一些关于平行近似Levenshtein搜索和跳过的工作,这可以适用于汉明搜索.请注意,位并行方法对查询字符串和字母大小的长度有限制,但是你的分别是25和4,所以没有问题.更新:跳过可能对字母大小为4的帮助不大.

当你谷歌进行汉明距离搜索时,你会注意到很多关于在硬件中实现它的东西,而不是软件.这是一个很大的暗示,也许你提出的任何算法应该用C或其他编译语言实现.

更新: 位并行方法的工作代码

我还提供了一种简单的方法来帮助进行正确性检查,并且我已经打包了一些Paul的重新编码以进行一些比较.请注意,使用re.finditer()会提供非重叠的结果,这会导致distance-1匹配阴影完全匹配; 看我上次的测试用例.

位并行方法具有以下特征:保证线性行为O(N),其中N是文本长度.注意天真方法是O(NM),因为正则表达式方法(M是模式长度).Boyer-Moore风格的方法将是最坏的情况O(NM)和预期的O(N).当输入必须被缓冲时,也可以很容易地使用位并行方法:它可以一次输入一个字节或一个兆字节; 没有预见,没有缓冲区边界的问题.最大的优势:用C编码时,简单的每输入字节代码的速度.

缺点:模式长度实际上限于快速寄存器中的位数,例如32或64.在这种情况下,模式长度为25; 没问题.它使用与模式中不同字符数成比例的额外内存(S_table).在这种情况下,"字母大小"仅为4; 没问题.

本技术报告的详细信息.该算法用于在Levenshtein距离内进行近似搜索.要转换为使用汉明距离,我只需(!)删除处理插入和删除的语句2.1.你会注意到很多带有"d"上标的"R"."d"是距离.我们只需要0和1.这些"R"成为下面代码中的R0和R1变量.

# coding: ascii

from collections import defaultdict
import re

_DEBUG = 0


# "Fast Text Searching with Errors" by Sun Wu and Udi Manber
# TR 91-11, Dept of Computer Science, University of Arizona, June 1991.
# http://citeseerx.ist.psu.edu/viewdoc/summary?doi=10.1.1.20.8854

def WM_approx_Ham1_search(pattern, text):
    """Generate (Hamming_dist, start_offset)
    for matches with distance 0 or 1"""
    m = len(pattern)
    S_table = defaultdict(int)
    for i, c in enumerate(pattern):
        S_table[c] |= 1 << i
    R0 = 0
    R1 = 0
    mask = 1 << (m - 1)
    for j, c in enumerate(text):
        S = S_table[c]
        shR0 = (R0 << 1) | 1
        R0 = shR0 & S
        R1 = ((R1 << 1) | 1) & S | shR0
        if _DEBUG:
            print "j= %2d msk=%s S=%s R0=%s R1=%s" \
                % tuple([j] + map(bitstr, [mask, S, R0, R1]))
        if R0 & mask: # exact match
            yield 0, j - m + 1
        elif R1 & mask: # match with one substitution
            yield 1, j - m + 1

if _DEBUG:

    def bitstr(num, mlen=8):
       wstr = ""
       for i in xrange(mlen):
          if num & 1:
             wstr = "1" + wstr
          else:
             wstr = "0" + wstr
          num >>= 1
       return wstr

def Ham_dist(s1, s2):
    """Calculate Hamming distance between 2 sequences."""
    assert len(s1) == len(s2)
    return sum(c1 != c2 for c1, c2 in zip(s1, s2))

def long_check(pattern, text):
    """Naively and understandably generate (Hamming_dist, start_offset)
    for matches with distance 0 or 1"""
    m = len(pattern)
    for i in xrange(len(text) - m + 1):
        d = Ham_dist(pattern, text[i:i+m])
        if d < 2:
            yield d, i

def Paul_McGuire_regex(pattern, text):
    searchSeqREStr = (
        '('
        + pattern
        + ')|('
        + ')|('.join(
            pattern[:i]
            + "[ACTGN]".replace(c,'')
            + pattern[i+1:]
            for i,c in enumerate(pattern)
            )
        + ')'
        )
    searchSeqRE = re.compile(searchSeqREStr)
    for match in searchSeqRE.finditer(text):
        locn = match.start()
        dist = int(bool(match.lastindex - 1))
        yield dist, locn


if __name__ == "__main__":

    genome1 = "TTTACGTAAACTAAACTGTAA"
    #         01234567890123456789012345
    #                   1         2

    tests = [
        (genome1, "ACGT ATGT ACTA ATCG TTTT ATTA TTTA"),
        ("T" * 10, "TTTT"),
        ("ACGTCGTAAAA", "TCGT"), # partial match can shadow an exact match
        ]

    nfailed = 0
    for genome, patterns in tests:
        print "genome:", genome
        for pattern in patterns.split():
            print pattern
            a1 = list(WM_approx_Ham1_search(pattern, genome))
            a2 = list(long_check(pattern, genome))
            a3 = list(Paul_McGuire_regex(pattern, genome))
            print a1
            print a2
            print a3
            print a1 == a2, a2 == a3
            nfailed += (a1 != a2 or a2 != a3)
    print "***", nfailed

Answer 2

Python 正则表达式库支持模糊正则表达式匹配.TRE的一个优点是它允许在文本中找到正则表达式的所有匹配(也支持重叠匹配).

import regex
m=regex.findall("AA", "CAG")
>>> []
m=regex.findall("(AA){e<=1}", "CAAG") # means allow up to 1 error
m
>>> ['CA', 'AG']

Answer 3

我搜索了"弓形虫寄生虫基因组",在网上找到了一些这些基因组文件.我发现了我认为很接近的文件,在http://toxodb.org上标题为"TgondiiGenomic_ToxoDB-6.0.fasta" ,大小约为158Mb.我使用以下pyparsing表达式提取基因序列,花了不到2分钟:

fname = "TgondiiGenomic_ToxoDB-6.0.fasta"
fastasrc = open(fname).read()   # yes! just read the whole dang 158Mb!

"""
Sample header:
>gb|scf_1104442823584 | organism=Toxoplasma_gondii_VEG | version=2008-07-23 | length=1448
"""
integer = Word(nums).setParseAction(lambda t:int(t[0]))
genebit = Group(">gb|" + Word(printables)("id") + SkipTo("length=") + 
                "length=" + integer("genelen") + LineEnd() + 
                Combine(OneOrMore(Word("ACGTN")),adjacent=False)("gene"))

# read gene data from .fasta file - takes just under a couple of minutes
genedata = OneOrMore(genebit).parseString(fastasrc)

(惊喜!一些基因序列包括'N'的运行!那到底是什么？!)

然后我把这个类写成了pyparsing Token类的子类,用于进行紧密匹配:

class CloseMatch(Token):
    def __init__(self, seq, maxMismatches=1):
        super(CloseMatch,self).__init__()
        self.name = seq
        self.sequence = seq
        self.maxMismatches = maxMismatches
        self.errmsg = "Expected " + self.sequence
        self.mayIndexError = False
        self.mayReturnEmpty = False

    def parseImpl( self, instring, loc, doActions=True ):
        start = loc
        instrlen = len(instring)
        maxloc = start + len(self.sequence)

        if maxloc <= instrlen:
            seq = self.sequence
            seqloc = 0
            mismatches = []
            throwException = False
            done = False
            while loc < maxloc and not done:
                if instring[loc] != seq[seqloc]:
                    mismatches.append(seqloc)
                    if len(mismatches) > self.maxMismatches:
                        throwException = True
                        done = True
                loc += 1
                seqloc += 1
        else:
            throwException = True

        if throwException:
            exc = self.myException
            exc.loc = loc
            exc.pstr = instring
            raise exc

        return loc, (instring[start:loc],mismatches)

对于每个匹配,这将返回包含匹配的实际字符串的元组和不匹配位置的列表.完全匹配当然会返回第二个值的空列表.(我喜欢这个课程,我想我会把它添加到下一个pyparsing版本中.)

然后我运行此代码来搜索从.fasta文件读取的所有序列中的"up-to-2-mismatch"匹配(回想一下,genedata是一个ParseResults组序列,每个组包含一个id,一个整数长度,以及一个序列字符串):

searchseq = CloseMatch("ATCATCGAATGGAATCTAATGGAAT", 2)
for g in genedata:
    print "%s (%d)" % (g.id, g.genelen)
    print "-"*24
    for t,startLoc,endLoc in searchseq.scanString(g.gene):
        matched, mismatches = t[0]
        print "MATCH:", searchseq.sequence
        print "FOUND:", matched
        if mismatches:
            print "      ", ''.join(' ' if i not in mismatches else '*' 
                            for i,c in enumerate(searchseq.sequence))
        else:
            print "<exact match>"
        print "at location", startLoc
        print
    print

我从一个基因位中随机取出了搜索序列,以确保我能找到一个完全匹配,只是出于好奇,看看有多少1和2元素不匹配.

这需要一段时间才能运行.45分钟后,我得到了这个输出,列出了每个id和基因长度,以及找到的任何部分匹配:

scf_1104442825154 (964)
------------------------

scf_1104442822828 (942)
------------------------

scf_1104442824510 (987)
------------------------

scf_1104442823180 (1065)
------------------------
...

我感到气馁,直到以下情况才看到任何比赛:

scf_1104442823952 (1188)
------------------------
MATCH: ATCATCGAATGGAATCTAATGGAAT
FOUND: ATCATCGAACGGAATCGAATGGAAT
                *      *        
at location 33

MATCH: ATCATCGAATGGAATCTAATGGAAT
FOUND: ATCATCGAATGGAATCGAATGGAAT
                       *        
at location 175

MATCH: ATCATCGAATGGAATCTAATGGAAT
FOUND: ATCATCGAATGGAATCGAATGGAAT
                       *        
at location 474

MATCH: ATCATCGAATGGAATCTAATGGAAT
FOUND: ATCATCGAATGGAATCGAATGGAAT
                       *        
at location 617

MATCH: ATCATCGAATGGAATCTAATGGAAT
FOUND: ATCATCGAATGGAATCGAATAGAAT
                       *   *    
at location 718

MATCH: ATCATCGAATGGAATCTAATGGAAT
FOUND: ATCATCGAATGGATTCGAATGGAAT
                    *  *        
at location 896

MATCH: ATCATCGAATGGAATCTAATGGAAT
FOUND: ATCATCGAATGGAATCGAATGGTAT
                       *     *  
at location 945

最后我完全匹配:

scf_1104442823584 (1448)
------------------------
MATCH: ATCATCGAATGGAATCTAATGGAAT
FOUND: ATCATCGAATGGACTCGAATGGAAT
                    *  *        
at location 177

MATCH: ATCATCGAATGGAATCTAATGGAAT
FOUND: ATCATCGAATGGAATCAAATGGAAT
                       *        
at location 203

MATCH: ATCATCGAATGGAATCTAATGGAAT
FOUND: ATCATCAAATGGAATCGAATGGAAT
             *         *        
at location 350

MATCH: ATCATCGAATGGAATCTAATGGAAT
FOUND: ATCATCGAATGGAATCGAATGGAAA
                       *       *
at location 523

MATCH: ATCATCGAATGGAATCTAATGGAAT
FOUND: ATCATCAAATGGAATCGAATGGAAT
             *         *        
at location 822

MATCH: ATCATCGAATGGAATCTAATGGAAT
FOUND: ATCATCGAATGGAATCTAATGGAAT
<exact match>
at location 848

MATCH: ATCATCGAATGGAATCTAATGGAAT
FOUND: ATCGTCGAATGGAGTCTAATGGAAT
          *         *           
at location 969

因此虽然没有设置任何速度记录,但我完成了工作,并且发现了一些2场比赛,以防它们可能有意义.

为了比较,这里是一个基于RE的版本,它只找到1个不匹配的匹配:

import re
seqStr = "ATCATCGAATGGAATCTAATGGAAT"
searchSeqREStr = seqStr + '|' + \
    '|'.join(seqStr[:i]+"[ACTGN]".replace(c,'') +seqStr[i+1:] 
             for i,c in enumerate(seqStr))

searchSeqRE = re.compile(searchSeqREStr)

for g in genedata:
    print "%s (%d)" % (g.id, g.genelen)
    print "-"*24
    for match in searchSeqRE.finditer(g.gene):
        print "MATCH:", seqStr
        print "FOUND:", match.group(0)
        print "at location", match.start()
        print
    print

(起初,我尝试搜索原始FASTA文件源本身,但很困惑为什么这么少的匹配与pyparsing版本相比.然后我意识到一些匹配必须跨越换行符,因为fasta文件输出包含在n字符).

因此,在第一次pyparsing传递以提取基因序列以匹配之后,这个基于RE的搜索者再花费大约1-1/2分钟来扫描所有未文本包装的序列,以找到所有相同的1-mismatch条目pyparsing解决方案做了.