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]
ps:事先说明,这一部分的代码这里存在一些问题。
至此为止,我们实现了以下功能:
但别忘了,这里我们始终都是对某个长度的序列进行计算的。开始的定义:某个长度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)
,可以得到对应输出吗?并输出它的位置信息?
这里的解法也很简单:
这里我并不想指定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 这两个就是互相挨着的序列,该如何避免这种情况呢?
再次强调,这一部分的代码这里存在一些问题。
上面的代码,我犯了非常严重的概念问题,我没有理解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