如何在4k条氨基酸序列中检测错义突变与移码突变

需求说明 #

我有近4000条同一蛋白的不同长度氨基酸序列，均为AAstrings类对象，需要检测其中的错义突变和移码突变（包括插入、缺失）。目前通过Stack Overflow的代码实现了错义突变检测，但该代码仅支持等长序列比对，无法识别导致序列长度差异的插入/缺失位点及具体突变信息。

期望输出格式 #

结果需包含突变ID、参考氨基酸、样本氨基酸及突变位置，示例如下：

#>      ID Reference_AA Sample_AA Pos
...
#>15 query9           F         -  16
...
#>30 query10          H         QN  8

若需将插入拆分为单独条目，格式如下：

#>30 query10          -         Q   8
#>31 query10          H         N   9

初始数据示例 #

seqs <- AAStringSet(c("FLGKIWPSYKGRPGNF", "FLGKIWPSHKGRPGNF", "FLGRIWPSHKGRPGNF", "FLGKIWPSHKGRPGNF", "FLGKIWPSHKGRPGNF", "FLGKVWPSHKGRPGNF", "FLGKVWPSHKGRPGNF", "FLGKIWPSHKGRPGNF", "FLGKIWPSHKGRPGN", "FLGKIWPSQNKGRPGNF"))
ref <- seqs[1]

解决方案代码 #

1. 加载依赖包 #

# 安装并加载Biostrings包
if (!requireNamespace("Biostrings", quietly = TRUE)) {
  if (!requireNamespace("BiocManager", quietly = TRUE)) {
    install.packages("BiocManager")
  }
  BiocManager::install("Biostrings")
}
library(Biostrings)

2. 预处理序列数据 #

# 给序列命名作为突变ID
names(seqs) <- paste0("query", 1:length(seqs))

3. 定义突变提取函数 #

合并插入版（对应第一个示例格式）

get_mutations <- function(sample_seq, ref_seq, sample_id) {
  # 全局比对，允许插入缺失，设置gap罚分
  align <- pairwiseAlignment(ref_seq, sample_seq, type = "global", gapOpening = 10, gapExtension = 1)
  
  # 提取比对后的序列并转为字符向量
  ref_aligned <- strsplit(as.character(pattern(align)), "")[[1]]
  sample_aligned <- strsplit(as.character(subject(align)), "")[[1]]
  
  mutations <- list()
  ref_pos <- 1  # 参考序列的实际位置（跳过gap）
  align_pos <- 1  # 比对后的索引位置
  
  while (align_pos <= length(ref_aligned)) {
    ref_char <- ref_aligned[align_pos]
    sample_char <- sample_aligned[align_pos]
    
    if (ref_char != sample_char) {
      # 缺失：参考有氨基酸，样本为gap
      if (ref_char != "-" && sample_char == "-") {
        mutations[[length(mutations)+1]] <- data.frame(
          ID = sample_id,
          Reference_AA = ref_char,
          Sample_AA = "-",
          Pos = ref_pos,
          stringsAsFactors = FALSE
        )
        ref_pos <- ref_pos + 1
        align_pos <- align_pos + 1
      } 
      # 插入：参考为gap，样本有连续氨基酸，合并显示
      else if (ref_char == "-" && sample_char != "-") {
        insert_chars <- c()
        while (align_pos <= length(ref_aligned) && ref_aligned[align_pos] == "-" && sample_aligned[align_pos] != "-") {
          insert_chars <- c(insert_chars, sample_aligned[align_pos])
          align_pos <- align_pos + 1
        }
        mutations[[length(mutations)+1]] <- data.frame(
          ID = sample_id,
          Reference_AA = "-",
          Sample_AA = paste(insert_chars, collapse = ""),
          Pos = ref_pos,
          stringsAsFactors = FALSE
        )
      } 
      # 错义突变：两者均非gap且不同
      else {
        mutations[[length(mutations)+1]] <- data.frame(
          ID = sample_id,
          Reference_AA = ref_char,
          Sample_AA = sample_char,
          Pos = ref_pos,
          stringsAsFactors = FALSE
        )
        ref_pos <- ref_pos + 1
        align_pos <- align_pos + 1
      }
    } else {
      # 序列匹配，移动位置
      if (ref_char != "-") ref_pos <- ref_pos + 1
      align_pos <- align_pos + 1
    }
  }
  
  # 合并结果，无突变则返回空表
  if (length(mutations) > 0) do.call(rbind, mutations) else {
    data.frame(ID = character(), Reference_AA = character(), Sample_AA = character(), Pos = integer(), stringsAsFactors = FALSE)
  }
}

拆分插入版（对应第二个示例格式）

若需将插入的每个氨基酸单独成行，修改函数中插入处理部分：

# 替换原函数中的插入处理代码
else if (ref_char == "-" && sample_char != "-") {
  # 拆分插入的每个氨基酸为单独条目
  while (align_pos <= length(ref_aligned) && ref_aligned[align_pos] == "-" && sample_aligned[align_pos] != "-") {
    mutations[[length(mutations)+1]] <- data.frame(
      ID = sample_id,
      Reference_AA = "-",
      Sample_AA = sample_aligned[align_pos],
      Pos = ref_pos,
      stringsAsFactors = FALSE
    )
    align_pos <- align_pos + 1
  }
}

4. 批量处理所有序列 #

# 循环处理所有样本序列，收集突变结果
all_mutations <- lapply(2:length(seqs), function(i) {
  get_mutations(seqs[i], ref, names(seqs)[i])
})

# 合并所有结果为数据框
all_mutations_df <- do.call(rbind, all_mutations)

# 查看结果
print(all_mutations_df)

代码说明 #

• 使用pairwiseAlignment进行全局序列比对，支持插入缺失，可通过调整gapOpening和gapExtension参数优化比对结果
• 自动区分错义突变、缺失（参考有AA/样本gap）、插入（参考gap/样本有AA）三种突变类型
• 两种格式可选：插入氨基酸合并显示或拆分为单独条目，满足不同分析需求

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