动态规划：基因序列比对

Posted 2021-04-18 程序员Hotel

基因序列比对（levenshtein算法）

中学学过生物学的人都会知道，生命的奥秘在DNA，它本质上是一个由碱基对组成的长长的序列，科学家们已经在20年前把智人的完整基因测定完毕了。自然，生命的奥秘远没有解开，测定只是第一步，相当于拿到了源代码，后面要查找子序列，寻找模式，互补反向变换，比对等等。大家都希望计算机程序能加速我们对生命的了解。
这些初步的工作，从编程的角度，就是文件和字符串操作。
我们先下载一个基因序列，我们去NCBI网站上下载一个基因序列，保存为一个文件如sequence.fasta。fasta是文本格式的，用于表示核酸序列或多肽序列的格式。其中核酸或氨基酸均以单个字母来表示，且允许在序列前添加序列名及注释。fasta格式首先以大于号“>”开头，接着是序列的标识符“gi|187608668|ref|NM_001043364.2|”，然后是序列的描述信息。换行后是序列信息，序列中允许空格，换行，空行，直到下一个大于号，表示该序列的结束。
我们可以看一下文件里面的内容（截取前十行）：

>NC_000006.12:c31170682-31164337 Homo sapiens chromosome 6, GRCh38.p13 Primary Assembly
GAGTAGTCCCTTCGCAAGCCCTCATTTCACCAGGCCCCCGGCTTGGGGCGCCTTCCTTCCCCATGGCGGG
ACACCTGGCTTCGGATTTCGCCTTCTCGCCCCCTCCAGGTGGTGGAGGTGATGGGCCAGGGGGGCCGGAG
CCGGGCTGGGTTGATCCTCGGACCTGGCTAAGCTTCCAAGGCCCTCCTGGAGGGCCAGGAATCGGGCCGG
GGGTTGGGCCAGGCTCTGAGGTGTGGGGGATTCCCCCATGCCCCCCGCCGTATGAGTTCTGTGGGGGGAT
GGCGTACTGTGGGCCCCAGGTTGGAGTGGGGCTAGTGCCCCAAGGCGGCTTGGAGACCTCTCAGCCTGAG
GGCGAAGCAGGAGTCGGGGTGGAGAGCAACTCCGATGGGGCCTCCCCGGAGCCCTGCACCGTCACCCCTG
GTGCCGTGAAGCTGGAGAAGGAGAAGCTGGAGCAAAACCCGGAGGAGGCAAGTGAGCTTCGACGGGGTTG
GGGTGTGGGGAGGTGGTCATGACAGGGCAGCCTGATGGGGAAGTGGTCACCTGCAGCTGCCCAGACCTGG
CACCCAGGAGAGGAGCAGGCAGGGTCAGCTGCCCTGGCCAGGGAGGGGTGTGTATCAACTGCTGGCAGCC
CTGGCAGGCAGGGGCCAGGTGGGAAGTGGAAGCTGGATTTCGAAGAGACAACTGCCGGTGAGGGCAGAGC

ATGC组成的序列就是生命的密码。我们现在全部拿到手上了，但是解密还任重道远。目前的阶段，还只能做一些基础性工作，进行零碎的解读。我相信，结合计算机科学，会大大加速我们对生命的理解。
有了这个文本的格式描述，我们可以先写一个函数从外部文件读取序列。程序如下：

def getgeneseq(filename):
    geneseq = ""
    file = open(filename, 'r')
    try:
        lines = file.readlines()
        for line in lines:
            if not line.startswith('>'):
                geneseq += line.rstrip()
    finally:
        file.close()

    return geneseq

程序以只读方式从外部打开一个文件，然后按行读取到一个列表中，通过首字符>号识别第一行的注释，之后把后面的行拼在一起成为一个大的字符串（rstrip()函数用于把右边的空格和换行去掉）
有了这个序列，我们可以计算一下核苷酸数量，程序如下：

# nucleotide counting
def ntcount(seq):
    counts = []
    for c in ['A', 'C', 'G', 'T']:
        counts.append(seq.count(c))
    return counts

这里，我们用到了count()方法，用于计算一个字符串中某字符出现的次数。
还可以将DNA表达翻译成RNA，程序如下：

# DNA -> RNA
def trans(seq):
    newseq = re.sub('T', 'U', seq)
    return newseq

这里，我们用到了正则表达式去替换T为U。
再进行一个反向变换，程序如下：

# reverse sequence
def reverse(seq):
    mapping = {
        "A":"T", "T":"A", "C":"G", "G":"C"
    }
    newseq = ""
    for c in seq:
        if c in mapping:
            newseq += mapping[c]
    return newseq

这里，我们用了一个字典，即ATGC的对应关系，程序对串中的每个字母，按照字典进行翻译，得到反向串。
有了基因序列，生物学上一个经常用到的功能就是比较两个序列之间的差异。
从计算机科学的角度，这个问题就转化成了两段字符串之间的距离。怎么理解“距离”这个词？我们可以规定一系列操作，比如增加一个字符、删除一个字符、替换一个字符，经过多次这样的操作后，一个串变成了另一个串，我们把操作的次数叫做“距离”。
如，对序列x：AGTAGTC和序列y：ATAGACC，我们可以这么操作：
把AGTAGTC的第二个字符G删除，变成ATAGTC；
把ATAGTC的倒数第二个字符T替换成A，变成ATAGAC；
把ATAGAC最后增加一个C，变成ATAGACC
经过三次操作，第一个序列变成了第二个序列。可以叫做距离为3。
当然，你也想到了，这种变换会有很多种可能性的组合，每种组合对应一个距离，我们把最小的距离叫做序列之间的差异，术语叫“Levenshtein距离”。一般来说，距离越小，两个序列的相似度越大。
我们看怎么求出这个Levenshtein距离。
猛地一拿到两个序列，想穷举一下，发现工作量实在太大了。我们人处理不了这么复杂的事情，那好，我们就想办法从小事情开始做起。
既然分析AGTAGTC和序列ATAGACC复杂，我们就分析小的序列。假定序列x只有一个字符A，序列y也只有一个字符A，这还是很简单的，两者一样，我们认为距离为0，不一样就认为距离为1。再看y的下一个字符T，跟x的A不一样，我们就在前面的距离基础上增加1。再继续看y的下一个字符A，跟x的一样，距离不变。
如果序列只有一个字符，那么分析都是很简单的。
更通用的分析，我们现在假定要比较的是x的i字符和y的j字符，怎么办呢？我们想象一下，积然我们是从第一个字符一个一个比较的，每次都是算出了距离，然后往下再比较，那到现在的位置上，只有三个变换的来源，第一个是从位置[i-1,j]变换来（相当于编辑增加），第二个是i,j-1，第三个是i-1,j-1。
我们画一个表格示意：

我们把序列x和序列y逐字符对应成一个大表格。对[i,j]这个格子，它来源于临近三个格子的变换：上边、左边和左上。如果Xi=Yj，那距离没变，仍然是[i-1,j-1]格子上的距离。如果Xi<>Yj，那就进行变换，从临近的三个格子变化而来，取最小值+1即可。
所以按照这个规则，[i,j]处要填上2。
按照这个思路编程：

def Levenshtein(X, Y):
    n = len(X)
    m = len(Y)
    D=[[0]*m for i in range(n)] #n*m 2-dimension array
    for j in range(m): #process 1st row
        if Y[j]==X[0]:
            if j>0:
                D[0][j] = D[0][j-1]
            else:
                D[0][j] = 0
        else:
            if j>0:
                D[0][j] = D[0][j-1]+1
            else:
                D[0][j] = 1
    for i in range(n): #process 1st column
        if Y[0]==X[i]:
            if i>0:
                D[i][0] = D[i-1][0]
            else:
                D[i][0] = 0 
        else:
            if i>0:
                D[i][0] = D[i-1][0]+1
            else:
                D[i][0] = 1

    for i in range(1,n): #process remaining rows
        for j in range(1,m): #process one row
            if X[i]==Y[j]:
                D[i][j] = D[i-1][j-1]
            else:
                D[i][j] = min(D[i][j-1],D[i-1][j],D[i-1][j-1])+1 #递推

    return D

我们测试一下：

strx="AGTAGTC"
stry="ATAGACC"
print(Levenshtein(strx, stry))

运行结果：

[[0, 1, 1, 2, 2, 3, 4], 
[1, 1, 2, 1, 2, 3, 4], 
[2, 1, 2, 2, 2, 3, 4], 
[2, 2, 1, 2, 2, 3, 4], 
[3, 3, 2, 1, 2, 3, 4], 
[4, 3, 3, 2, 2, 3, 4], 
[5, 4, 4, 3, 3, 2, 3]]

输出了表格中的每一格。我们只要看最后那个数，就知道了这两个序列之间的距离为3。
其实，Python提供了现成的Levenshtein包，你直接import Levenshtein就可以了。
实际上，生物信息学中经常要进行序列比对，很多时候用到Needleman-Wunsch 算法和Smith-Waterman算法，思路也是类似的。
对基因序列的解读才刚开始，解读生命源代码的工作才迈出第一步，生命的奥秘还隐藏在幽暗之中，期待下一个牛顿的降生，为人类带来光明。
上面介绍的这个算法应用比较广，unix里面有一个命令diff就是用的同样的动态规划的思路。

点击“阅读原文”,可进入原创书籍，阅读目录