bioinfo08-算法04-复制起点你在哪？

北野茶缸子

发布于 2022-05-19 11:26:54

6790

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

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

1-像聚合酶一样思考

1.1-做个成熟的聚合酶

DNA 都是反向互补的：

当DNA 的双螺旋结构打开后，其会形成两股复制叉，且这两股复制叉会朝着各自的方向沿着染色体（两条单链）进行复制：

但这种想法过于天真，DNA 聚合酶只能够从5' 到3' 复制，因此我们只能够往一个方向复制了QAQ，这样的话一条链只有一半可以复制了，所以如果你是DNA 聚合酶，你会怎么去复制呢？

DNA 的复制被发现只能由5' 到3'，而上面的复制规则却是双向的。为了实现另一条链复制，反向的DNA复制只能以一段一段（冈崎片段）的形式出现，因为DNA解旋的方向与复制方向相反。

冈崎片段的发现者：

DNA 复制过程有以下注意事项：

DNA 的复制是边解链边复制的；
在每条链的复制起点（oriC），会结合两个DNA 聚合酶，沿着相反的方向复制，各自完成整条单链的一半复制工作；
DNA 的复制是半不连续的；

而在复制的结果中，正向复制的5'-3' 片段称为reverse，因为其相比起被复制的片段来说，它是反向的；而另一段则为forward：

完整过程如下：

forward与reverse二者的区别，除了在于其复制特征不同，它们的速度也有所差别。由于forward 需要不断将其冈崎片段进行拼接，以形成完整的链，其速度相比reverse 要慢得多。因此大部分情况下，forward 都是处于双链状态的（double-stranded），而reverse 则是单链状态（single-stranded）。这种单双链的差别非常重要，因为单链DNA 相比双链DNA，有高得多的突变率。

而因为正反链碱基比例上存在的差异，也就为我们寻找复制起点提供了一个新方法。

1.2-寻找起点的新方法

比如观察下表：

不难发现，反链和正链碱基计数值C、G 分别为+11617、-9973。而互补配对原则告诉我们，C=G。

其中一种解释是，正链中的C 由于脱氨基反应（deamination），突变成了T，从而产生了GT错配，因此导致正链中C 的计数，低于G 的计数。另外，GT 错配还可能进一步在DNA 修补步骤变成AT，因此导致反链中互补的G 减少。

反链中的G 相对较少，C 相对较多，我们可以利用这个特点进行判断：如果我们持续计算GC 差值，如果差值增大，就说明我们可能在正链上，反之亦然。

制定一个规则，遍历序列，遇到C 减一，遇到G 加一：

如果遍历整个基因组：

通过这个图的峰和谷我们就可以得到复制的终点和起点。如果此时能够确定链的正反，比如是正链，则方向是从5' 到3'，如上图所示，正链的G-C数值首先逐渐增大，到达复制终点后，转为负链，数值不断减小。因此，这个峰图中的最低点，即是我们关心的复制起点的位置！

比如一段输入：

TAAAGACTGCCGAGAGGCCAACACGAGTGCTAGAACGAGGGGCGTAAACGCGGGTCCGAT

我们通过代码实现一下。

核心思想有两个：

将字符串转成列表，再将这个列表通过字典转成ATCG 对应的数值，比如C 为-1，AT 都为0；
将这个列表计算累加，每个索引下的数值是其与之前全部数值之和。

import sys

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

def pattern2Count(Text):
  count_list = list()
  text_list = list(Text.lower())
  dict = {'a':0, 't':0, 'c':-1, 'g':1}
  count_list = [dict[i] for i in text_list]
  return count_list

def findMinLocation(Text):
  cumcount_list = list()
  count_list = pattern2Count(Text)
  a = 0
  for i in count_list:
    a += i
    cumcount_list.append(a)
  # max_count = max(cumcount_list)
  min_count = min(cumcount_list)
  # max_location = [x for x,y in enumerate(cumcount_list) if y == max_count
  min_location = [x+1 for x,y in enumerate(cumcount_list) if y == min_count]  
  result = {"min":{"count":min_count, "location":min_location}}
  return result


def main():
  result = findMinLocation(text)
  print(result)

if __name__ == "__main__":
  main()

输出：

% python3 4-find_OriC/step09_find_GC_max.py 4-find_OriC/Input/genome2.txt
{'min': {'count': -411, 'location': [61733, 61734, 61735, 61737, 61739]}}

习惯ggplot 绘图的我，也来尝试一下matplotlib 的基本绘图：

# 首先将cumcount_list 结果转换为 pandas 对象
pd_df = pd.DataFrame(enumerate(cumcount_list))
>>> pd_df.head()
   0  1
0  0 -1
1  1  0
2  2  1
3  3  1
4  4  0

import matplotlib.pyplot as plt

# Create a Figure and an Axes with plt.subplots
fig, ax = plt.subplots()
ax.plot(pd_df[1])
# Call the show function
plt.show()

这个图应该是课程考核制作的随机生成的序列，所以不太正常也是理所当然的啦。

2-隐藏线索真难琢磨

2.1-小小错配大奥妙

通过上面的方法，我们发现E.coli 基因组在 3923620 位置有最小的GC 差值，正如上面介绍的那样，这是一个潜在的复制起点。

接下来，我们在这个复制起点的周围看看能否发现DnaA，即9-mer：

aatgatgatgacgtcaaaaggatccggataaaacatggtgattgcctcgcataacgcggt atgaaaatggattgaagcccgggccgtggattctactcaactttgtcggcttgagaaaga cctgggatcctgggtattaaaaagaagatctatttatttagagatctgttctattgtgat ctcttattaggatcgcactgccctgtggataacaaggatccggcttttaagatcaacaac ctggaaaggatcattaactgtgaatgatcggtgatcctggaccgtataagctgggatcag aatgaggggttatacacaactcaaaaactgaacaacagttgttctttggataactaccgg ttgatccaagcttcctgacagagttatccacagtagatcgcacgatctgtatacttattt gagtaaattaacccacgatcccagccattcttctgccggatcttccggaatgtcgtgatc aagaatgttgatcttcagtg

令人大失所望，竟然没有任何的9-mer 出现次数大于或等于三次。

让我们回到最开始已知的那个Vibrio cholerae 的oriC：

我们发现有两个片段ATGATCAAC 与 CATGATCAT，其和我们关注的k-mer，仅仅有一个碱基的差异。

这种情况我们称之为，DNA 碱基的错配（mismatch）。

如果有两段我们认为应该是相同的k-mer，p,q，对应碱基位置有p1 … pk 及 q1 … qk，但如果pi ≠ qi，则说明i 位置发生了错配。这个错配的计数，称为Hamming distance。

比如：

GGGCCGTTGGT
GGACCGTTGAC

有三个碱基错配，则Hamming distance 为3 。

尝试撸个代码吧。我的想法也很简单，循环比一下，不对就加一：

import sys

text_path1 = sys.argv[1]
text_path2 = sys.argv[2]
with open(text_path1) as a:
  text1 = a.read()
with open(text_path2) as a:
  text2 = a.read()


def calHammingDistance(Text1, Text2):
  a = 0
  if len(Text1) != len(Text2):
    return "Length of Text1 must be equal with Text2!"
  for i in range(len(Text1)):
    if Text1[i] != Text2[i]:
      a += 1
  return a


def main():
  result = calHammingDistance(text1, text2)
  print(result)

if __name__ == "__main__":
  main()

输出：

% python3 4-find_OriC/step10_cal_Hamming_distance.py ./4-find_OriC/Input/text1.txt 4-find_OriC/Input/text2.txt  
854

接下来就是新的问题了，我们能否在一段长度L 的序列中，找到Hamming distance 在某个容忍数值大小的全部k-mer 呢？

比如输入k-mer，L 以及某个Hamming distance：

ATTCTGGA
CGCCCGAATCCAGAACGCATTCCCATATTTCGGGACCACTGGCCTCCACGGTACGGACGTCAATCAAAT
3

输出全部的小于容忍数值大小的对应k-mer 的位置，比如输入3，则输出的对应位置的pattern 和比较的pattern 的Hamming distance 都是小于3 的。

import sys

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


def calHammingDistance(Text1, Text2):
  a = 0
  if len(Text1) != len(Text2):
    return "Length of Text1 must be equal with Text2!"
  for i in range(len(Text1)):
    if Text1[i] != Text2[i]:
      a += 1
  return a

def findHammingKmerLocation(Text, Pattern, d):
  location = str()
  for i in range(len(Text) - len(Pattern) + 1):
    pattern = Text[i:i+len(Pattern)]
    hd_score = calHammingDistance(pattern, Pattern)
    if hd_score <= d:
      location += str(i);
      location += " "
  return location

def main():
  result = findHammingKmerLocation(text, pattern, d)
  print(result)

if __name__ == "__main__":
  main()

% python3 4-find_OriC/step11_findHamming_kmer.py 4-find_OriC/Input/text3.txt CGGAATCCTAAC 4
147 560 642 919 1468 1635 2140 2163 2355 2366 2497 2774 3861 4154 4311 4493 4860 5033 5267 5354 5484 6092 6572 6632 6786 7132 7941 8021 8031 8766 8978 8997 9021 9902 9991 10427 10683 10853 11077 11190 12359 12529 12562 12568 13067 13227 13407 13475 13784 13970 14001 14170 14433 15286 15978 16024 16530 16875 16999

为了适应题目糟糕的输入要求，不得不把输出修改为了字符串。

这里另外说个小技巧，也就是序列的解包：

也能满足输出。

除了知道位置，我们也可以加上计数的信息，这里太简单的我就不写了。

2.2-重置重复序列查找函数

回顾一下之前[[04-算法01-频繁出现的秘密]] 中的函数：

$ python3 02-k_mer_pattern.py ./Input/pattern01.txt 5
{'GTAGC', 'TAGCA'}

这里可否额外考虑错配信息呢？加上参数d，计算所有匹配的片段以及它们的mismatch，输出加上mismatch 后，计数最大的那个k-mer。

这里按照我的理解，所谓的包含一定容忍的k，比如指定d 的大小，指的是，先找到一些合适的k-mer，再从这些k-mer 池里，计算每个k-mer 及其容忍的k，最终合并容忍k 的总数记为各个k-mer 的计数。

而按照网站提供的伪代码提示，其思路为，遍历每个pattern，并找到这个pattern 对应的若干个符合d 容忍的一切pattern，并对这些pattern 全部计数：

不过经过思考后，我也发现，后者也存在其合理之后，也许存在某个k-mer，其全部的计数都是mismatch的。也就是说，某个pattern 的全部mismatch，都有可能是这个pattern 本来对应的序列。

至于代码怎么写，下篇再写吧。

参考资料

[1]

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

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

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

数据分析

编程算法

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

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

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

数据分析

编程算法

登录后参与评论

0 条评论

热度