目录

  • 序列对比过程中的罚分规则
  • 选择的序列
  • 名称
  • 具体的序列info
  • DNA的dotplot实现
  • 采用蛋白质进行dotplot
  • 使用矩阵进行打分（积分+罚分）
  • BLOSUM62的规则
  • 空位罚分
  • 最优化（optimization）
  • 使用needle软件进行在线global对比。
  • 本地实现打分运算
  • 调用BIO库进行本地运算
  • 手动计算
  • 方式
  • 一些局限性
  • 手动计算过程
  • 局部的序列的次优比对
  • 局部次优比对的运算结果
  • 附录
  • 相关引用

Reference URL
https://biopython.org/DIST/docs/api/Bio.pairwise2-module.html （Bio官方文档）
https://www.biostars.org/p/307285/ （论坛关于Bio中Align部分的问题）
http://www.ryxxff.com/91470.html （Biopython序列比对简要概括）

序列对比过程中的罚分规则

python 双序列比对 双序列比对回溯_html

选择的序列

名称

选择的序列是冠状病毒envelop的gene序列

  • NC_045512.2:26245-26472 Severe acute respiratory syndrome coronavirus 2 isolate Wuhan-Hu-1, envelop genome
  • NC_002645.1:24750-24983 Human coronavirus 229E, envelop genome

具体的序列info

>NC_002645.1:24750-24983 Human coronavirus 229E, envelop genome
ATGTTCCTTAAGCTAGTGGATGATCATGCTTTGGTTGTTAATGTACTACTCTGGTGTGTGGTGCTTATAG
TGATACTACTAGTGTGTATTACAATAATTAAACTAATTAAGCTTTGTTTCACTTGCCATATGTTTTGTAA
TAGAACAGTTTATGGCCCCATTAAAAATGTGTACCACATTTACCAATCATATATGCACATAGACCCTTTC
CCTAAACGAGTTATTGATTTCTAA
>NC_045512.2:26245-26472 Severe acute respiratory syndrome coronavirus 2 isolate Wuhan-Hu-1, envelop genome
ATGTACTCATTCGTTTCGGAAGAGACAGGTACGTTAATAGTTAATAGCGTACTTCTTTTTCTTGCTTTCG
TGGTATTCTTGCTAGTTACACTAGCCATCCTTACTGCGCTTCGATTGTGTGCGTACTGCTGCAATATTGT
TAACGTGAGTCTTGTAAAACCTTCTTTTTACGTTTACTCTCGTGTTAAAAATCTGAATTCTTCTAGAGTT
CCTGATCTTCTGGTCTAA

DNA的dotplot实现

  1. 在线的server端实现
  2. 本地化（未完成）
URL：http://www.bioinformatics.nl/emboss-explorer/
网站概览：

python 双序列比对 双序列比对回溯_生物信息学_02

设定“word size = 2”，得到的结果如下，并不能观察到直观的关系。

python 双序列比对 双序列比对回溯_python 双序列比对_03


设定“word size = 3”，得到的结果如下，仍然不能观察到直观的关系。但是已经比之前要稀疏很多了。

python 双序列比对 双序列比对回溯_生物信息学_04


小结：实际上如果只用gene序列的话，可能出现太宽泛的问题，因为很多的蛋白都有不止一个codon，所以我们进一步采用对蛋白质进行探究。

采用蛋白质进行dotplot

使用Biopython库对原序列进行翻译，得到结果，重新进行dotplot。阈值设定为2。

python 双序列比对 双序列比对回溯_python 双序列比对_05

使用矩阵进行打分（积分+罚分）

BLOSUM62的规则

python 双序列比对 双序列比对回溯_html_06


两条序列通过BLOSUM的打分政策，通过积分和罚分的方式，最终得到两条序列的最终打分，中间的打分过程则需要用到最优化（optimal processing）。

空位罚分

由于不可能所有的匹配都是按顺序的，所以中间需要适时进行空位罚分，从而进行一些Global内部的Local化。以下是两种罚分的方式。

python 双序列比对 双序列比对回溯_python_07

最优化（optimization）

使用needle软件进行在线global对比。

python 双序列比对 双序列比对回溯_最优化_08

本地实现打分运算

以上的结果都是基于在线端实现的，下面我们使用了本地化的Bio库进行了运算，可以弥补批量运算的缺陷。但是经过实地测试，如果出现长序列的pairwise运算，对内存的压力依旧很大。

我们一开始使用了两端全序列进行运算，而不是选择中间的一部分Gene进行运算，内存消耗到一定的量之后程序中断，递出报错。

python 双序列比对 双序列比对回溯_python 双序列比对_09

调用BIO库进行本地运算

经过查询我们发现对于比对，已经有很多人提供了大量的库可以使用，这里我们选择了使用范围较广的biopython库。

biopython中提供了大量的矩阵和选项供使用。相关的代码请见附录。

python 双序列比对 双序列比对回溯_html_10


相关的简要使用说明

"""
The names of the alignment functions in this module follow the
convention
<alignment type>XX
where <alignment type> is either "global" or "local" and XX is a 2
character code indicating the parameters it takes.  The first
character indicates the parameters for matches (and mismatches), and
the second indicates the parameters for gap penalties.

The match parameters are::

    CODE  DESCRIPTION
    x     No parameters. Identical characters have score of 1, otherwise 0.
    m     A match score is the score of identical chars, otherwise mismatch
          score.
    d     A dictionary returns the score of any pair of characters.
    c     A callback function returns scores.

The gap penalty parameters are::

    CODE  DESCRIPTION
    x     No gap penalties.
    s     Same open and extend gap penalties for both sequences.
    d     The sequences have different open and extend gap penalties.
    c     A callback function returns the gap penalties.
"""

这里我们选择globaldd方法：

# arguments: seq1, seq2, matrix, gap_penalty_of_seq1, extendpenalty_of_seq1, gap_penalty_of_seq2, extendpenalty_of_seq2 
pairwise2.align.globaldd(str(seq1).strip("*"), str(seq2).strip("*"), matrix, -10, -.5, -10, -.5)

所得结果和线上接近（可能是因为删除终止符号的原因？），如下所示：

MFLKLVDDHA-LVVNVLLWCVVLIVILLVCITIIKLIKLCFTCHMFCNRTVYGPIKNVYHIYQSYMHIDPFPKRVIDF--
|......... |.||..|......|.|||...|.....||..|   ||.......|.....|.......  ..||.|.  
MYSFVSEETGTLIVNSVLLFLAFVVFLLVTLAILTALRLCAYC---CNIVNVSLVKPSFYVYSRVKNLN--SSRVPDLLV
  Score=61

手动计算

方式

  1. 扩列
  2. 带有方向性的顺序计算
  3. 回溯计算（traceback）

一些局限性

查询字典的过程中需要三个步骤：

  1. 得到需求（查询矩阵）
  2. 翻阅字典（BLOSUM62）
  3. 得到结果（填充矩阵）

平均三个过程至少需要2s，一共耗时6s。假设我们的有m×n大小的对比矩阵，如果m和n的大小差距不大，则假定中间值为x。

查询BLOSUM62的过程中，需要的耗时约为6*x(s)。蛋白质最少需要大约20个左右的残基才能组成一个具有生物功能的3D结构，如果我们忽略初始的行列，最少的矩阵大小为20×20 ，那么至少需要 400×6 = 2400(s) = 40(min)才能完成矩阵的运算。

以上的步骤只是包含了查询的过程，实际的计分过程中还涉及到一些逻辑运算（简单的比大小工作）。所以需要消耗的时间可能远超过1个小时。

python 双序列比对 双序列比对回溯_python_11


但是手动的计算能够加深对整体公式的了解，所以我把原来的序列截断一部分，仅仅留下其中匹配较好的一部分加以计算和最优化。

手动计算过程

以下过程是通过全局比对的方式进行，但是回溯的路径上符合两个条件：

  1. 没有负值
  2. 最后的值：17是所有值中的最大值

所以可以认为此结果也是局部比对（Waterman&Smith Algorithm）的结果。

python 双序列比对 双序列比对回溯_最优化_12

局部的序列的次优比对

局部序列的最优比对原理

python 双序列比对 双序列比对回溯_python_13

局部次优比对的运算结果

计算过程仿照之前，但是由于我们的提取序列和之前的序列不太合适，如果直接盲目采用Local运算，并且取次优结果的话，得到的结果会非常不可靠（匹配度非常理想的两条序列往往不能采取次优运算）。

附录

# !/usr/bin/env/python3
# -*- coding:UTF-8 -*-
"""
# Time:2020年3月31日11:01:16
# Author:zhengjm
"""

import os
import Bio
import Bio.SeqIO
import Bio.Align
from Bio import pairwise2
from Bio.pairwise2 import format_alignment
from Bio.SubsMat import MatrixInfo as matlist


matrix = matlist.blosum62
path = os.getcwd()

file = path + '/compare.fasta'
compare = Bio.SeqIO.parse(file, format='fasta')
seq_list = []
for i in compare:
    seq_list.append(i.seq)
    # print(i.seq)

seq1 = seq_list[0]
seq2 = seq_list[1]


# caution: end sign (e.g."*" ) should not appear in seq.
# penalty should be negative
# arguments(seq1, seq2, matrix, gap_penalty_of_seq1, extendpenalty_of_seq1, gap_penalty_of_seq2, extendpenalty_of_seq2 )
test_new = pairwise2.align.globaldd(str(seq1).strip("*"), str(seq2).strip("*"), matrix, -10, -.5, -10, -.5)
for test in test_new:
    print(format_alignment(*test))

相关引用

[1]佚名. Different alignment results between Emboss Needle and Biopython pairwise2EB/OL[2020–04–02]. https://www.biostars.org/p/307285/.
[2]陈铭. 生物信息学第二版陈铭.pdf[M]. .