在滑动窗口中查找 k-mers

Posted

技术标签:

【中文标题】在滑动窗口中查找 k-mers【英文标题】:Finding k-mers in a sliding window 【发布时间】:2014-12-24 15:30:52 【问题描述】:

我正在尝试解决这个生物信息学问题:https://stepic.org/lesson/An-Explosion-of-Hidden-Messages-4/step/1?course=Bioinformatics-Algorithms-2&unit=8

具体问题在上面链接的第5个窗口,问题是:大肠杆菌基因组中有多少个不同的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)

【问题讨论】:

错误 403:未订阅课程的人无法使用问题链接。 什么是 (500,3)-团块? 你的代码有什么问题?错误信息? (然后复制它)或错误的输出? (然后复制它,以及预期的输出) 我看到很多索引、计数器和 for 循环。 Python 不是 Matlab 或 C。请查看 Python 教程! 对不起,这里是问题的解释:dropbox.com/s/qcb8mrc7fab2ra5/… 【参考方案1】:

有趣的问题。 I've put up an implementation with a few tests on github here。继续阅读以获得一些解释。

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)-clump 的 3-mers:“CAT”。还有一些其他的(“AAA”代表一个),但这个例子说明了 k、L 和 t 的作用。

现在,进入算法。让我们进一步简化这个问题,以便我们可以可视化我们如何解析和存储它:让我们看一个简单的 3-mers (5,3)-clump。

括号表示我们这里宽度为 5 的滑动窗口。我们可以看到,在我们的窗口中,我们的 3-mer 分解为 ATATAAAAA。当我们向右滑动一个窗口时,ATA 退出,我们获得第二个AAA。当我们再次将窗口向右滑动时,现在TAA 退出,我们获得了第三个AAA - 我们发现了AAAs 的(5,3) 块。

显然是微不足道的,但对于弄清楚我们如何处理更大的团块很有用 - 重要的是,当我们移动窗口时,我们不会丢弃整个前一个窗口的数据;我们只是丢弃第一个 k-mer 并将新的 k-mer 附加到我们窗口的末尾。下一个见解是,我们可以使用哈希支持结构(在 python 中,dicts)在我们的窗口内对 k-mers 进行计数。这消除了对我们的数据结构进行线性搜索以找出其中有多少特定 k-mer 的需要。

所以这两个要求——记住插入顺序和哈希支持的数据结构——意味着我们应该创建一个自定义类来维护你拥有的每个 kmer 的 list - 或者更好,deque -你的窗口,还有一个dict - 或者更好的Counter - 来跟踪你的双端队列中每个kmer的频率。请注意,OrderedDict 几乎可以为您完成所有这些工作,但并不完全;如果其计数大于 1,则弹出最旧的 kmer 将是错误的。

真正应该用来简化代码的另一件事是正确的sliding window iterator。

把它们放在一起:

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 here。随意分叉等。

【讨论】:

【参考方案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()')

【讨论】:

以上是关于在滑动窗口中查找 k-mers的主要内容,如果未能解决你的问题,请参考以下文章

力扣

1030 完美数列 (25 分)(滑动窗口)

“智能”滑动滑动窗口? [关闭]

滑动窗口9:239. 滑动窗口最大值和480滑动窗口中位数

LeetCode-滑动窗口滑动窗口最大值

算法总结之滑动窗口