寻找k聚体的滑动窗口

Question

我想通过这个问题，生物信息学工作：https://stepic.org/lesson/An-Explosion-of-Hidden-Messages-4/step/1?course=Bioinformatics-Algorithms-2&unit=8

的具体问题是在上面的链接的第5窗口，问题是：多少不同9- mers在大肠杆菌基因组中形成（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

有趣的问题。 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 

---

这里这个问题（这是在大数据共同）其实就是要选择合适的数据结构，并使得一些时间/空间权衡。如果您选择正确，您可以按照您的基因组长度和空间线性与滑动窗口的长度成线性关系。但我正在超越自己。让我们直观地解释问题（主要是这样，我可以理解它:-)）。

有3聚体在此窗口中的（20,3）-clump： “CAT”。还有其他一些（一个是“AAA”），但这个例子说明了k，L和t在做什么。

现在，在算法上。让我们进一步减少这个问题，以便我们可以想象我们将如何解析和存储它：让我们看看一个简单的（5,3）3个团体。

支架在这里表示我们的滑动宽度5的窗口。我们可以看到，在我们的窗口中，我们的3人分解为ATA,TAA和AAA。当我们向右滑动窗口时，ATA退出，我们获得第二个AAA。当我们将窗口向右滑动另一次时，现在TAA退出并且我们获得第三个AAA - 并且我们已经找到了一个（5,3）团块AAA s。

琐碎，显然，但对于弄清楚我们如何处理更大的团块很有用 - 重要的是，我们不会在我们移动窗口时丢弃整个先前窗口的数据;我们只丢弃第一个k-mer并在窗口的末尾添加新的k-mer。下一个见解是我们可以使用散列支持的结构（在python中，dict s）来做计数我们窗口内k-mers的。这不需要线性搜索我们的数据结构来找出有多少特定的k-mer在那里。

所以说起来这两个要求 - 插入和哈希支持的数据结构的记忆顺序 - 意味着我们应该维护一个list自定义类 - 或者更好，deque - 您在窗口各有k聚体，和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。随意拔叉等

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()')