在滑动窗口中找到k-mers

Question

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

具体问题在上面链接的第5个窗口,问题是: 大肠杆菌基因组中有多少个不同的9聚体形成(500,3)个团块？(换句话说,不要多次计算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)

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-mers的团块:"CAT".还有一些其他的(一个"AAA"),但这个例子说明了k,L和t正在做什么.

现在,关于算法.让我们进一步减少问题所以我们可以想象我们如何解析和存储它:让我们看看一个简单的(5,3) - 3-mers的团队.

括号在这里表示我们的宽度为5的滑动窗口.我们可以看到,在我们的窗口我们的3聚物分解到ATA,TAA和AAA.当我们向右滑动窗口时,ATA退出并获得一秒钟AAA.当我们再向右滑动窗口时,现在TAA退出并获得第三个AAA- 我们找到了一个(5,3)AAAs的团块.

很明显,这对于弄清楚我们如何处理更大的团块很有用 - 重要的是,当我们移动窗口时,我们不会丢弃整个先前窗口的数据 ; 我们只是丢弃第一个k-mer并在窗口末尾添加新的k-mer.下一个见解是我们可以使用哈希支持的结构(在python,dicts中)来计算窗口内的k-mers.这消除了线性搜索我们的数据结构以确定其中有多少特定k聚体的需要.

所以说起来这两个要求-插入的记忆顺序和哈希支持的数据结构-意味着我们应该做的是保持一个自定义类list-或更好,deque-您在窗口各有k聚体,和dict-或更好,Counter - 跟踪你的双端队列中每个kmer的频率.请注意,OrderedDict接近为您完成所有这些,但并不完全; 如果计数大于1,弹出最老的kmer将是错误的.

你真正应该用来简化你的代码的另一件事是一个合适的滑动窗口迭代器.

把它们放在一起:

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 这里.随意叉,等