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端實現
- 本地化（未完成）
設定“word size = 2”，得到的結果如下，並不能觀察到直觀的關系。
設定“word size = 3”，得到的結果如下，仍然不能觀察到直觀的關系。但是已經比之前要稀疏很多了。
小結：實際上如果只用gene序列的話，可能出現太寬泛的問題，因為很多的蛋白都有不止一個codon，所以我們進一步采用對蛋白質進行探究。
采用蛋白質進行dotplot
使用Biopython庫對原序列進行翻譯，得到結果，重新進行dotplot。閾值設定為2。
使用矩陣進行打分（積分+罰分）
BLOSUM62的規則
兩條序列通過BLOSUM的打分政策，通過積分和罰分的方式，最終得到兩條序列的最終打分，中間的打分過程則需要用到最優化（optimal processing）。
空位罰分
由於不可能所有的匹配都是按順序的，所以中間需要適時進行空位罰分，從而進行一些Global內部的Local化。以下是兩種罰分的方式。
最優化（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))
相關引用
[1]佚名. Different alignment results between Emboss Needle and Biopython pairwise2EB/OL[2020–04–02]. https://www.biostars.org/p/307285/.
[2]陳銘. 生物信息學第二版陳銘.pdf[M]. .