主页 / 问题 / 860336
Accepted
Tony Stark
Asked:2020-07-26 04:40:49 +0000 UTC2020-07-26 04:40:49 +0000 UTC

查找字符串中给定长度 k 的最频繁子串

  • 772

在 coursera 上学习了生物信息学课程。有一个我无法解决的问题。关键是我们有两个变量：

  1. 包含重复字母组合的长单词中的字母列表（例如 ACGTTGCATGTCGCATGATGCATGAGAGCT 和 CATG GCAT）
  2. 一个数字，指定重复字母的数量（在上面的示例中，给出了一个数字

  3. 我们需要找出哪些组合是最重复的。我找到了一个解决方案，但无法理解一行：

    Count = CountDict(Text, k)

    我抛出一个错误，据我所知，它的含义是创建字典。

    # Input:  A string Text and an integer k
# Output: A list containing all most frequent k-mers in Text
def FrequentWords(Text, k):
    FrequentPatterns = []
    Count = CountDict(Text, k)
    m = max(Count.values())
    for i in Count:
        if Count[i] == m:
            FrequentPatterns.append(Text[i:i+k]
    FrequentPatternsNoDuplicates = remove_duplicates(FrequentPatterns)
    return FrequentPatternsNoDuplicates
python

5 个回答

  1. Best Answer

    我认为生物信息学意味着一个有效的解决方案，因为。基因组字符串可能很长。

    而一种有效的解决方案，尤其是对于大k，在这种情况下是使用一个后缀数组结合确定两个相邻后缀（LCP）的最大公共前缀。通过简单的方式构建后缀数组具有复杂性O(NlogN)或O(Nlog^2(N)). 还有一些算法更难在线性时间内实现。

    编辑：
    在代码中注释掉了一个带有错误渐近的简短 Python 实现。从这里
    替换为 Manber-Myers 算法的实现（渐近线应该是 O (NlogN)，但是，据我了解，该实现不为任意输入数据提供此功能）

    LCP 是在线性时间中构建的——我给出了一个高效且简单的 Kasai 算法。

    构建完LCP后，需要从中分离出最长的序列，所有的元素都不小于给定的长度k（我认为这是在Python中的一行完成的）。时间是线性的。

    最大序列的长度对应于（更准确地说，减一）所需长度的最频繁子串的重复次数。要获取此子字符串本身，您需要从后缀数组的相应元素的索引开始获取一个子字符串。

    #медленная реализация
#def get_suffix_array(s):
#    return sorted(range(len(s)), key=lambda i: s[i:])

from collections import defaultdict

def sort_bucket(s, bucket, order):
    d = defaultdict(list)
    for i in bucket:
        key = s[i + order // 2:i + order]
        d[key].append(i)
    result = []
    for k, v in sorted(d.items()):
        if len(v) > 1:
            result += sort_bucket(s, v, 2 * order)
        else:
            result.append(v[0])
    return result


def suffix_array_ManberMyers(s):
    return sort_bucket(s, range(len(s)), 1)

def lcp_kasai(s, suffarr):
    n = len(suffarr)
    k = 0
    lcp = [0] * n
    rank = [0] * n
    for i in range(n):
        rank[suffarr[i]] = i

    for  i in range(n):
        if (k>0):
            k -= 1
        if(rank[i]==n-1):
             k = 0
             continue
        j = sa[rank[i]+1]
        while((i+k<n) & (j+k<n) & (s[i+k]==s[j+k])):
            k += 1
        lcp[rank[i]] = k
    return lcp

sa = suffix_array_ManberMyers("ACGTTGCATGTCGCATGATGCATGAGAGCT$")
print(sa)
lc = lcp_kasai("ACGTTGCATGTCGCATGATGCATGAGAGCT$", sa)
print(lc)

    后缀数组输出和 lcp：

    [30, 0, 24, 26, 21, 14, 17, 7, 20, 13, 6, 11, 1, 28, 23, 25, 16, 19, 12, 5, 
 27, 9, 2, 29, 10, 22, 15, 18, 4, 8, 3]
[0, 1, 2, 1, 4, 3, 3, 0, 5, 4, 1, 2, 1, 0, 3, 2, 1, 6, 5, 2, 1, 2, 0, 1, 1,
 3, 2, 6, 2, 1, 0]

    在 lcp 中，我们看到两个长度为 2 的块（意味着两个子字符串出现 3 次），其值 >=4：5, 4 и 6, 5从位置 8 和 17 开始。后缀数组中的这些位置包含原始字符串 20 和 19 中的索引，它们对应于子字符串CATG和GCAT

    • 8

  2. 一个简单直接的算法：计算每个长度为 [重叠] 子串的出现次数k，然后输出重复次数最多的子串：

    #!/usr/bin/env python
"""Find most common substring of length k."""
from collections import Counter

text = "ACGTTGCATGTCGCATGATGCATGAGAGCT"
k = 4

# count occurrences of all k-mers: words of length k in the text
words = Counter(text[i:i+k] for i in range(len(text) - k + 1))  # O(n*k)
print(words.most_common(1)[0][0])  #-> GCAT


    <script src="https://cdn.rawgit.com/brython-dev/brython/3.4.0/www/src/brython.js"></script><script src="https://cdn.rawgit.com/brython-dev/brython/3.4.0/www/src/brython_stdlib.js"></script><body onload="brython()"><script type="text/python">
from collections import Counter
from browser import document

@document["mybutton"].bind("click")
def on_click(event):
    text = document['text'].value
    k = int(document['k'].value)
    words = Counter(text[i:i+k] for i in range(len(text) - k + 1))
    print(words.most_common()[0][0])

on_click('dummy on start')
</script><label for="k">k&nbsp;=</label>&nbsp;<input id="k" value="4"> <label for="text">text&nbsp;=</label>&nbsp;<input id="text" value="ACGTTGCATGTCGCATGATGCATGAGAGCT"> <button id="mybutton">Найти k-мер</button></body>

    • 4 
  3. from collections import defaultdict, Counter

text = 'ACGTTGCATGTCGCATGATGCATGAGAGCT'
key_len = 4

accumulator = Counter(text[i: i + key_len] for i in range(len(text) - key_len + 1))

print(accumulator)  # {'A': 7, 'C': 6, 'G': 9, 'T': 7,

max_items = defaultdict(list)

for k, v in accumulator.items():
    max_items[v].append(k)

print(max_items)  # ... , 2: ['TGCA', 'ATGA'], 3: ['GCAT', 'CATG']})

# Находим ключ с максимальным значением
max_key = max(max_items)
print(max_key, max_items[max_key])  # 3 ['GCAT', 'CATG']
    • 2

  4. 尝试类似：

    def CountDict(Text, k):
    # Создаём полный список всех подстрок длиной k
    if k > len(Text): return None
    if k == len(Text): return {Text:1}

    list_substr = list()   # Список подстрок длиной k
    for j in range(len(Text) - k + 1):
        list_substr.append( Text[j:j+k])

    # Список превращаем в словарь, с подсчётом частот
    RetMap = dict()
    for j in range(len(list_substr)):
        key_val = list_substr[j]
        RetMap[key_val] = RetMap.get(key_val, 0) + 1

    return RetMap

dnk_map = CountDict('ACGTTGCATGTCGCATGATGCATGAGAGCT', 4)
print(dnk_map)
    • 0

  5. 另一种选择（我发现很难估计算法的复杂性）：

    from itertools import product
text = "ACGTTGCATGTCGCATGATGCATGAGAGCT"

cnt = [[i,text.count(''.join(i))] for i in product('ACGT',repeat=4)]
m = max(cnt,key=lambda x: x[1])[1]
patterns = [''.join(i) for i,_ in filter(lambda x: x[1]==m,cnt)]

print(patterns) # ['CATG', 'GCAT']
    • -1

相关问题