在线性区间中查找局部坐标：基于基因转录本区间坐标数据框

嘿，我来帮你搞定这个转录本线性坐标查找的问题！看起来你需要实现基因组坐标与转录本拼接后的线性局部坐标之间的转换——不管是把基因组上的位置转成转录本里的局部位置，还是反过来，下面是基于你的data.frame的具体实现方案：

首先你的原始数据里没有转录本ID，这是分组的关键。我先假设可以按数据的结构给每个转录本分配ID（比如每2行对应一个转录本，你可以根据实际数据调整规则），然后把每个转录本的外显子按基因组位置排序：

# 给每个转录本分配唯一ID（示例规则：每2行一组）
df$transcript_id <- rep(paste0("trans_", 1:(nrow(df)/2)), each=2)

# 用dplyr按转录本分组，然后按基因组start位置排序外显子
library(dplyr)
df <- df %>% 
  group_by(transcript_id) %>% 
  arrange(start, .by_group = TRUE) %>% 
  ungroup()

接下来，我们要给每个外显子计算它在转录本拼接后的线性序列里的起始和结束位置：

df <- df %>%
  group_by(transcript_id) %>%
  # 先算每个外显子的长度（注意基因组坐标是闭区间，所以要+1）
  mutate(exon_length = end - start + 1) %>%
  # 计算该外显子在转录本中的起始位置：前面所有外显子长度之和 +1
  mutate(transcript_start = cumsum(lag(exon_length, default = 0)) + 1) %>%
  # 计算该外显子在转录本中的结束位置
  mutate(transcript_end = transcript_start + exon_length - 1) %>%
  ungroup()

现在你的df里就多了transcript_start和transcript_end两列，代表这个外显子在转录本线性序列中的位置区间。

写一个简单的函数，输入基因组的染色体、位置和转录本ID，就能返回对应的局部坐标：

genome_to_local <- function(df, seqname, pos, transcript_id) {
  # 筛选出目标转录本中包含该基因组位置的外显子
  target_exon <- df %>%
    filter(seqnames == seqname, transcript_id == transcript_id, start <= pos, end >= pos)
  
  if(nrow(target_exon) == 0) {
    stop("这个坐标不在目标转录本的任何外显子里哦！")
  }
  
  # 计算局部坐标：外显子内的偏移量 + 转录本中的起始偏移
  local_pos <- (pos - target_exon$start + 1) + (target_exon$transcript_start - 1)
  return(local_pos)
}

# 举个例子：查chr1上3670600在trans_1里的局部坐标
genome_to_local(df, "chr1", 3670600, "trans_1")

反过来，如果你有转录本的局部坐标，想转成基因组坐标，用这个函数：

local_to_genome <- function(df, transcript_id, local_pos) {
  # 筛选出目标转录本中包含该局部坐标的外显子
  target_exon <- df %>%
    filter(transcript_id == transcript_id, transcript_start <= local_pos, transcript_end >= local_pos)
  
  if(nrow(target_exon) == 0) {
    stop("这个局部坐标不在目标转录本的线性区间里哦！")
  }
  
  # 计算基因组坐标：外显子起始位置 + 局部坐标在该外显子内的偏移
  genome_pos <- target_exon$start + (local_pos - target_exon$transcript_start)
  return(list(seqname = target_exon$seqnames, pos = genome_pos))
}

# 举个例子：查trans_1里局部坐标100对应的基因组坐标
local_to_genome(df, "trans_1", 100)

如果你的数据里同时有完整外显子和编码序列（CDS）的区间，只需要给数据加个type标记（比如区分是exon还是cds），然后单独对CDS分组计算线性坐标就行：

# 假设每两行里第一行是完整外显子，第二行是CDS区间（根据你的数据结构调整）
df$type <- rep(c("exon", "cds"), nrow(df)/2)

# 单独处理CDS的线性坐标
cds_df <- df %>%
  filter(type == "cds") %>%
  group_by(transcript_id) %>%
  arrange(start, .by_group = TRUE) %>%
  mutate(cds_length = end - start + 1) %>%
  mutate(cds_start = cumsum(lag(cds_length, default = 0)) + 1) %>%
  mutate(cds_end = cds_start + cds_length - 1) %>%
  ungroup()

之后就可以用类似上面的函数处理CDS的坐标转换啦！

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