原創: hxj7
一般而言,運用動態規劃演算法進行序列比對對記憶體空間的要求是 o(mn) 階的,本文介紹了一種線性空間要求的序列比對方法。前文如《序列比對(一)全域性比對needleman-wunsch演算法》所介紹的運用動態規劃演算法進行序列比對時,對記憶體空間的要求是o(mn) 階的。如果只是要求最大得分而不要求最大得分所對應的序列(也就是不進行回溯)的話,其實只需要保留上一行的得分就夠了,這一點可以從計算得分矩陣的迭代公式中看出來。
引自但是如果要求回溯呢,是否有一種線性空間演算法來進行序列比對呢?前人已經給出了多種演算法。其中一種在《生物序列分析》一書中給出,描述如下:
如圖中所說,關鍵點就是找到v值,然後通過不斷的分劃,最終得到全部的比對序列。本文給出了這種演算法的一種**實現。
**的關鍵在於終止條件的設定以及必要時巧妙地顛倒行列。與 o(mn) 階的演算法相比,這種演算法只能得到其中一種最佳比對方式,而無法得到所有的可能。
具體的**如下:
(由於**中運用了「引用(具體地就是指**中的 int &n 這一用法)」這一方法,所以是以cpp檔案編譯的)
#include
#include
#include
#define maxseq 1000
#define gap_char '-'
#define max2(a, b) ((a) > (b) ? (a) : (b))
#define exchange(a, b, c)
// 對空位的罰分是線性的
struct unit
;typedef
struct unit *punit;
void
strupper
(char
*s);
float
getfscore
(char a,
char b)
;float
dividealign
(punit*
* a,
int ti,
int tj,
int bi,
int bj,
char
* saln,
char
* raln,
int&n)
;void
align
(void);
char
* sraw;
char
* rraw;
float gap =
-2.5
;// 對空位的罰分
float tm[3]
;int
main()
printf
("the 1st seq: ");
scanf
("%s"
, sraw)
;printf
("the 2nd seq: ");
scanf
("%s"
, rraw)
;align()
;free
(sraw)
;free
(rraw)
;return0;
}void
strupper
(char
*s) s++;}
}// 替換矩陣:match分值為5,mismatch分值為-4
// 陣列下標是兩個字元的ascii碼減去65之後的和
float fmatrix=
;float
getfscore
(char a,
char b)
float
dividealign
(punit*
* a,
int ti,
int tj,
int bi,
int bj,
char
* saln,
char
* raln,
int&n)
if(ti == bi)
return
(bj - tj)
* gap;}if
(tj == bj)
return
(bi - ti)
* gap;
}// 為什麼有時要顛倒行列?
// 因為在bi - ti = 1這種情況下有可能會出現u = ti, v = tj從而導致無限迴圈
// 如果顛倒行列,可以將分治進行下去,因為此時不可能出現u = ti, v = tj的情況。
flag =0;
if(bi - ti <=1)
u =(ti + bi)/2
;// 中間行/列
for(j = tj; j <= bj; j++
) a[0]
[j]->m =
(j - tj)
* gap;
k =1;
// 當前行
for(i = ti +
1; i <= bi; i++
, k =
1- k)}}
v = a[
1- k]
[bj]
->v;
maxscore = a[
1- k]
[bj]
->m;
if(flag)
dividealign
(a, ti, tj, u, v, saln, raln, n)
;dividealign
(a, u, v, bi, bj, saln, raln, n)
;return maxscore;
}void
align
(void
)for
(i =
0; i <=
1; i++
)for
(j =
0; j <= r; j++)}
}if((salign =
(char*)
malloc
(sizeof
(char)*
(m + n)))
==null
|| \
(ralign =
(char*)
malloc
(sizeof
(char)*
(m + n)))
==null
)// 將字串都變成大寫
strupper
(sraw)
;strupper
(rraw)
;// 遞迴比對
l =0;
maxscore =
dividealign
(aunit,0,
0, m, n, salign, ralign, l)
;// 列印比對結果
printf
("max score: %f\n"
, maxscore)
;for
(i =
0; i < l; i++
)printf
("\n");
for(i =
0; i < l; i++
)printf
("\n");
// 釋放記憶體
free
(salign)
;free
(ralign)
;for
(i =
0; i <=
1; i++
)free
(aunit)
;}
2序列比對問題
str1 abdad str2 bacd 兩字串進行序列比對,定義乙個可以用來衡量比對效能的得分函式 令f x,y 表示x與y比對的得分。假設x和y都是字元,如果x與y相同,那麼f x,y 2,如果x與y不同,那麼f x,y 1,如果x或y是 那麼f x,y 1。str1和str2的2序列比對問題是...
生物序列區域性比對之Blast演算法
演算法基本原理 blast演算法是 1990 年由altschul 等人提出的兩序列區域性比對演算法,採用了一種短片段匹配演算法和一種有效的統計模型來找出目的序列和資料庫之間的最佳區域性比對效果。blast 演算法是一種基於區域性序列比對的序列比對演算法。廣泛被使用在蛋白質 dna序列的分析問題中,...
雙序列比對的基礎之PAM矩陣
pam矩陣的記分方法是基於蛋白序列中單點可接受 point accepted mutation,pam 的概念,通過對蛋白質進化模式的研究而建立的。pam矩陣是由dayhoff等人構建了與71個家族的序列關聯的假想系統發育樹,其中每對序列間的差異不超過它們殘基總數的15 用簡約法建樹,統計相似序列比...