简单双序列比对实现

本文转载自查看原文 2020-03-31 00:43 935 生物信息学

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序列比对简要概括）
https://www.jianshu.com/p/f9179c517d4e（山东大学生物信息学教程）

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

两种基本罚分手段

选择的序列

名称

选择的序列是冠状病毒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实现

在线的server端实现
本地化（未完成）

URL：http://www.bioinformatics.nl/emboss-explorer/
网站概览：

设定“word size = 2”，得到的结果如下，并不能观察到直观的关系。
size值为2的结果
设定“word size = 3”，得到的结果如下，仍然不能观察到直观的关系。但是已经比之前要稀疏很多了。
在这里插入图片描述
小结：实际上如果只用gene序列的话，可能出现太宽泛的问题，因为很多的蛋白都有不止一个codon，所以我们进一步采用对蛋白质进行探究。

采用蛋白质进行dotplot

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

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

BLOSUM62的规则

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

空位罚分

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

最优化（optimization）

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

参考https://www.cnblogs.com/zhengjm/p/12602522.html

本地实现打分运算

以上的结果都是基于在线端实现的，下面我们使用了本地化的Bio库进行了运算，可以弥补批量运算的缺陷。但是经过实地测试，如果出现长序列的pairwise运算，对内存的压力依旧很大。
我们一开始使用了两端全序列进行运算，而不是选择中间的一部分Gene进行运算，内存消耗到一定的量之后程序中断，递出报错。
示意图

调用BIO库进行本地运算

经过查询我们发现对于比对，已经有很多人提供了大量的库可以使用，这里我们选择了使用范围较广的biopython库。
biopython中提供了大量的矩阵和选项供使用。相关的代码请见附录。
各种矩阵的示意
相关的简要使用说明

"""
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

手动计算

方式

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

一些局限性

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

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

平均三个过程至少需要2s，一共耗时6s。假设我们的有m×n大小的对比矩阵，如果m和n的大小差距不大，则假定中间值为x。
查询BLOSUM62的过程中，需要的耗时约为6*x(s)。蛋白质最少需要大约20个左右的残基才能组成一个具有生物功能的3D结构，如果我们忽略初始的行列，最少的矩阵大小为20×20 ，那么至少需要 400×6 = 2400(s) = 40(min)才能完成矩阵的运算。
以上的步骤只是包含了查询的过程，实际的计分过程中还涉及到一些逻辑运算（简单的比大小工作）。所以需要消耗的时间可能远超过1个小时。

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

手动计算过程

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

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

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

局部的序列的次优比对

局部序列的最优比对原理

局部次优比对的运算结果

计算过程仿照之前，但是由于我们的提取序列和之前的序列不太合适，如果直接盲目采用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))