含间隔的DNA序列对应氨基酸起始/终止位置自动计算方法

问题背景 #

现有带间隔的DNA序列myseq，每3个有效碱基对应氨基酸序列myaa中的一个氨基酸（共13个）。无间隔时可以用seq()快速生成每个碱基组的起始(st_pos)和终止(en_pos)位置，进而构建包含位置与对应氨基酸的mydf数据框，但序列存在间隔时无法直接依赖seq()，需要自动生成位置向量的方法。

间隔案例示例 #

案例1 #

• 带间隔DNA序列：myseq <- "ATG-TCG--AGCCTG-TAGAAGCTT-GCT"
• 对应氨基酸序列：myaa <- c("M", "S", "A", "L", "*", "K", "L", "A", "X", "Y", "P", "Q", "R")（共13个）
• 期望生成的mydf核心结构：

  amino_acid	st_pos	en_pos
  M	1	3
  S	5	7
  A	10	12
  L	13	15
  *	17	19
  ...	...	...

案例2 #

• 带间隔DNA序列：myseq <- "---ATGTCGAGC--CTGTAG---AAG--CTTGCT"
• 对应氨基酸序列同案例1
• 期望生成的mydf核心结构：

  amino_acid	st_pos	en_pos
  M	4	6
  S	7	9
  A	10	12
  L	15	17
  *	18	20
  ...	...	...

自动生成位置向量的方法 #

用R语言实现的自定义函数，可适配任意间隔符号（以下以-作为间隔为例，可按需修改）：

generate_pos_df <- function(myseq, myaa) {
  # 将DNA序列拆分为单个字符
  seq_chars <- strsplit(myseq, "")[[1]]
  # 筛选出非间隔碱基的原始位置
  valid_base_pos <- which(seq_chars != "-")
  # 每3个有效碱基为一组，对应一个氨基酸
  base_groups <- split(valid_base_pos, ceiling(seq_along(valid_base_pos)/3))
  # 提取每组的起始和终止位置
  st_pos <- sapply(base_groups, min)
  en_pos <- sapply(base_groups, max)
  # 构建最终数据框
  data.frame(
    amino_acid = myaa,
    st_pos = st_pos,
    en_pos = en_pos,
    stringsAsFactors = FALSE
  )
}

函数说明 #

1. 先拆分序列为单个字符，定位所有非间隔碱基的原始位置；
2. 按每3个有效碱基为一组，对应myaa中的一个氨基酸；
3. 每组的最小位置即为该氨基酸对应的碱基起始位置，最大位置为终止位置；
4. 若间隔符号不是-，只需修改seq_chars != "-"中的符号即可，比如换成.或其他标记。

测试示例 #

# 案例1测试
myseq1 <- "ATG-TCG--AGCCTG-TAGAAGCTT-GCT"
myaa1 <- c("M", "S", "A", "L", "*", "K", "L", "A", "X", "Y", "P", "Q", "R")
mydf1 <- generate_pos_df(myseq1, myaa1)

# 案例2测试
myseq2 <- "---ATGTCGAGC--CTGTAG---AAG--CTTGCT"
mydf2 <- generate_pos_df(myseq2, myaa1)

运行后即可得到与手动编码一致的mydf数据框。

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