AI复现Nature biotechnology同款，只需1min

pythonic生物人2026年08月21日约 6 分钟

大家好，我是pythonic生物人。

文献：McCreary, J., Hemez, C.F., Raymond, M.H.et al. Evolution of botulinum neurotoxin serotype X proteases to induce inflammatory cell death in cancer cells. Nat Biotechnol(2026). https://doi.org/10.1038/s41587-026-03243-9.

原文中的Fig. 3子图bc展示的内容：

子图b,野生型BoNT/X 在VAMP1来源的17-mer肽段上的底物富集图谱，该肽段已比对至 procaspase-1的底物序列背景中。黄色线标示预测剪切位点。
子图c,突变体 X(PC) 在用于进化的17-mer procaspase-1底物上的底物富集图谱。黄色线标示预测剪切位点，黄色框标示天然氨基酸。星号表示终止密码子。

其中，这两张图本质上是饱和突变结果的二维投影：

X轴，底物肽链的17个位置（Position）。
Y轴，每个位置被替换成的20种天然氨基酸 + 1个终止密码子（*）。
Fill填充色，log2(Enrichment Score)。

深红，该突变体在筛选后完全丢失，意味着蛋白酶完全不接受这个突变。
深蓝，该突变体被显著富集，意味着蛋白酶偏好这个突变。
白色，中性突变，不影响切割。

AI + base r实现

AI提示词：模拟数据，复现s41587-026-03243-9中的Fig. 3子图bc，使用base r实现。

模拟数据

# 1. 参数设置 ---------------------------------------------------------------
set.seed(3243)

amino_acids <- c("*", "A", "C", "D", "E", "F", "G", "H", "I", "K",
"L", "M", "N", "P", "Q", "R", "S", "T", "V", "W", "Y")
vamp1_sequence <- strsplit("LERDQKLSELDDRADAL", "")[[1]]
caspase1_sequence <- strsplit("NLSLPTTEEFEDDAIKK", "")[[1]]

# 数值范围为 log2(富集倍数)；更负表示该替换更不易被蛋白酶接受。
score_limits <- c(-6, 1)
heat_colors <- colorRampPalette(c("#df1644", "#ed7091", "#f7cbd5",
"#ffffff", "#c7d8e9", "#719bc7"))(256)
native_color <- "#f4c62b"
cleavage_color <- "#f2d400"

# 2. 模拟数据 ---------------------------------------------------------------
simulate_profile <- function(sequence, sensitivity) {
  n_rows <- length(amino_acids)
  n_cols <- length(sequence)
  score <- matrix(NA_real_, nrow = n_rows, ncol = n_cols,
                  dimnames = list(amino_acids, sequence))

# sensitivity 越大，该位点的非天然替换越容易呈红色（负富集）。
for (column in seq_len(n_cols)) {
    score[, column] <- 0.65 -
      sensitivity[column] * runif(n_rows, 0.75, 1.35) +
      rnorm(n_rows, sd = 0.28)
  }

# 终止密码子通常导致最强负富集；天然氨基酸用黄色方框标记。
  score["*", ] <- -6
for (column in seq_len(n_cols)) {
    score[sequence[column], column] <- runif(1, -0.25, 0.35)
  }
  score[score < score_limits[1]] <- score_limits[1]
  score[score > score_limits[2]] <- score_limits[2]
  score
}

# BoNT/X 主要限制切割位点 R|A 周围；X(PC) 对 E|D 两侧约六个位点更严格。
vamp1_sensitivity <- c(0.8, 0.7, 0.35, 0.55, 0.30, 0.25, 0.80, 0.65,
0.35, 0.75, 1.00, 1.40, 4.20, 3.80, 3.00, 0.75, 0.30)
caspase1_sensitivity <- c(0.10, 0.15, 0.20, 0.25, 0.35, 0.35, 0.50, 1.20,
2.20, 3.20, 4.80, 4.60, 4.00, 3.00, 1.70, 0.55, 0.25)

vamp1_score <- simulate_profile(vamp1_sequence, vamp1_sensitivity)
caspase1_score <- simulate_profile(caspase1_sequence, caspase1_sensitivity)

# 加入少量具有生化合理性的"可容忍替换"，形成与原图相近的局部浅色方格。
vamp1_score[c("A", "C", "G", "S"), 14] <- c(0.25, 0.05, 0.10, -0.15)
vamp1_score[c("R", "K"), 13] <- c(0.30, -0.45)
caspase1_score[c("D", "E", "Q"), 11] <- c(-0.25, 0.30, -0.60)
caspase1_score[c("D", "E", "N"), 12] <- c(0.25, -0.25, -0.70)

绘制热图

# 4. 绘图函数 ---------------------------------------------------------------
score_to_color <- function(values) {
  scaled <- (values - score_limits[1]) / diff(score_limits)
  index <- 1 + floor(scaled * (length(heat_colors) - 1))
  heat_colors[pmax(1, pmin(length(heat_colors), index))]
}

draw_profile <- function(score, sequence, panel_label, heading,
                         substrate_label, cleavage_after) {
  n_rows <- nrow(score)
  n_cols <- ncol(score)
  plot.new()
  plot.window(xlim = c(-3.7, n_cols + 3.5),
              ylim = c(0.1, n_rows + 3.1), xaxs = "i", yaxs = "i")

# 热图主体：第一行（终止密码子）显示在最上方。
for (column in seq_len(n_cols)) {
for (row in seq_len(n_rows)) {
      y <- n_rows - row + 1
      rect(column - 0.5, y - 0.5, column + 0.5, y + 0.5,
           col = score_to_color(score[row, column]),
           border = "#66717a", lwd = 0.45)
    }
  }

# 天然残基黄色框和预测切割位点黄色竖线。
for (column in seq_len(n_cols)) {
    native_row <- match(sequence[column], amino_acids)
    y <- n_rows - native_row + 1
    rect(column - 0.48, y - 0.48, column + 0.48, y + 0.48,
         border = native_color, lwd = 1.5)
  }
  abline(v = cleavage_after + 0.5, col = cleavage_color, lwd = 2.4)

  text(0.05, n_rows:1, labels = amino_acids, adj = 1, cex = 0.72)
  text(seq_len(n_cols), n_rows + 0.78, labels = sequence, cex = 0.72)
  text(-3.45, n_rows + 2.55, panel_label, font = 2, cex = 1.35, adj = 0)
  text(n_cols / 2 + 0.5, n_rows + 2.42, heading, cex = 0.76)
  text(-3.1, n_rows + 0.78, substrate_label, adj = 0, cex = 0.75)
  text(-2.15, n_rows / 2, "Substitutions", srt = 90, cex = 0.74)

  arrows(cleavage_after + 0.5, n_rows + 2.15,
         cleavage_after + 0.5, n_rows + 1.25,
         length = 0.08, lwd = 1.2, col = cleavage_color)

# 独立色标，范围与论文图一致。
  legend_x <- n_cols + 1.45
  legend_breaks <- seq(score_limits[1], score_limits[2],
                       length.out = length(heat_colors) + 1)
  legend_y <- seq(1, n_rows, length.out = length(legend_breaks))
for (index in seq_along(heat_colors)) {
    rect(legend_x, legend_y[index], legend_x + 1.05, legend_y[index + 1],
         col = heat_colors[index], border = NA)
  }
  rect(legend_x, 1, legend_x + 1.05, n_rows, border = "#66717a", lwd = 0.6)
  tick_values <- -6:1
  tick_y <- 1 + (tick_values - score_limits[1]) /
    diff(score_limits) * (n_rows - 1)
  segments(legend_x + 1.05, tick_y, legend_x + 1.28, tick_y,
           col = "#555555", lwd = 0.6)
  text(legend_x + 1.45, tick_y, labels = tick_values, adj = 0, cex = 0.66)
}

draw_figure <- function() {
  old_par <- par(no.readonly = TRUE)
  on.exit(par(old_par))
  par(mfrow = c(1, 2), mar = c(1.0, 0.8, 1.0, 0.8),
      oma = c(0, 0, 0, 0), family = "sans")
  draw_profile(vamp1_score, vamp1_sequence, "b", "WT BoNT/X cleavage",
"VAMP1:", cleavage_after = 13)
  draw_profile(caspase1_score, caspase1_sequence, "c", "X(PC) cleavage",
"Caspase-1:", cleavage_after = 11)
}

# 5. 保存图片 ---------------------------------------------------------------
draw_figure()