序列比对正文

一前言

DNA蛋白质序列比对是生物信息学中的基本手段，深入了解比对算法的核心内容是学习掌握生物信息学的必需内容。本专题通过 Perl 程序分别实现动态规划算法与模糊匹配的方式（即字符串相似度算法），处理比较多条序列，从而比较两种算法的区别和优劣。进而通过实际操作（程序运行）了解动态规划算法与其它算法之间的区别，掌握动态规划算法的特点和优势。

二本论

2.1动态规划算法

动态规划解决序列比对问题的基本思想:使用迭代法计算出两个序列的相似分值, 并存入一个得分矩阵中, 根据得分矩阵回溯寻找最优的比对序列.

序列比对是动态规划的一个重要应用。序列比对问题通常是使用编辑操作（替换、插入、删除一个要素等）进行序列转换。每次操作对应不同成本，目标是找到编辑序列的最低成本。可以很自然地想到使用递归解决这个问题，序列A到B的最优编辑通过以下措施之一实现：插入B的第一个字符，对A和B的剩余序列进行最优比对；删去A的第一个字符，对A和B进行最优比对；用B的第一个字符替换A的第一个字符，对A的剩余序列和B进行最优比对。局部比对可在矩阵中列表表示，单元(i,j)表示A[1..i]到b[1..j]最优比对的成本[1]。单元(i,j)的成本计算可通过累加相邻单元的操作成本并选择最优解实现。

两条序列分别作为横坐标和纵坐标放置，组成一个路径矩阵，即得分矩阵, 矩阵元素(i,j)值为比对的得分值。在得分矩阵中到达位置为(i,j)的某一个元素有三种可能的路径：通过位置i- 1,j- 1 的对角方向, 没有空位罚分：通过列 j 的垂直方向, 通过行i 的水平方向, 空位罚分的值取决于插入空格的个数。序列比对序列同源(homology)指的是序列来自相同的祖先, 意味着这些序列具有相同的进化历史, 而序列的相似性(similarity)指的是两序列在某参数条件下的相像, 它可以用相同残基的百分比或是其他的方法来表示[2]。列之间的相似度是可以量化的参数，列是否同源需要有进化事实的验证, 显著的相似性通常意味着同源是运用某种特定的数学模型或算法, 找出两个或多个序列之间的最大匹配碱基或残基数, 比对算法的结果在很大程度上反映了序列之间的相似性程度以及它们的生物学特征。

2.2 动态规划序列比对算法

动态规划算法比较序列相似度程序如下：

open(INFILE, 'C:\Users\Administrator\Desktop\生信比对算法应用专题比对序列.txt) or die "Can't open\生信比对算法应用专题比对序列.txt $!";

o pen(OUTFILE, '>C:\Users\Administrator\Desktop\guozhanqiang.txt') or die "Can't open guozhanqiang.txt: $!";

while()

{ if ( $_ =~ / SEQUENCE1:/ )

{ $_=~s/ SEQUENCE1://g;

@case1=split("",$_); }

if ( $_ =~ / SEQUENCE2:/ )

{ $_=~s/ SEQUENCE2://g;

@case2=split("",$_); }

} close (INFILE);

%HoH=(); #初始化"双维关联数组"

$blosumFile='C:\Users\Administrator\Desktop\BLOSUM62.txt';

&BLOSUM; #从blosum文件，赋值"双维关联数组"%HoH

$lenght=$i; # 得出两个序列的长度

sub dynamicprogramming{ # 动态规划函数

@compare=@_; #将传入的数组赋给compare

for ($t=0; $t<=$lenght; $t++)

{$guo[0][$t]=(-4)*$t; }

for ($j=0; $j<=$lenght; $j++)

{$guo[$j][0]=(-4)*$j;} #初始化得分矩阵

for ($i = 1; $i <= $lenght; $i++)

{ for ($j = 1; $j <=$lenght; $j++)

{ if ($case1[$i-1] eq $compare[$j-1]) #如果两个元素匹配，则斜着走

{ $guo[$i][$j]=$guo[$i-1][$j-1]+$HoH{$case1[$i]}{$compare[$j]};

$b[$i][$j]=1; }

esif ($guo[$i-1][$j-1]>= $guo[$i-1][$j]-4 && $guo[$i-1][$j-1]>= $guo[$i][$j-1]-4 )

{ $guo[$i][$j]=$guo[$i-1][$j-1]+$HoH{$case1[$i]}{$compare[$j]};

$b[$i][$j]=1; }

elsif ($guo[$i-1][$j]-4>=$guo[$i][$j-1]-4)

{ $guo[$i][$j]=$guo[$i-1][$j]-4;

$b[$i][$j]=2; }

else { $guo[$i][$j]=$guo[$i][$j-1]+$HoH{$case1[$i]}{$compare[$j]}-4;

$b[$i][$j]=3; } } }

$guo[$lenght][$lenght];

}

sub BLOSUM{

open(BLOSUM, $blosumFile) or die "Can't open data.txt: $!";

while () {

next if /^#/; #忽略所有以 # 开始的行

#然后以空格为标准，劈开每一行，放入双维关联数组

($AA, $A, $R, $N, $D, $C, $Q, $E, $G, $H, $I, $L, $K, $M, $F, $P, $S, $T, $W, $Y, $V, $B, $Z, $X, $m ) = split( ' ' );

$HoH{ $AA }{ 'A' } = $A; $HoH{ $AA }{ 'R' } = $R;

$HoH{ $AA }{ 'N' } = $N; $HoH{ $AA }{ 'D' } = $D; $HoH{ $AA }{ 'C' } = $C; $HoH{ $AA }{ 'Q' } = $Q; $HoH{ $AA }{ 'E' } = $E; $HoH{ $AA }{ 'G' } = $G; $HoH{ $AA }{ 'H' } = $H; $HoH{ $AA }{ 'I' } = $I; $HoH{ $AA }{ 'L' } = $L; $HoH{ $AA }{ 'K' }