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