【R语言双序列比对】全局比对Needleman-Wunsch算法&局部比对Smith-Waterman算法原理及代码实现

这篇具有很好参考价值的文章主要介绍了【R语言双序列比对】全局比对Needleman-Wunsch算法&局部比对Smith-Waterman算法原理及代码实现。希望对大家有所帮助。如果存在错误或未考虑完全的地方,请大家不吝赐教,您也可以点击"举报违法"按钮提交疑问。

目录

全局比对算法(Needleman-Wunsch)

原理

R代码实现

局部比对算法(Smith-Waterman)

原理

R代码实现

总结


全局比对算法(Needleman-Wunsch)

原理

其实这个跟数据结构学过的最短路径问题很像,核心思想就是依次寻求重复子问题的最优子结构。Needleman-Wunsch算法是一种全局联配算法,从整体上分析两个序列的关系,即考虑序列总长的整体比较,用类似于使整体相似最大化的方式,对序列进行联配。两个不等长度序列的联配分析必须考虑在一个序列中一些碱基的删除,即在另一序列做空位(Gap)处理。

【R语言双序列比对】全局比对Needleman-Wunsch算法&局部比对Smith-Waterman算法原理及代码实现,r语言,算法,开发语言

R代码实现
#全局比对(Needleman-Wunsch)
# 定义匹配、不匹配、gap的分数
library(stringr)
match_score <- 9  #匹配得分
mismatch_score <- -6  #不匹配得分
gap_penalty <- -2  #gap罚分
distance <- c(1,2,4)   #记录数据来源方向,分别代表对角、从上到下和从左到右
# 两个DNA序列
seq1 <- c("AACGTACTCAAGTCT")
seq2 <- c("TCGTACTCTAACGAT")

# 创建一个二维矩阵来保存比对得分
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)))  #创建一个得分矩阵,初始化全为0
route <- matrix(0, nrow = n+1, ncol = m+1, dimnames = list(c("0", seq1_n), c("0", seq2_n)))  #创建一个路径矩阵,初始化全为0

# 初始化得分矩阵(第一行和第一列,依次加gap罚分)
for (i in 2:(n + 1)) {
  score_matrix[i, 1] <- (i-1)*gap_penalty
  route[i, 1] <- distance[2]   #记录方向为从上到下
}
for (j in 2:(m + 1)) {
  score_matrix[1, j] <- (j-1)*gap_penalty
  route[1, j] <- distance[3]   #记录方向为从左到右
}

# 计算得分矩阵
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,    #走对角匹配
                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))  #看三种路径哪个值最大,存储这个方向来源
    for(z in index){
      route[i, j] <- route[i, j] + distance[z]  #把得到最大值的方向来源加和作为路径值
    }
  }
}
result <- list(score_matrix,route) #返回一个列表,包含得分矩阵和路径矩阵

#绘制得分矩阵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记录当前位置
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{
    #来源有对角,可能是对角(1),对角+上(3),对角+左(5),对角只用考虑是否匹配,不引入空缺
    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)
    }

    #来源有上方,可能是只有上方来的(2),或者上+对角(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]  #没有走对角,引入空缺,故下那个碱基仍然是该行碱基
      }
      seq_m2[y+1] <- "-" #由于是从上到下来源,所以第二列序列引入空缺
      path <- back_location(seq_m1, seq_m2, route, x - 1, y, path)
    }

    #来源有左方,同理,可能是左方(4),左+对角(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)
    }
  }
  #因为用了递归的方式,所以path更新后要传给上一层
  return(path)
}

p <- back_route(result[[2]])  #result[[2]]是list的第二个数据,即路径得分
print(p)

【R语言双序列比对】全局比对Needleman-Wunsch算法&局部比对Smith-Waterman算法原理及代码实现,r语言,算法,开发语言

 上图为全局比对的得分矩阵热图,可以看到右下角分数为83,即全局比对最优结果为83分,下面看看回溯结果。一共有15个最优比对结果。

【R语言双序列比对】全局比对Needleman-Wunsch算法&局部比对Smith-Waterman算法原理及代码实现,r语言,算法,开发语言【R语言双序列比对】全局比对Needleman-Wunsch算法&局部比对Smith-Waterman算法原理及代码实现,r语言,算法,开发语言

局部比对算法(Smith-Waterman)

原理

史密斯-沃特曼算法(Smith-Waterman algorithm)是一种进行局部序列比对(相对于全局比对)的算法,用于找出两个核苷酸序列或蛋白质序列之间的相似区域。该算法的目的不是进行全序列的比对,而是找出两个序列中具有高相似度的片段

该算法由坦普尔·史密斯(Temple F. Smith)和迈克尔·沃特曼(Michael S. Waterman)于1981年提出。Smith-Waterman算法是Needleman-Wunsch算法的一个变体,二者都是动态规划算法。这一算法的优势在于可以在给定的打分方法下找出两个序列的最优的局部比对(打分方法使用了置换矩阵和空位罚分)。该算法和Needleman-Wunsch算法的主要区别在于该算法不存在负分(负分被替换为零),因此局部比对成为可能。回溯从分数最高的矩阵元素开始,直到遇到分数为零的元素停止。分数最高的局部比对结果在此过程中产生。在实际运用中,人们通常使用该算法的优化版本。创建初始矩阵时初始化其首行和首列均为0。

R代码实现
##局部比对(Smith-Waterman)
# 定义匹配、不匹配、gap的分数
library(stringr)
match_score <- 9  #匹配得分
mismatch_score <- -6  #不匹配得分
gap_penalty <- -2  #gap罚分

# 两个DNA序列
seq1 <- c("AACGTACTCAAGTCT")
seq2 <- c("TCGTACTCTAACGAT")

# 创建一个二维矩阵来保存比对得分
n <- nchar(seq1)  #字符串长度
m <- nchar(seq2)
seq1_n <- strsplit(seq1, "")[[1]]  #得到单个碱基(character)
seq2_n <- strsplit(seq2, "")[[1]]
score_matrix <- matrix(0, n + 1, m + 1, dimnames = list(c("0", seq1_n), c("0", seq2_n)))  #创建一个得分矩阵,初始化全为0
#计算得分矩阵
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,  #走对角匹配
                  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)

# 回溯过程,递归
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  #定义当前位置
  aligned_seq1 <- ""  #定义两个空串
  aligned_seq2 <- ""  #定义两个空串

  while (score_matrix[current[1], current[2]] != 0) {
move <- choose_move(score_matrix, current[1], current[2])
#从对角来
    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))
}

# 选择移动方向
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")  #这里需要注意,回溯对角来源时要看是否匹配,否则返回结果除首尾匹配中间的碱基是平移的,即中间错配
  } else if (score_matrix[i, j] == (score_matrix[i - 1, j - 1] + mismatch_score) && seq1_n[i - 1] != seq2_n[j - 1]) {
    return("diagonal")  #不匹配也可能从对角来
  } else if (score_matrix[i, j] == (score_matrix[i - 1, j] + gap_penalty)) {
    return("up")
  } else {
    return("left")
  }
}
# 调用回溯函数
result <- backtrack(score_matrix, max_indices)
aligned_seq1 <- result[[1]]
aligned_seq2 <- result[[2]]

# 输出比对结果
print(aligned_seq1)
print(aligned_seq2)

【R语言双序列比对】全局比对Needleman-Wunsch算法&局部比对Smith-Waterman算法原理及代码实现,r语言,算法,开发语言

上图为局部比对的得分矩阵热图,可以看到得分最大值为93,即局部比对最优结果为93分,下面看看回溯结果。

【R语言双序列比对】全局比对Needleman-Wunsch算法&局部比对Smith-Waterman算法原理及代码实现,r语言,算法,开发语言


 除此之外,R自带函数也可以比对,调用也很容易,更改type="global"或者"local"即可实现全局和局部比对,代码如下:

install.packages("Biostrings")  #安装R包,它提供了用于序列操作、比对、模式匹配、转录组分析等功能的工具集合。
library(Biostrings)  #加载R包
seq1 <- DNAString("AACGTACTCAAGTCT")  #输入DNA序列
seq2 <- DNAString("TCGTACTCTAACGAT")   #输入DNA序列

#设置参数,匹配得分9分,不匹配得分-6分
mat <- nucleotideSubstitutionMatrix(match = 9, mismatch = -6, baseOnly = TRUE)  
View(mat)  #查看碱基匹配矩阵
#执行全局序列比对
Alignment <- pairwiseAlignment(seq1, seq2, substitutionMatrix = mat,gapOpening = 0, gapExtension = -2, type = “global”)
Alignment

比对结果如下,不同的是,回溯结果都只输出一个。

【R语言双序列比对】全局比对Needleman-Wunsch算法&局部比对Smith-Waterman算法原理及代码实现,r语言,算法,开发语言

【R语言双序列比对】全局比对Needleman-Wunsch算法&局部比对Smith-Waterman算法原理及代码实现,r语言,算法,开发语言

总结

全局比对和局部比对的显著区别是局部比对的得分矩阵中的分值均≥0,即在考虑打分方式时,局部比对要多考虑一种。在初始化得分矩阵的第一行和第一列时,两者有不同。对于全局比对,需要依次加上gap罚分值,而对于局部比对,则只需要都初始化为0。全局比对和局部比对的核心思想是依次寻求重复子问题的最优子结构,通过序列比对,可以发现序列之间的相似性,也可以判别序列之间的同源性,推断序列之间的进化关系,辨别序列之间的差别,寻找遗传变异。文章来源地址https://www.toymoban.com/news/detail-758114.html

到了这里,关于【R语言双序列比对】全局比对Needleman-Wunsch算法&局部比对Smith-Waterman算法原理及代码实现的文章就介绍完了。如果您还想了解更多内容,请在右上角搜索TOY模板网以前的文章或继续浏览下面的相关文章,希望大家以后多多支持TOY模板网!

本文来自互联网用户投稿,该文观点仅代表作者本人,不代表本站立场。本站仅提供信息存储空间服务,不拥有所有权,不承担相关法律责任。如若转载,请注明出处: 如若内容造成侵权/违法违规/事实不符,请点击违法举报进行投诉反馈,一经查实,立即删除!

领支付宝红包 赞助服务器费用

相关文章

  • 算法分析:C语言实现动态规划之最长公共子序列

    最长公共子序列问题:          下面的简单问题说明了动态规划的基本原理。在字母表一∑上,分别给出两个长度为n和m的字符串A和B,确定在A和B中最长公共子序列的长度。这里,A = a₁a₂...an。的子序列是一个形式为a₁ka₂k...aik的字符串,其中每个i都在1和n之间,并且

    2023年04月21日
    浏览(35)
  • 金豺(GJO)优化算法、matlab代码实现以及与PSO、GWO、SO算法比对

    提示:文章写完后,目录可以自动生成,如何生成可参考右边的帮助文档 金豺(GJO)优化算法是2022年由Nitish Chopra 等人提出,GJO 的灵感来自金豺 (Canis aureus) 的协作狩猎行为。算法的三个基本步骤是 猎物搜索、包围和突袭 。 该算法主要是由雄雌豺狼带领各个豺狼对猎物进行

    2024年02月02日
    浏览(91)
  • 一键AI高清换脸——基于InsightFace、CodeFormer实现高清换脸与验证换脸后效果能否通过人脸比对、人脸识别算法

    1、项目简介 AI换脸是指利用基于深度学习和计算机视觉来替换或合成图像或视频中的人脸。可以将一个人的脸替换为另一个人的脸,或者将一个人的表情合成到另一个人的照片或视频中。算法常常被用在娱乐目上,例如在社交媒体上创建有趣的照片或视频,也有用于电影制作

    2024年02月08日
    浏览(45)
  • C语言——局部和全局变量

    局部变量 定义在函数内部的变量称为局部变量(Local Variable) 局部变量的作用域(作用范围)仅限于函数内部, 离开该函数后是无效的 离开该函数后,局部变量自动释放 示例代码: 全局变量 在所有函数外部定义的变量称为全局变量(Global Variable),它的作用域默认是整个程

    2024年02月10日
    浏览(37)
  • 汇编调用C语言定义的全局变量

    在threadx移植中,系统的systick通过了宏定义的方式定义,很难对接库函数的时钟频率,不太利于进行维护 所以在C文件中自己定义了一个systick_Div的变量,通过宏定义方式设定systick的时钟频率 在汇编下要加载这个systick分频系数 方法: 总结:对汇编指令需要进一步熟悉。

    2024年02月15日
    浏览(42)
  • C语言变量和全局变量能否重名?

            局部变量和全局变量能否重名?         全局变量和局部变量是按照变量的作⽤域划分的。         简单地说,局部变量是定义在函数内部的变量;全局变量是定义在函数之外的变量。         全局变量可以为本⽂件中其他函数所共⽤。         局

    2024年01月20日
    浏览(40)
  • 算法题:摆动序列(贪心算法解决序列问题)

    这道题是一道贪心算法题,如果前两个数是递增,则后面要递减,如果不符合则往后遍历,直到找到符合的。(完整题目附在了最后) 代码如下: 完整题目: 376. 摆动序列 如果连续数字之间的差严格地在正数和负数之间交替,则数字序列称为  摆动序列 。 第一个差(如果

    2024年02月07日
    浏览(60)
  • c语言全局变量和局部变量问题汇总

    ✅作者简介:嵌入式领域优质创作者,博客专家 ✨个人主页:咸鱼弟 🔥系列专栏:单片机设计专栏  1、static的作用是什么?  定义静态变量  2、static有什么用途?(请至少说明两种)  1).限制变量的作用域(在程序的整个运行期间都不释放)  2).设置变量的存储域(存

    2024年02月06日
    浏览(41)
  • 【路径规划】全局路径规划算法——蚁群算法(含python实现)

    路径规划与轨迹跟踪系列算法 蚁群算法原理及其实现 蚁群算法详解(含例程) 图说蚁群算法(ACO)附源码 蚁群算法Python实现 蚁群算法(Ant Colony Algorithm, ACO) 于1991年首次提出,该算法模拟了自然界中蚂蚁的觅食行为。蚂蚁在寻找食物源时, 会在其经过的路径上释放一种信息

    2024年01月18日
    浏览(49)
  • 启发式搜索算法:A算法(全局、局部择优算法)+A*算法 解决八数码问题

    参考博客:人工智能搜索策略:A*算法 在图搜索算法中,如果能在搜索的每一步都利用估价函数f(n)=g(n)+h(n)对Open表中的节点进行排序,则该搜索算法为 A算法 。由于估价函数中带有问题自身的启发性信息,因此,A算法又称为启发式搜索算法。 对启发式搜索算法,又可根据搜

    2024年02月10日
    浏览(37)

觉得文章有用就打赏一下文章作者

支付宝扫一扫打赏

博客赞助

微信扫一扫打赏

请作者喝杯咖啡吧~博客赞助

支付宝扫一扫领取红包,优惠每天领

二维码1

领取红包

二维码2

领红包