Sequence Alignment :: Similarity - dynamic programming

其實這是通識課作業XD

用來比對DNA、mRNA、cDNA序列相似度的方法，計分方式為相同的鹼基+1分、不同的-1、對上空白則-2，兩個序列不改變序列順序而允許插入空白(gap;null)的情況下做排列組合比較取出最高分者。計算出來的分數可以代表這兩個序列的相似度，越高分相似度就越高，例如ACG對上AC可以是(以-代表gap) (1)ACG對AC-：1+1-2=0、(2)ACG對A-C：1-2-1=-2、(3)ACG對-AC：-2-1-1=-4。

這好像叫做Needleman-Wunsch algorithm，但本來他是用暴力算法，就是同時列出兩個序列並包含gap的可能組合做計分，後來發現利用prefix計分來推算分數可以降低運算複雜度。

	null	A	G
null	(0,0)	(0,1)	(0,2)
A	(1,0)	(1,1)	(1,2)
C	(2,0)	(2,1)	(2,2)
G	(3,0)	(3,1)	(3,2)

	null	A	G
null	0	-2	-4
A	-2	1	-1
C	-4	-1	0
G	-6	-3	0

比較序列s：ACG與序列t：AG；例如在(2,1)就是比較(null)A與(null)AC得分為-1。p,s null對上null為0分 比較快的推算法就是說，(i,j)可以從(i-1,j)、(i-1,j-1)、(i,j-1)推來，因為(i-1,j)、(i-1,j-1)、(i,j-1)都已經是最佳解了，所以只要取他們中推算到(i,j)最高的那個分數。

以PHP實做的程式碼，另外還有個豪華物件版 XD可以輸出prefix觀察比對過程 <?php function Similarity_csv_lite($s1,$s2) { $datas=array(); $csv=', '; //建立初始值 $datas[0][0]=0; for ($i=1;$i<=strlen($s1);$i++) { $datas[$i][0]=-2*$i; } for ($j=1;$j<=strlen($s2);$j++) { $datas[0][$j]=-2*$j; $csv.=','.$s2[$j-1]; } //開始推算 for ($i=1;$i<=strlen($s1);$i++) { for ($j=1;$j<=strlen($s2);$j++) { $tmp_value=0; $tmp_p=($s1[$i-1]==$s2[$j-1])?1:-1; $datas[$i][$j]=$datas[$i-1][$j]-2;; //(i-1,j) $tmp_value=$datas[$i-1][$j-1]+$tmp_p; //(i-1,j-1) $datas[$i][$j]=($tmp_value>$datas[$i][$j])?$tmp_value:$datas[$i][$j]; $tmp_value=$datas[$i][$j-1]-2; //(i,j-1) $datas[$i][$j]=($tmp_value>$datas[$i][$j])?$tmp_value:$datas[$i][$j]; } } $csv.="\r\n ,".join(',',$datas[0]); for ($i=1;$i<=strlen($s1);$i++) { $csv.="\r\n".$s1[$i-1].','.join(',',$datas[$i]); } return $csv; } echo Similarity_csv_lite('GACGGATTAG','GATCGGAATAG'); ?>