R语言中如何匹配变异肽段序列对应的参考氨基酸序列？

核心逻辑 #

针对每个变异肽段，遍历所有参考序列，计算参考序列中所有与变异肽段等长的子串（含正向、反向两种匹配方向）与变异肽段的匹配字符数，取匹配数最高的参考序列作为匹配结果。

实现代码 #

# 辅助函数：计算两个等长字符串的匹配字符数
count_match <- function(str1, str2) {
  str1_chr <- strsplit(str1, "")[[1]]
  str2_chr <- strsplit(str2, "")[[1]]
  sum(str1_chr == str2_chr)
}

# 主匹配流程
pep_length <- nchar(Variants$AA_seq[1])
Variants$Seq_name <- sapply(Variants$AA_seq, function(pep) {
  # 计算当前肽段与所有参考序列的最高匹配得分
  ref_scores <- sapply(Ref_seq$AA_seq, function(ref_seq) {
    ref_length <- nchar(ref_seq)
    max_score <- 0
    # 滑动窗口遍历参考序列的所有等长子串
    for (i in 1:(ref_length - pep_length + 1)) {
      sub_fwd <- substr(ref_seq, i, i + pep_length - 1)
      score_fwd <- count_match(pep, sub_fwd)
      # 计算反向匹配得分
      sub_rev <- paste(rev(strsplit(sub_fwd, "")[[1]]), collapse = "")
      score_rev <- count_match(pep, sub_rev)
      current_max <- max(score_fwd, score_rev)
      if (current_max > max_score) max_score <- current_max
    }
    return(max_score)
  })
  # 返回得分最高的参考序列名称
  Ref_seq$Seq_name[which.max(ref_scores)]
})

# 查看匹配结果
print(Variants)

运行结果 #

和你给出的预期输出完全一致：

peptideID     AA_seq Seq_name
1      Pep1 QEISALVKYF     Ref1
2      Pep2 HTGERGNLVT     Ref5
3      Pep3 NKMTTSVLIK     Ref2
4      Pep4 SMNLKNDYPD     Ref9
5      Pep5 NEPGYSQSTI     Ref9
6      Pep6 NPQDVIMVKL     Ref6
7      Pep7 MAAKFNKMTL     Ref2
8      Pep8 RRQKDPSSGT     Ref3
9      Pep9 QQQWTELFSV     Ref7

优化提示 #

如果待匹配的序列量级较大，可直接调用Biostrings包的pairwiseAlignment函数实现批量快速比对，无需自行实现滑动窗口逻辑，运行效率更高。

内容的提问来源于stack exchange，提问作者Douglas