如何编写高效R代码找到最小化MSE的最优种子序列

给定基础数据与函数 #

参考序列与种子模式 #

ref_seq      <- "MGHQQLYWSHPRKFGQGSRSCRVTSNRHGLIRKYGLNMSRQSFR"
seed_pattern <- "FKDHKHIDVKDRHRTRHLAK??????????"

种子模式包含10个通配符?，需替换为标准氨基酸。

氨基酸计数函数 #

# 归一化氨基酸计数（频率）
aa_count_normalized <- function(x) {
  AADict <- c(
    "A", "R", "N", "D", "C", "E", "Q", "G", "H",
    "I", "L", "K", "M", "F", "P", "S", "T", "W", "Y", "V"
  )
  AAC <- summary(factor(strsplit(x, split = "")[[1]], levels = AADict),
                 maxsum = 21
  ) / nchar(x)
  AAC
}

# 原始氨基酸计数（数量）
aa_count <- function(x) {
  AADict <- c(
    "A", "R", "N", "D", "C", "E", "Q", "G", "H",
    "I", "L", "K", "M", "F", "P", "S", "T", "W", "Y", "V"
  )
  AAC <- summary(factor(strsplit(x, split = "")[[1]], levels = AADict),
                 maxsum = 21
  ) 
  AAC
}

归一化参考序列氨基酸计数 #

为匹配种子模式的长度，将参考序列的氨基酸频率按种子长度缩放，得到目标计数：

# 按种子模式长度归一化参考序列氨基酸计数
refseq_aa_content <- aa_count_normalized(ref_seq) * nchar(seed_pattern) 
refseq_aa_content
# 输出结果：
#         A         R         N         D         C         E 
# 0.0000000 4.7727273 1.3636364 0.0000000 0.6818182 0.0000000 
#         Q         G         H         I         L         K 
# 2.7272727 3.4090909 2.0454545 0.6818182 2.0454545 1.3636364 
#         M         F         P         S         T         W 
# 1.3636364 1.3636364 0.6818182 4.0909091 0.6818182 0.6818182 
#         Y         V 
# 1.3636364 0.6818182  

需求说明 #

需保留种子模式中的非通配符，用AADict中的氨基酸替换10个通配符，使得最终序列的氨基酸计数与refseq_aa_content的均方误差（MSE）最小。MSE计算函数如下：

mse <- function (ref, new_seq) {
  return(mean((ref - new_seq)^2))
}

已知最优替换方案为：3个G、2个Q、1个R、3个S、1个Y，最终序列为FKDHKHIDVKDRHRTRHLAKRQQGGGSSSY。

高效R实现代码 #

# 1. 统计种子模式中非通配符的氨基酸数量
seed_fixed_count <- aa_count(gsub("\\?", "", seed_pattern))

# 2. 计算每个氨基酸需要补充的目标数量（参考计数 - 现有计数）
target_add <- refseq_aa_content - seed_fixed_count
# 补充数量不能为负（无法减少已有氨基酸）
target_add[target_add < 0] <- 0

# 3. 贪心分配10个通配符：每次选添加后MSE最小的氨基酸
remaining <- 10
add_counts <- rep(0, length(target_add))
names(add_counts) <- names(target_add)

while (remaining > 0) {
  # 模拟添加每个氨基酸后的MSE
  mse_values <- sapply(names(target_add), function(aa) {
    temp_add <- add_counts
    temp_add[aa] <- temp_add[aa] + 1
    new_counts <- seed_fixed_count + temp_add
    mse(refseq_aa_content, new_counts)
  })
  # 选MSE最小的氨基酸添加
  best_aa <- names(which.min(mse_values))
  add_counts[best_aa] <- add_counts[best_aa] + 1
  remaining <- remaining - 1
}

# 4. 生成替换字符串并替换通配符
replace_str <- paste(rep(names(add_counts), add_counts), collapse = "")
final_seq <- gsub("\\?", replace_str, seed_pattern)

# 输出结果
cat("最优序列：", final_seq, "\n")
cat("添加的氨基酸计数：\n")
print(add_counts[add_counts > 0])

代码说明 #

• 先统计种子模式现有氨基酸数量，计算每个氨基酸的目标补充量；
• 用贪心算法逐步分配通配符，每次选择能最小化MSE的氨基酸添加；
• 小数量分配场景下，该方法高效且能得到最优解。

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