在滑动窗口中寻找k-mer

Question

我正在尝试解决这个生物信息学问题：https://stepic.org/lesson/An-Explosion-of-Hidden-Messages-4/step/1?course=Bioinformatics-Algorithms-2&unit=8 具体问题在上面链接的第5个窗口中，问题是：在 E. coli 基因组中有多少个不同的 9-mers 形成了（500，3）clumps？（换句话说，不要重复计算一个 9-mer。） 以下是我的代码。它是错误的，我很想知道为什么以及如何改进它（显然 O 效率很差，但我几天前才开始编写 Python...）非常感谢！

genome = '' #insert e. Coli genome here
k = 4 #length of k-mer
L = 50 #size of sliding window
t = 3 #k-mer appears t times
counter = 0
Count = []


for i in range(0,len(genome)-L): #slide window down the genome
    pattern = genome[i:i+k] #given this k-mer
    for j in range(i,i+L): #calculate k-mer frequency in window of len(L)
        if genome[j:j+k] == pattern:
            counter = counter + 1
    Count.append(counter)
    counter = 0 #IMPORTANT: reset counter after each i

Clump = []
for i in range(0,len(Count)):
    if Count[i] == t: #figure out the window that has k-mers of frequency t
        Clump.append(i)

Output = []
for i in range(0,len(Clump)):
    Output.append(genome[Clump[i]:Clump[i]+k])
print " ".join(list(set(Output))) #remove duplicates if a particular k-mer is found more than once
print len(Output)
print len(list(set(Output))) #total number of Clump(k,L,t)

- epicode

---

错误403：未订阅课程的用户无法访问链接。 - Hugues Fontenelle

什么是（500,3）-clump？ - Hugues Fontenelle

你的代码有什么问题？出现了错误信息吗？（然后复制它）或者是输出结果有误吗？（那么请同时复制实际输出和期望输出） - Hugues Fontenelle

我看到了很多的索引、计数器和for循环。Python不是Matlab或C。看一下Python教程吧！ - Hugues Fontenelle

抱歉，这里有一个问题的解释：https://www.dropbox.com/s/qcb8mrc7fab2ra5/Screenshot%202014-10-29%2012.51.11.png?dl=0 - epicode

一个(500,3)-clump指的是给定长度为500的滑动窗口和k-mer频率为3，你能找到多少个长度为k的唯一k-mer？我的代码有效，对于给定k=9，L=500和t=3的大肠杆菌基因组，它输出“872”（大肠杆菌基因组在此处：https://stepic.org/media/attachments/lessons/4/E-coli.txt）根据Google的说法，正确答案是1904。是的 - 我正在努力使它更具Python风格！ - epicode

Answer 1

有趣的问题。我在这里放了一个带有一些测试的实现，可以在github上查看。继续阅读以获取更多解释。

ben@nixbox:~/bin$ time python kmers.py ../E-coli.txt 9 500 3
(500, 3)-clumps of 9-mers found in that file: 1904

real    0m15.510s
user    0m14.241s
sys     0m0.956s

---

这个问题（在大数据中很常见）实际上归结于选择正确的数据结构，并进行一些时间/空间权衡。如果你选择正确，你可以在线性时间内完成，与你的滑动窗口长度成线性关系的空间。但是我有点超前了。让我们通过可视化方式来解释这个问题（主要是为了让我理解它:-)）。

在这个窗口中存在一个(20,3)的3-mer族群："CAT"。还有一些其他的（例如"AAA"），但这个例子说明了k、L和t是如何起作用的。

现在，让我们进入算法。为了能够可视化地解析和存储它，让我们进一步简化问题：看看一个简单的(5,3)的3-mer族群。

括号表示我们这里的滑动窗口宽度为5。我们可以看到，在我们的窗口中，我们的3-mers分解为ATA，TAA和AAA。当我们将窗口向右滑动一个时，ATA消失了，我们获得了第二个AAA。当我们再向右滑动窗口时，现在TAA消失了，我们获得了第三个AAA - 我们找到了一个(5,3)的AAA块。

显然，这很简单，但对于确定如何处理更大的块非常有用 - 重要的是，当我们移动窗口时，我们不会丢弃整个先前窗口的数据；我们只是丢弃第一个k-mer，并将新的k-mer附加到我们窗口的末尾。下一个见解是，我们可以使用哈希支持的结构（在Python中，dict）来进行我们窗口内k-mer的计数。这消除了通过线性搜索我们的数据结构来确定其中特定k-mer数量的需要。

因此，这两个要求 - 记住插入的顺序和基于哈希的数据结构 - 意味着我们应该制作一个自定义类来维护您窗口中每个kmer的deque列表（或更好）和Counter字典，以跟踪您的deque中每个kmer的频率。请注意，OrderedDict接近为您完成所有这些工作，但并非完全如此；如果最旧的kmer计数大于1，则弹出它是错误的。

另一件真正应该在此处用于简化代码的事情是适当的滑动窗口迭代器。

将它们全部组合起来：

def get_clumps(genome, k, L, t):
    kmers = KmerSequence(L-k, t)

    for kmer in sliding_window(genome, k):
        kmers.add(kmer)

    return kmers.clumps

class KmerSequence(object):
    __slots__ = ['order', 'counts', 'limit', 'clumps', 't']

    def __init__(self, limit, threshold):
        self.order = deque()
        self.counts = Counter()
        self.limit = limit
        self.clumps = set()
        self.t = threshold

    def add(self, kmer):
        if len(self.order) > self.limit:
            self._remove_oldest()
        self._add_one(kmer)

    def _add_one(self,kmer):
        self.order.append(kmer)
        new_count = self.counts[kmer] + 1
        self.counts[kmer] = new_count

        if new_count == self.t:
            self.clumps.add(kmer)

    def _remove_oldest(self):
        self.counts[self.order.popleft()] -= 1

使用方法：

with open(genomefile) as f:
    genome = f.read()

k = 9
L = 500
t = 3

clumps = get_clumps(genome, k,L,t)

正如顶部所提到的，完整的代码 - 包括一些辅助函数和处理脚本直接作为__main__运行时的情况 - 可在github 这里找到。请随意fork等操作。

- roippi

Answer 2

又一个实现

def findClumps(genome, k, L, t):
    length = len(genome) - k + 1
    indexes = defaultdict(list)
    result = set()

    for i in range(length):
        kterm = genome[i:i + k]

        # remove positions with distance > L
        while indexes[kterm] and i + k - indexes[kterm][0] > L:
            indexes[kterm].pop(0)

        # add current kterm position
        indexes[kterm].append(i)

        if len(indexes[kterm]) == t:
            result.add(kterm)

    return result

filename = 'E-coli.txt'
file = open(filename)
genome = file.read()
k=9
L=500
t=3

def test():
    print findClumps(genome, k, L, t)

import cProfile

cProfile.run('test()')

- Arnaldo P. Figueira Figueira