【ROSALIND】【练Python，学生信】21 确定限制酶酶切位点
如果第一次阅读本系列文档请先移步阅读【ROSALIND】【练Python，学生信】00 写在前面 谢谢配合~
题目：
确定限制酶酶切位点（Locating Restriction Sites）
Given: A DNA string of length at most 1 kbp in FASTA format.
所给：一条长度不超过1kb的DNA序列。
Return: The position and length of every reverse palindrome in the string having length between 4 and 12. You may return these pairs in any order.
需得：这条序列上长度在4-12nt间的回文序列的位置和长度，以任意顺序输出。
测试数据
>Rosalind_24
TCAATGCATGCGGGTCTATATGCAT
测试输出
4 6
5 4
6 6
7 4
17 4
18 4
20 6
21 4
背景
限制性核酸内切酶是可以识别并附着特定的脱氧核苷酸序列，并在每条链中特定部位的两个脱氧核糖核苷酸之间的磷酸二酯键进行切割的一类酶，简称限制酶。限制酶识别DNA序列中的回文序列，回文序列是一种旋转对称结构，在轴的两侧序列相同而反向，特点是在该段的碱基序列的互补链之间正读反读都相同。
限制酶的生理作用是限制酶降解外源DNA，维护宿主遗传稳定的保护机制。限制酶在几乎所有细菌的种、属中都有发现，是细菌防御噬菌体侵染的机制。
思路
本题只需抓住回文序列的特点即可，即互补链反向阅读结果与正链相同。因此挨个扫描序列中长度在4-12之间的片段，与反向互补序列比较，相同即为酶切位点序列。
代码
def readfasta(lines):
seq = []
index = []
seqplast = ""
numlines = 0
for i in lines:
if '>' in i:
index.append(i.replace("\n", "").replace(">", ""))
seq.append(seqplast.replace("\n", ""))
seqplast = ""
numlines += 1
else:
seqplast = seqplast + i.replace("\n", "")
numlines += 1
if numlines == len(lines):
seq.append(seqplast.replace("\n", ""))
seq = seq[1:]
return index, seq
f = open(rosalind_revp.txt', 'r')
lines = f.readlines()
f.close()
[index, seq] = readfasta(lines)
seq = seq[0]
c = ""
for i in seq:
if i == 'A':
c = c + 'T'
elif i == 'G':
c = c + 'C'
elif i == 'T':
c = c + 'A'
elif i == 'C':
c = c + 'G'
rseq = c # rseq中存储原序列的互补序列
i = 0
p = []
l = []
while i < len(seq)-4+1: # 逐个扫描序列，因长度最小为4，所以到倒数第四个即可
j = 4
if len(seq) - i >= 12:
while j <= 12: # 当未扫描序列还很长时，只需扫描到开始位点后的12位
s1 = seq[i:i+j]
s2 = rseq[i:i+j]
s2 = s2[::-1] # 序列反向
if s1 == s2: # 满足回文序列
p.append(i+1)
l.append(j)
j += 1
else:
while j <= len(seq) - i: 当未扫描序列长度不足12，扫到序列结尾
s1 = seq[i:i+j]
s2 = rseq[i:i+j]
s2 = s2[::-1]
if s1 == s2:
p.append(i+1)
l.append(j)
j += 1
i += 1
f = open('output.txt', 'a')
i = 0
while i < len(p):
f.write(str(p[i]) + " " + str(l[i]) + '\n') # 按格式输出到文件
i += 1
f.close()
