使用 python 查找并删除公共子字符串

find and remove common substring with python

所以我有两个引物序列:

五素数 = "GATTCGAAGTCCACTATC"
三素数 = "TGAGTAGGACGGCACTATC"

我需要做的是当我有另一个序列时我需要检查它是否包含一个部分 这些引物序列中的一个,如果需要,我需要删除匹配部分,留下不匹配部分以供进一步分析。

例如我的序列:CACTATCAAAAAAA 里面有 fiveprimer 的一部分。我需要找到这个公共子字符串,然后将其删除,只留下 AAAAA

我 运行 遇到的问题是当你有这样的序列时 CACTATCGAAG 其中 GAAG 也在引物中,但不是引物序列的一部分它仍然被删除。我试图通过确保找到的结构位于引物的左侧并在右侧使用三引物来解决此问题,例如:

对于 CACTATCGAAG,我们有 2 个通用结构 CACTATCGAAG 现在我可以将 CACTATC 与 fiveprime GATTCGAAGTCCACTATC 的末尾进行比较,并在匹配时告诉它是引物的一部分,然后将其删除。因此,当我们将 GAAG 与五引物的最后一端的长度进行比较时,它会给我们这个 GATTCGAAGTCCACTATC ,它不匹配所以 GAAG可以继续进一步分析。

出于某种原因,我的脚本出现错误或无法正常工作。对于这个问题还有其他的解决方法或者建议吗?

def longestSubstringFinder(string1, string2):
answer = ""
len1, len2 = len(string1), len(string2)
for i in range(len1):
    match = ""
    for j in range(len2):
        if i + j < len1 and string1[i + j] == string2[j]:
            match += string2[j]
        else:
            if len(match) > len(answer): answer = match
            match = ""
return answer

def get_sequence(test, fiveprime, threeprime):
    if test == fiveprime:
        pass
    elif test == threeprime:
        pass
    elif test in fiveprime:
        pass
    elif test in threeprime:
        pass

        # find out if there is a matching part between the primers and the found
        # single stranded region, then calculates what part that is, and removes it
        # from the single stranded region
    else:
        overlap = longestSubstringFinder(fiveprime, test)
        l = len(overlap)

        if fiveprime[-l:] == test[:l]:
            check = test[l:]
        else:
            check = test

        overlap2 = longestSubstringFinder(check, threeprime)
        l = len(overlap2)

        if threeprime[:l] == check[-l:]:
            check2 = check[:-l]
            structure.append(check2)
        else:
            structure.append(check)

return structure

我认为如果您 select 适当的数据结构来表示您要查找的数据,您的问题会更容易处理。我能想到的最好的是 trie.

这种结构的好处是,它允许您表示给定初始字母序列的所有可能匹配项,因此如果您有序列 AABAB,它将允许从初始 A 遍历到 A 和 B,但不是从 A 到 G 或 T 的遍历。这使得查找部分匹配变得高效,因为 trie 中的任何点都代表那么多字母的匹配。

这个数据结构类似于:

class Trie(object):
    def __init__(self):
        self.children = {}

    def add_child(self, letter):
        if letter in self.children:
            return self.children[letter]
        else:
            child = Trie()
            self.children[letter] = child
            return child

    def traverse(self, letter):
        return self.children.get(letter, None)

然后您可以像这样填充它:

root = Trie()
current_positions = []
for letter in letters:
    current_positions = [
        position.add_child(letter)
        for position in current_positions
    ]
    current_positions.append(root.add_child(letter))

一旦设置了所有这些,您就应该能够遍历此结构,直到遍历 returns null。这将指示最长的当前匹配。字母的初始化将每个字母视为匹配的潜在起点,您也应该如此。

然后您可以像这样搜索最长的子串匹配:

class TrieSearch(object):
    def __init__(self, trie, starting_index):
        self.trie = trie
        self.starting_index = starting_index
        self.ending_index = starting_index + 1

    def update(self, letter):
        """ This returns a boolean indicating
            if the search can accept the letter """
        self.trie = self.trie.traverse(letter)
        if self.trie is not None:
            self.ending_index = self.ending_index + 1
            return True
        return False

    def get_match(self, letters):
        return letters[self.starting_index:self.ending_index]

def find_matches(root, letters):
    completed_matches = []
    current_matches = []

    for index, letter in enumerate(letters):
        new_current = []

        for current in current_matches:
            if current.update(letter):
                new_current.append(current)
            else:
                completed_matches.append(current)

        new_search_trie = root.traverse(letter)
        if new_search_trie is not None:
            new_current.append(TrieSearch(new_search_trie, index))

        current_matches = new_current

    all_matches = completed_matches + current_matches
    return [match.get_match(letters) for match in all_matches]

我已将所有这些放在 a gist 中,当使用 threeprimefiveprime 值初始化 trie 并且输入数据为 CACTATCAAAAAAA 结果时是:

['CACTATC', 'ACTATC', 'CTATC', 'TATC', 'ATC', 'TC', 'CA', 'AA', 'AA', 'AA', 'AA', 'AA', 'AA', 'A']

由于您无疑要处理大量字符串,因此您可能需要复习一些更有效的通用字符串替换算法。 Aho-Corasick algorithm is an extension of the trie approach outlined here. There is also the Knuth-Morris-Pratt algorithm 使用表而不是 trie。这两种算法都是线性复杂度算法,将大大改进您的 longestSubstringFinder 方法使用的二次方法。

问题是你的 longestSubstringFinder 坏了。如果条件 if len(match) > len(answer): answer = matchFalse,则您只检查 if len(match) > len(answer): answer = match 。如果 for j in range(len2): 循环结束而没有发生这种情况,您只需丢弃 match.

修复很简单:在 for j in range(len2): 循环之后添加 for j in range(len2):

def longestSubstringFinder(string1, string2):
    answer = ""
    len1, len2 = len(string1), len(string2)
    for i in range(len1):
        match = ""
        for j in range(len2):
            if i + j < len1 and string1[i + j] == string2[j]:
                match += string2[j]
            else:
                if len(match) > len(answer): answer = match
                match = ""
        if len(match) > len(answer): answer = match # this was missing
    return answer