目錄
全局比對(duì)算法（Needleman-Wunsch）
原理
R代碼實(shí)現(xiàn)
局部比對(duì)算法（Smith-Waterman）
原理
R代碼實(shí)現(xiàn)
總結(jié)
全局比對(duì)算法（Needleman-Wunsch）
原理
其實(shí)這個(gè)跟數(shù)據(jù)結(jié)構(gòu)學(xué)過(guò)的最短路徑問(wèn)題很像，核心思想就是依次尋求重復(fù)子問(wèn)題的最優(yōu)子結(jié)構(gòu)。Needleman-Wunsch算法是一種全局聯(lián)配算法，從整體上分析兩個(gè)序列的關(guān)系，即考慮序列總長(zhǎng)的整體比較，用類似于使整體相似最大化的方式，對(duì)序列進(jìn)行聯(lián)配。兩個(gè)不等長(zhǎng)度序列的聯(lián)配分析必須考慮在一個(gè)序列中一些堿基的刪除，即在另一序列做空位(Gap)處理。
R代碼實(shí)現(xiàn)
#全局比對(duì)（Needleman-Wunsch）
# 定義匹配、不匹配、gap的分?jǐn)?shù)
library(stringr)
match_score <- 9 #匹配得分
mismatch_score <- -6 #不匹配得分
gap_penalty <- -2 #gap罰分
distance <- c(1,2,4) #記錄數(shù)據(jù)來(lái)源方向,分別代表對(duì)角、從上到下和從左到右
# 兩個(gè)DNA序列
seq1 <- c("AACGTACTCAAGTCT")
seq2 <- c("TCGTACTCTAACGAT")
# 創(chuàng)建一個(gè)二維矩陣來(lái)保存比對(duì)得分
n <- nchar(seq1)
m <- nchar(seq2)
seq1_n <- strsplit(seq1, "")[[1]]
seq2_n <- strsplit(seq2, "")[[1]]
score_matrix <- matrix(0, n + 1, m + 1, dimnames = list(c("0", seq1_n), c("0", seq2_n))) #創(chuàng)建一個(gè)得分矩陣，初始化全為0
route <- matrix(0, nrow = n+1, ncol = m+1, dimnames = list(c("0", seq1_n), c("0", seq2_n))) #創(chuàng)建一個(gè)路徑矩陣，初始化全為0
# 初始化得分矩陣（第一行和第一列，依次加gap罰分）
for (i in 2:(n + 1)) {
score_matrix[i, 1] <- (i-1)*gap_penalty
route[i, 1] <- distance[2] #記錄方向?yàn)閺纳系较?}
for (j in 2:(m + 1)) {
score_matrix[1, j] <- (j-1)*gap_penalty
route[1, j] <- distance[3] #記錄方向?yàn)閺淖蟮接?}
# 計(jì)算得分矩陣
for (i in 2:(n + 1)) {
for (j in 2:(m + 1)) {
match_value <- ifelse(substr(seq1, i - 1, i - 1) == substr(seq2, j - 1, j - 1), match_score, mismatch_score)
scores <- c(score_matrix[i - 1, j - 1] + match_value, #走對(duì)角匹配
score_matrix[i - 1, j] + gap_penalty, #從上到下引入空缺
score_matrix[i, j - 1] + gap_penalty) #從左到右引入空缺
score_matrix[i, j] <- max(scores) #取三種路徑最大值
index <- which(scores == max(scores)) #看三種路徑哪個(gè)值最大，存儲(chǔ)這個(gè)方向來(lái)源
for(z in index){
route[i, j] <- route[i, j] + distance[z] #把得到最大值的方向來(lái)源加和作為路徑值
}
}
}
result <- list(score_matrix,route) #返回一個(gè)列表，包含得分矩陣和路徑矩陣
#繪制得分矩陣retype
library(pheatmap)
pheatmap(score_matrix,
color = colorRampPalette(c("lightpink", "red","purple"))(256),
display_numbers = score_matrix,
number_color = "white",
cluster_rows = FALSE,
cluster_cols = FALSE,
show_rownames = TRUE,
show_colnames = TRUE,
fontsize = 13,
legend = TRUE,
fontsize_number = 13)
#回溯算法
back_route <- function(route){
seq1 <- rownames(route)
seq2 <- colnames(route)
path <- list(0)
#從最右下角的值開始回溯
path <- back_location(seq1, seq2, route, length(seq1), length(seq2), path)
path[[1]] <- NULL
return(path)
}
#x和y記錄當(dāng)前位置
back_location <- function(seq1, seq2, route, x, y, path){
#回溯的一種匹配模式，將其加入path
if(x == 1 & y == 1){
#匹配到左上角
seq1 <- str_c(c( seq1[2:length(seq1)]), collapse = "")
seq2 <- str_c(c( seq2[2:length(seq2)]), collapse = "")
path[[length(lengths(path))+1]] <- rbind(seq1, seq2)
}
else{
#來(lái)源有對(duì)角，可能是對(duì)角（1），對(duì)角+上（3），對(duì)角+左（5），對(duì)角只用考慮是否匹配，不引入空缺
if(route[x, y] == 1 || route[x, y] == 3 || route[x, y] == 5){
seq_m1 <- seq1
seq_m2 <- seq2
path <- back_location(seq_m1, seq_m2, route, x - 1, y - 1, path)
}
#來(lái)源有上方，可能是只有上方來(lái)的（2），或者上+對(duì)角（3），或者上+左（6）
if(route[x, y] == 2 || route[x, y] == 3 || route[x, y] == 6){
seq_m1 <- seq1
seq_m2 <- seq2
#在seq2的y+1位置插入“-”
for(i in length(seq_m2):min((y+1), length(seq_m2))){
seq_m2[i+1] <- seq_m2[i] #沒有走對(duì)角，引入空缺，故下那個(gè)堿基仍然是該行堿基
}
seq_m2[y+1] <- "-" #由于是從上到下來(lái)源，所以第二列序列引入空缺
path <- back_location(seq_m1, seq_m2, route, x - 1, y, path)
}
#來(lái)源有左方，同理，可能是左方（4），左+對(duì)角（5）、左+上（6）
if(route[x, y] == 4 || route[x, y] == 5 || route[x, y] == 6){
seq_m1 <- seq1
seq_m2 <- seq2
#在seq1的x+1位置插入“-”
for(i in length(seq_m1):min((x+1), length(seq_m1))){
seq_m1[i+1] <- seq_m1[i]
}
seq_m1[x+1] <- "-"
path <- back_location(seq_m1, seq_m2, route, x, y - 1, path)
}
}
#因?yàn)橛昧诉f歸的方式，所以path更新后要傳給上一層
return(path)
}
p <- back_route(result[[2]]) #result[[2]]是list的第二個(gè)數(shù)據(jù)，即路徑得分
print(p)
?上圖為全局比對(duì)的得分矩陣熱圖，可以看到右下角分?jǐn)?shù)為83，即全局比對(duì)最優(yōu)結(jié)果為83分，下面看看回溯結(jié)果。一共有15個(gè)最優(yōu)比對(duì)結(jié)果。
局部比對(duì)算法（Smith-Waterman）
原理
史密斯-沃特曼算法（Smith-Waterman algorithm）是一種進(jìn)行局部序列比對(duì)（相對(duì)于全局比對(duì)）的算法，用于找出兩個(gè)核苷酸序列或蛋白質(zhì)序列之間的相似區(qū)域。該算法的目的不是進(jìn)行全序列的比對(duì)，而是找出兩個(gè)序列中具有高相似度的片段。
該算法由坦普爾·史密斯（Temple F. Smith）和邁克爾·沃特曼（Michael S. Waterman）于1981年提出。Smith-Waterman算法是Needleman-Wunsch算法的一個(gè)變體，二者都是動(dòng)態(tài)規(guī)劃算法。這一算法的優(yōu)勢(shì)在于可以在給定的打分方法下找出兩個(gè)序列的最優(yōu)的局部比對(duì)（打分方法使用了置換矩陣和空位罰分）。該算法和Needleman-Wunsch算法的主要區(qū)別在于該算法不存在負(fù)分（負(fù)分被替換為零），因此局部比對(duì)成為可能?；厮輳姆?jǐn)?shù)最高的矩陣元素開始，直到遇到分?jǐn)?shù)為零的元素停止。分?jǐn)?shù)最高的局部比對(duì)結(jié)果在此過(guò)程中產(chǎn)生。在實(shí)際運(yùn)用中，人們通常使用該算法的優(yōu)化版本。創(chuàng)建初始矩陣時(shí)初始化其首行和首列均為0。
R代碼實(shí)現(xiàn)
##局部比對(duì)（Smith-Waterman）
# 定義匹配、不匹配、gap的分?jǐn)?shù)
library(stringr)
match_score <- 9 #匹配得分
mismatch_score <- -6 #不匹配得分
gap_penalty <- -2 #gap罰分
# 兩個(gè)DNA序列
seq1 <- c("AACGTACTCAAGTCT")
seq2 <- c("TCGTACTCTAACGAT")
# 創(chuàng)建一個(gè)二維矩陣來(lái)保存比對(duì)得分
n <- nchar(seq1) #字符串長(zhǎng)度
m <- nchar(seq2)
seq1_n <- strsplit(seq1, "")[[1]] #得到單個(gè)堿基（character）
seq2_n <- strsplit(seq2, "")[[1]]
score_matrix <- matrix(0, n + 1, m + 1, dimnames = list(c("0", seq1_n), c("0", seq2_n))) #創(chuàng)建一個(gè)得分矩陣，初始化全為0
#計(jì)算得分矩陣
for (i in 2:(n + 1)) {
for (j in 2:(m + 1)) {
match_value <- ifelse(substr(seq1, i - 1, i - 1) == substr(seq2, j - 1, j - 1), match_score, mismatch_score)
scores <- c(score_matrix[i - 1, j - 1] + match_value, #走對(duì)角匹配
score_matrix[i - 1, j] + gap_penalty, #從上到下引入空缺
score_matrix[i, j - 1] + gap_penalty, #從左到右引入空缺
0) #還有一種0
score_matrix[i, j] <- max(scores) #取四種路徑最大值
}
}
#繪制得分矩陣熱圖
library(pheatmap)
pheatmap(score_matrix,
color = colorRampPalette(c("lightpink", "red","purple"))(256),
display_numbers = score_matrix,
number_color = "white",
cluster_rows = FALSE,
cluster_cols = FALSE,
show_rownames = TRUE,
show_colnames = TRUE,
fontsize = 13,
legend = TRUE,
fontsize_number = 13)
# 回溯過(guò)程，遞歸
max_score <- max(score_matrix) #找到得分矩陣最大值位置
max_indices <- which(score_matrix == max_score, arr.ind = TRUE)
backtrack <- function(score_matrix, max_indices) {
nrow <- nrow(score_matrix)
ncol <- ncol(score_matrix)
current <- max_indices #定義當(dāng)前位置
aligned_seq1 <- "" #定義兩個(gè)空串
aligned_seq2 <- "" #定義兩個(gè)空串
while (score_matrix[current[1], current[2]] != 0) {
move <- choose_move(score_matrix, current[1], current[2])
#從對(duì)角來(lái)
if (move == "diagonal") {
aligned_seq1 <- paste(seq1_n[current[1] - 1], aligned_seq1, sep = "")
aligned_seq2 <- paste(seq2_n[current[2] - 1], aligned_seq2, sep = "")
current <- c(current[1] - 1, current[2] - 1)
}
#從上到下
else if (move == "up") {
aligned_seq1 <- paste(seq1_n[current[1] - 1], aligned_seq1, sep = "")
aligned_seq2 <- paste("-", aligned_seq2, sep = "")
current <- c(current[1] - 1, current[2])
}
#從左到右
else {
aligned_seq1 <- paste("-", aligned_seq1, sep = "")
aligned_seq2 <- paste(seq2_n[current[2] - 1], aligned_seq2, sep = "")
current <- c(current[1], current[2] - 1)
}
}
return(list(aligned_seq1, aligned_seq2))
}
# 選擇移動(dòng)方向
choose_move <- function(score_matrix, i, j) {
if (score_matrix[i, j] == (score_matrix[i - 1, j - 1] + match_score) && seq1_n[i - 1] == seq2_n[j - 1]) {
return("diagonal") #這里需要注意，回溯對(duì)角來(lái)源時(shí)要看是否匹配，否則返回結(jié)果除首尾匹配中間的堿基是平移的，即中間錯(cuò)配
} else if (score_matrix[i, j] == (score_matrix[i - 1, j - 1] + mismatch_score) && seq1_n[i - 1] != seq2_n[j - 1]) {
return("diagonal") #不匹配也可能從對(duì)角來(lái)
} else if (score_matrix[i, j] == (score_matrix[i - 1, j] + gap_penalty)) {
return("up")
} else {
return("left")
}
}
# 調(diào)用回溯函數(shù)
result <- backtrack(score_matrix, max_indices)
aligned_seq1 <- result[[1]]
aligned_seq2 <- result[[2]]
# 輸出比對(duì)結(jié)果
print(aligned_seq1)
print(aligned_seq2)
上圖為局部比對(duì)的得分矩陣熱圖，可以看到得分最大值為93，即局部比對(duì)最優(yōu)結(jié)果為93分，下面看看回溯結(jié)果。
?除此之外，R自帶函數(shù)也可以比對(duì)，調(diào)用也很容易，更改type="global"或者"local"即可實(shí)現(xiàn)全局和局部比對(duì)，代碼如下：
install.packages("Biostrings") #安裝R包，它提供了用于序列操作、比對(duì)、模式匹配、轉(zhuǎn)錄組分析等功能的工具集合。
library(Biostrings) #加載R包
seq1 <- DNAString("AACGTACTCAAGTCT") #輸入DNA序列
seq2 <- DNAString("TCGTACTCTAACGAT") #輸入DNA序列
#設(shè)置參數(shù)，匹配得分9分，不匹配得分-6分
mat <- nucleotideSubstitutionMatrix(match = 9, mismatch = -6, baseOnly = TRUE)
View(mat) #查看堿基匹配矩陣
#執(zhí)行全局序列比對(duì)
Alignment <- pairwiseAlignment(seq1, seq2, substitutionMatrix = mat,gapOpening = 0, gapExtension = -2, type = “global”)
Alignment比對(duì)結(jié)果如下，不同的是，回溯結(jié)果都只輸出一個(gè)。
文章來(lái)源：http://www.zghlxwxcb.cn/news/detail-758114.html
總結(jié)
全局比對(duì)和局部比對(duì)的顯著區(qū)別是局部比對(duì)的得分矩陣中的分值均≥0，即在考慮打分方式時(shí)，局部比對(duì)要多考慮一種。在初始化得分矩陣的第一行和第一列時(shí)，兩者有不同。對(duì)于全局比對(duì)，需要依次加上gap罰分值，而對(duì)于局部比對(duì)，則只需要都初始化為0。全局比對(duì)和局部比對(duì)的核心思想是依次尋求重復(fù)子問(wèn)題的最優(yōu)子結(jié)構(gòu)，通過(guò)序列比對(duì)，可以發(fā)現(xiàn)序列之間的相似性，也可以判別序列之間的同源性，推斷序列之間的進(jìn)化關(guān)系，辨別序列之間的差別，尋找遺傳變異。文章來(lái)源地址http://www.zghlxwxcb.cn/news/detail-758114.html
到了這里，關(guān)于【R語(yǔ)言雙序列比對(duì)】全局比對(duì)Needleman-Wunsch算法&局部比對(duì)Smith-Waterman算法原理及代碼實(shí)現(xiàn)的文章就介紹完了。如果您還想了解更多內(nèi)容，請(qǐng)?jiān)谟疑辖撬阉鱐OY模板網(wǎng)以前的文章或繼續(xù)瀏覽下面的相關(guān)文章，希望大家以后多多支持TOY模板網(wǎng)！
