其實這是通識課作業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'); ?>