06-算法02-激动人心的新线索

北野茶缸子

发布于 2022-05-19 11:47:55

2890

发布于 2022-05-19 11:47:55

文章被收录于专栏：北野茶缸子的专栏

参考：
- 生物信息学 | Coursera[1]

1-正去反来都一样

DNA 都是反向互补的：

因此AGTCGCATAGT 与ACTATGCGACT 是互补序列。为什么是ACTATGCGACT呢？因为DNA 的复制规则是单向的，其只允许从5' 到3' 复制。所以其互补序列可能是出现在基因组的反链上。

p：其实我自己也有点不太理解为何要计算互补序列。虽然互补序列在复制之后会得到原有链一样的片段（方向相反），但这个假设前提是二者存在于oriC 即复制起点的两侧，对应在forward 与reverse 的区域。如果AGTCGCATAGT 与ACTATGCGACT 两个片段就在彼此的附近，显然复制的结果是相反的。而且有意思的是，课程后面计算clump 时，并没有强调要计算互补，为了严谨考虑，我后面的代码只有部分考虑了互补（最终并未考虑），以及设计了对应的输出。毕竟互补计算只是多了一个判断，多了一个转换，并不影响效率。

给定一个序列，我们可以获得其互补序列吗？

这里非常简单：

def reverse_complement_1(pattern):
    complement = []
    for i in pattern:
        if i == 'A':
            complement.append('T')
        elif i == 'T':
            complement.append('A')
        elif i == 'C':
            complement.append('G')
        elif i == 'G':
            complement.append('C')
        else :
            print('warning')
            break
    complement.reverse()
    re_complement = "".join(complement)
    return re_complement

对每个碱基进行判断，并取其互补碱基，最后将得到的列表取反即可。

更为简单的操作是使用字典，根据输入片段作为键，去查找其互补的值，最后再对字符串取反即可：

def reverse_complement_2(pattern):
    re_pattern = ''
    dict = {'A':'T', 'T':'A', 'C':'G', 'G':'C'}
    for i in pattern:
        re_pattern += dict[i]
    return re_pattern[::-1]

字母大小写直接强制转换会好一点吗？

def reverse_complement_1(pattern):
    complement = []
    for i in pattern:
        if i == 'a':
            complement.append('t')
        elif i == 't':
            complement.append('a')
        elif i == 'c':
            complement.append('g')
        elif i == 'g':
            complement.append('c')
        else :
            print('warning')
            break
    complement.reverse()
    re_complement = "".join(complement)
    return re_complement

def reverse_complement_2(pattern):
    re_pattern = ''
    dict = {'a':'t', 't':'a', 'c':'g', 'g':'c'}
    for i in pattern:
        re_pattern += dict[i]
    return re_pattern[::-1]

现在我们发现，对于Vibrio cholerae 的复制起始区域来说：

atcaatgatcaacgtaagcttctaagcatgatcaaggtgctcacacagtttatccacaacctgagtggatgacatcaagataggtcgttgtatctccttcctctcgtactctcatgaccacggaaagatgatcaagagaggatgatttcttggccatatcgcaatgaatacttgtgacttgtgcttccaattgacatcttcagcgccatattgcgctggccaaggtgacggagcgggattacgaaagcatgatcatggctgttgttctgtttatcttgttttgactgagacttgttaggatagacggtttttcatcactgactagccaaagccttactctgcctgacatcgaccgtaaattgataatgaatttacatgcttccgcgacgatttacctcttgatcatcgatccgattgaagatcttcaattgttaattctcttgcctcgactcatagccatgatgagctcttgatcatgtt tccttaaccctctattttttacggaagaatgatcaagctgctgctcttgatcatcgtttc

其输出的9-mers 中：

$ python3 03-k_mer_pattern_freq_array.py Input/genome02.txt 9
['atgatcaag', 'ctcttgatc', 'tcttgatca', 'cttgatcat']

atgatcaag 和cttgatcat 是互补的。这绝对不是巧合。

为什么会出现这种情况呢？我们先按下不表。不过如果加上这两个互补序列，总共的计数就是6了，自然atgatcaag/cttgatcat 是起始段中重复最多的序列了。

接下来，我们还可以尝试获取这些片段在基因上的位置。

import sys

sys.path.append('../2-exciting_new_clues')
from step04_reverse_complement import reverse_complement_2 

text_path = sys.argv[1]

with open(text_path) as a:
  text = a.read()
    
pattern = sys.argv[2]

def patternLocalizing(Text, Pattern):
 rev_pattern = reverse_complement_2(Pattern)
 count = []
 for i in range(len(Text) - len(Pattern) - 1):
  if Pattern == Text[i:i+len(Pattern)] or rev_pattern == Text[i:i+len(Pattern)]:
   count.append(i)
 return count

def main():
 print(patternLocalizing(text, pattern))

if __name__ == "__main__":
 main()

这个非常简单，滑动整个序列，看是否有正反片段与其相等，返回位置信息即可。

% python3 step05_pattern_localization.py ./Input/genome03.txt atgatcaag
[27, 129, 403, 475, 516, 533]

2-成功揭秘真不错

ps：事先说明，这一部分的代码这里存在一些问题。

至此为止，我们实现了以下功能：

正反链的互换；
获得指定长度序列下的k-mer；
获得指定片段在序列中的位置；

但别忘了，这里我们始终都是对某个长度的序列进行计算的。开始的定义：某个长度L区域内，k长度的碱基重复t次，并且t 是该长度下重复出现最大的次数，为k-mer (L,t)-clump。这个k-mer 则是重复出现长度为k 的片段，即在k 长度下，该片段重复次数最多。

如果我们想要在一段更长的序列，比如基因组上，得到 k-mer (L,t)-clump 呢？

比如下面的序列：

gatcagcataagggtccCTGCAATGCATGACAAGCCTGCAGTtgttttac

TGCA 就是 4-mer(25,3) 的clump。

如果我们给一个更长的序列三个参数4-mer(25,3)，可以得到对应输出吗？并输出它的位置信息？

这里的解法也很简单：

对整个大的序列按照L 进行滑动；
每次滑动内，计算滑动结果下的frequency_array，字典记录frequency_pattern 及其count。整个代码思路可以直接使用

这里我并不想指定t，t 直接为输出的最大值即可。当然代码修改也很简单。

import sys

text_path = sys.argv[1]

with open(text_path) as a:
  text = a.read()
    
k = int(sys.argv[2])
L = int(sys.argv[3])


def reverse_complement_2(pattern):
    re_pattern = ''
    dict = {'a':'t', 't':'a', 'c':'g', 'g':'c'}
    for i in pattern:
        re_pattern += dict[i]
    return re_pattern[::-1]


def fasterFrequentWords(Text, k, frequency_array):
  for i in range(len(Text) - k + 1):
    pattern = Text[i:i+k]
    rev_pattern = reverse_complement_2(pattern)
    if pattern in frequency_array:
      frequency_array[pattern] += 1
    elif rev_pattern in frequency_array:
      frequency_array[rev_pattern] += 1
    else:
      frequency_array[pattern] = 1
  return frequency_array


def patternLocalizing(Text, Pattern):
  rev_pattern = reverse_complement_2(Pattern)
  count = []
  for i in range(len(Text) - len(Pattern) - 1):
    if Pattern == Text[i:i+len(Pattern)] or rev_pattern == Text[i:i+len(Pattern)]:
      count.append(i)
  return count


def clumpFinding(genome, k, L):
  frequency_array = {}
  frequent_patterns_location = {}
  for i in range(len(genome) - L):
    gene = genome[i:i+L]
    frequency_array = fasterFrequentWords(gene, k, frequency_array)
  max_count = max(frequency_array.values())
  frequent_patterns = [k for k,v in frequency_array.items() if v==max_count]
  for i in frequent_patterns:
    location = patternLocalizing(genome, i)
    frequent_patterns_location[i] = location
  return(frequent_patterns_location)


def main():
  print(clumpFinding(text, k, L))

if __name__ == "__main__":
  main()

% python3 step06_find_clumps.py ./Input/genome03.txt 9 50
{'atgatcaag': [27, 127, 397, 468, 508, 525]}

尝试一下更大的 E. coli genome，2h 都没跑出来。

不难发现，这里的运行速度慢了许多，毕竟是一个长度四百万的基因组。

这里的clumpFinding 算法，还有更大的优化空间吗？

不过这里也存在一个小小的bug，比如在k 较小时，比如：

atcaatgatcaacgtaag

tgatc 和gatca 这两个就是互相挨着的序列，该如何避免这种情况呢？

再次强调，这一部分的代码这里存在一些问题。

3-正确理解clump的概念

上面的代码，我犯了非常严重的概念问题，我没有理解clump 的意义。

某个长度L区域内，k长度的碱基重复t次，并且t 是该长度下重复出现最大的次数，为k-mer (L,t)-clump。这个k-mer 则是重复出现长度为k 的片段，即在k 长度下，该片段重复次数最多。

既然是指定了某个长度L区域内，那如果跳到了下个长度L 的区域，之前array 里计算的那些序列就都不算了。因此需要重新计数，而不是累加：

  for i in range(len(genome) - L):
    gene = genome[i:i+L]
    frequency_array = fasterFrequentWords(gene, k, frequency_array)

因此，正确的代码如下，每一次滑动，获得滑动组别中的最大值及pattern，最后再将全部的最大组进行比较，获得最大组的最大组。我还加入了一些计算时间的功能：

import sys
import time

text_path = sys.argv[1]

with open(text_path) as a:
  text = a.read()
    
k = int(sys.argv[2])
L = int(sys.argv[3])


def reverse_complement_2(pattern):
    re_pattern = ''
    dict = {'a':'t', 't':'a', 'c':'g', 'g':'c'}
    for i in pattern:
        re_pattern += dict[i]
    return re_pattern[::-1]


def fasterFrequentWords(Text, k):
  frequency_array = {}
  for i in range(len(Text) - k + 1):
    pattern = Text[i:i+k]
    rev_pattern = reverse_complement_2(pattern)
    if pattern in frequency_array:
      frequency_array[pattern] += 1
    elif rev_pattern in frequency_array:
      frequency_array[rev_pattern] += 1
    else:
      frequency_array[pattern] = 1
  max_count = max(frequency_array.values())
  max_frequency_array_tmp = {k:v for k,v in frequency_array.items() if v==max_count}
  return max_frequency_array_tmp


def patternLocalizing(Text, Pattern):
  rev_pattern = reverse_complement_2(Pattern)
  count = []
  for i in range(len(Text) - len(Pattern) - 1):
    if Pattern == Text[i:i+len(Pattern)] or rev_pattern == Text[i:i+len(Pattern)]:
      count.append(i)
  return count


def clumpFinding(genome, k, L):
  genome = genome.lower()
  frequent_patterns_location = {}
  max_frequency_array = {}
  for i in range(len(genome) - L):
    gene = genome[i:i+L]
    max_frequency_array_tmp = fasterFrequentWords(gene, k)
    max_frequency_array.update(max_frequency_array_tmp)
  max_count = max(max_frequency_array.values())
  frequent_patterns = [k for k,v in max_frequency_array.items() if v==max_count]
  for i in frequent_patterns:
    location = patternLocalizing(genome, i)
    frequent_patterns_location[i] = {"location" : location, 
                                     "rev_pattern" : reverse_complement_2(i)}
  frequent_patterns_location["counts"] = max_count
  return(frequent_patterns_location)


def main():
  time_start=time.time()
  print(clumpFinding(text, k, L))
  time_end=time.time()
  print('totally cost %f seconds' % (time_end-time_start))

if __name__ == "__main__":
  main()

输出：

 % python3 step06_find_clumps.py ./Input/genome04.txt 9 500
{'gtgcgggct': {'location': [244, 278, 293], 'rev_pattern': 'agcccgcac'}, 'counts': 3}
totally cost 27.875074 seconds

% python3 step06_find_clumps.py ./Input/genome03.txt 8 100
{'cttgatca': {'location': [28, 128, 397, 468, 509, 525], 'rev_pattern': 'tgatcaag'}, 'ttgatcat': {'location': [4, 27, 127, 398, 469, 508, 526], 'rev_pattern': 'atgatcaa'}, 'counts': 3}
totally cost 0.069806 seconds

很快，我又发现了一处错误：

  for i in range(len(genome) - L):
    gene = genome[i:i+L]
    max_frequency_array_tmp = fasterFrequentWords(gene, k)
    max_frequency_array.update(max_frequency_array_tmp)

update 语法形如dict1.update(dict2)，因为字典特性规定了键需要唯一，因此，如果dict2 中存在和dict1 中一样的键，则前者中的值会被后者“更新”。也就是说，如果开始clump 找到了某个max_count，其可能会被后面clump 中相对较小的max_count 替代。

因此这里需要更正，其实也就是增加一个判断，如果dict2 中的maxcount 更大，则dict2 替换dict1，反之亦然。反正我们关心的也是出现频数最多的k-mer。

  for i in range(len(genome) - L + 1):
    gene = genome[i:i+L]
    max_frequency_array_tmp = fasterFrequentWords(gene, k)
    if i > 0:
      if max(max_frequency_array_tmp.values()) > max(max_frequency_array.values()):
        max_frequency_array.update(max_frequency_array_tmp)
      else:
        max_frequency_array_tmp.update(max_frequency_array)
        max_frequency_array = max_frequency_array_tmp
    else:
      max_frequency_array.update(max_frequency_array_tmp)

但这里同样存在问题，如果刚好是：

>>> a = {"a":3, "b":2}
>>> b = {"a":3, "b":3}

虽然我们关心的是maxcount k-mer，但并不代表只有一个maxcount k-mer。

好家伙，难道我就要一个个循环比较吗？

  for i in range(len(genome) - L + 1):
    gene = genome[i:i+L]
    max_frequency_array_tmp = fasterFrequentWords(gene, k)
    if i > 0:
      tmp_pattern = max_frequency_array_tmp.keys()
      array_pattern = max_frequency_array.keys()
      intersect_pattern = [x for x in tmp_pattern if x in array_pattern]
      diff_pattern = tmp_pattern - array_pattern
      diff_array = {k:v for k,v in max_frequency_array_tmp.items() if k in diff_pattern}
      for x in intersect_pattern:
        if max_frequency_array_tmp[x] > max_frequency_array[x]:
          max_frequency_array[x] = max_frequency_array_tmp[x] 
      max_frequency_array.update(diff_array)
    else:
      max_frequency_array.update(max_frequency_array_tmp)

交集则循环比较，差集则添加到之前的array 中：

import sys
import time

text_path = sys.argv[1]

with open(text_path) as a:
  text = a.read()
    
k = int(sys.argv[2])
L = int(sys.argv[3])


def reverse_complement_2(pattern):
    re_pattern = ''
    dict = {'a':'t', 't':'a', 'c':'g', 'g':'c'}
    for i in pattern:
        re_pattern += dict[i]
    return re_pattern[::-1]


def fasterFrequentWords(Text, k):
  frequency_array = {}
  for i in range(len(Text) - k + 1):
    pattern = Text[i:i+k]
    # rev_pattern = reverse_complement_2(pattern)
    if pattern in frequency_array:
      frequency_array[pattern] += 1
    # elif rev_pattern in frequency_array:
    #   frequency_array[rev_pattern] += 1
    else:
      frequency_array[pattern] = 1
  max_count = max(frequency_array.values())
  max_frequency_array_tmp = {k:v for k,v in frequency_array.items() if v==max_count}
  return max_frequency_array_tmp



def patternLocalizing(Text, Pattern):
  # rev_pattern = reverse_complement_2(Pattern)
  count = []
  for i in range(len(Text) - len(Pattern) - 1):
    if Pattern == Text[i:i+len(Pattern)] :
    # if Pattern == Text[i:i+len(Pattern)] or rev_pattern == Text[i:i+len(Pattern)]:
      count.append(i)
  return count


def clumpFinding(genome, k, L):
  genome = genome.lower()
  frequent_patterns_location = {}
  max_frequency_array = {}
  for i in range(len(genome) - L + 1):
    gene = genome[i:i+L]
    max_frequency_array_tmp = fasterFrequentWords(gene, k)
    if i > 0:
      tmp_pattern = max_frequency_array_tmp.keys()
      array_pattern = max_frequency_array.keys()
      intersect_pattern = [x for x in tmp_pattern if x in array_pattern]
      diff_pattern = tmp_pattern - array_pattern
      diff_array = {k:v for k,v in max_frequency_array_tmp.items() if k in diff_pattern}
      for x in intersect_pattern:
        if max_frequency_array_tmp[x] > max_frequency_array[x]:
          max_frequency_array[x] = max_frequency_array_tmp[x] 
      max_frequency_array.update(diff_array)
    else:
      max_frequency_array.update(max_frequency_array_tmp)
  max_count = max(max_frequency_array.values())
  frequent_patterns = [k for k,v in max_frequency_array.items() if v==max_count]
  # for i in frequent_patterns:
  #   location = patternLocalizing(genome, i)
  #   frequent_patterns_location[i] = {"location" : location, 
  #                                    "rev_pattern" : reverse_complement_2(i)}
  # frequent_patterns_location["counts"] = max_count
  # return(frequent_patterns_location)
  frequent_patterns2 = {max_count : frequent_patterns, "times" : len(frequent_patterns)}
  return(frequent_patterns2)


def main():
  time_start=time.time()
  result = clumpFinding(text, k, L)
  print(result)
  time_end=time.time()
  print('totally cost %f seconds' % (time_end-time_start))

if __name__ == "__main__":
  main()

% python3 step06_find_clumps.py ./Input/genome04.txt 8 500                               
{4: ['ggcggcgg'], 'times': 1}
totally cost 32.157041 seconds

上面的代码里，我取消了互补片段的判断，即认为互补片段也是一个独立的片段。对应课程中的问题。

再来看看更新array 部分代码：

  for i in range(len(genome) - L + 1):
    gene = genome[i:i+L]
    max_frequency_array_tmp = fasterFrequentWords(gene, k)
    if i > 0:
      tmp_pattern = max_frequency_array_tmp.keys()
      array_pattern = max_frequency_array.keys()
      intersect_pattern = [x for x in tmp_pattern if x in array_pattern]
      diff_pattern = tmp_pattern - array_pattern
      diff_array = {k:v for k,v in max_frequency_array_tmp.items() if k in diff_pattern}
      for x in intersect_pattern:
        if max_frequency_array_tmp[x] > max_frequency_array[x]:
          max_frequency_array[x] = max_frequency_array_tmp[x] 
      max_frequency_array.update(diff_array)
    else:
      max_frequency_array.update(max_frequency_array_tmp)

不难发现，我们每次得到的tmp_pattern，其大部分的内容都是没有改变的。但我们还是重复计算了其k-mer 及计数，获得tmp_array，并用这个tmp_array 重复比较之前的array。

你能想到代码优化的方案吗？

对了，发现公众号的代码格式是默认不换行的，比如我上面的存放碱基的代码块，只能通过滑动浏览全部，感觉很不好看，你们有什么解决方案吗？

参考资料

[1]

生物信息学 | Coursera: https://www.coursera.org/specializations/bioinformatics#courses

本文参与腾讯云自媒体同步曝光计划，分享自微信公众号。

原始发表：2022-04-19，如有侵权请联系 cloudcommunity@tencent.com 删除

数据分析

本文分享自北野茶缸子微信公众号，前往查看

如有侵权，请联系 cloudcommunity@tencent.com 删除。

本文参与腾讯云自媒体同步曝光计划，欢迎热爱写作的你一起参与！

数据分析

登录后参与评论

0 条评论

热度