如何用tidyverse替换基因数据中Upstream与Downstream间的整行

问题需求 #

针对gene_id列中的每个唯一值，完成以下操作：

• 将位于exon_identity列中"Upstream"和"Downstream"之间的整行，替换为同组内exon_identity为"Event"的行
• 移除原"Event"行
• 对于不存在需要替换行的基因（如示例中的"B"），保留原数据结构

测试数据集 #

test_df <- data.frame(
  start = c(2, 9, 13, 19, 13, 20, 25, 35, 39),
  end = c(8, 12, 18, 24, 16, 24, 30, 38, 45),
  gene_id = c("A", "A", "A", "A", "A", "B", "B", "B", "B"),
  exon_identity = c(NA, "Upstream", NA, "Downstream", "Event", NA, "Upstream", "Downstream", NA)
)

初始数据集输出：

> test_df
  start end gene_id exon_identity
1     2   8       A          <NA>
2     9  12       A      Upstream
3    13  18       A          <NA>
4    19  24       A    Downstream
5    13  16       A         Event
6    20  24       B          <NA>
7    25  30       B      Upstream
8    35  38       B    Downstream
9    39  45       B          <NA>

尝试的代码及问题 #

使用以下tidyverse代码时出现长度不匹配警告，输出结果不符合预期：

library(tidyverse)

test_replace <- test_df %>% 
  group_by(gene_id) %>% 
  mutate(start = replace(start, row_number() > which(exon_idnetity == "Upstream") & row_number() < which(exon_idnetity == "Downstream"), start[exon_idnetity == "Event"]),
         end = replace(end, row_number() > which(exon_idnetity == "Upstream") & row_number() < which(exon_idnetity == "Downstream"), end[exon_idnetity == "Event"]),
         exon_idnetity = replace(exon_idnetity, row_number() > which(exon_idnetity == "Upstream") & row_number() < which(exon_idnetity == "Downstream"), "Event")
         )

警告信息：

Warning message:
There were 2 warnings in `mutate()`.
The first warning was:
ℹ In argument: `start = replace(...)`.
ℹ In group 1: `gene_id = "A"`.
Caused by warning in `x[list] <- values`:
! number of items to replace is not a multiple of replacement length
ℹ Run dplyr::last_dplyr_warnings() to see the 1 remaining warning. 

错误输出：

> test_replace
# A tibble: 9 × 4
# Groups:   gene_id [2]
  start   end gene_id exon_idnetity
  <dbl> <dbl> <chr>   <chr>        
1     2     8 A       NA           
2     9    12 A       Upstream     
3    NA    NA A       Event        
4    19    24 A       Downstream   
5    13    16 A       Event        
6    20    24 B       NA           
7    25    30 B       Upstream     
8    35    38 B       Downstream   
9    39    45 B       NA     

期望输出 #

> desired_outcome 
  start end gene_id exon_idnetity
1     2   8       A          <NA>
2     9  12       A      Upstream
3    13  16       A         Event
4    19  24       A    Downstream
5    20  24       B          <NA>
6    25  30       B      Upstream
7    35  38       B    Downstream
8    39  45       B          <NA>

正确解决方案 #

基于tidyverse的实现代码如下：

library(tidyverse)

test_df %>%
  group_by(gene_id) %>%
  mutate(
    # 定位Upstream和Downstream的行号
    up_pos = which(exon_identity == "Upstream"),
    down_pos = which(exon_identity == "Downstream"),
    # 标记需要替换的行（Upstream与Downstream之间的行）
    to_replace = row_number() > up_pos & row_number() < down_pos,
    # 提取Event行的start和end值（无Event行则为NA）
    event_start = ifelse(any(exon_identity == "Event"), start[exon_identity == "Event"], NA),
    event_end = ifelse(any(exon_identity == "Event"), end[exon_identity == "Event"], NA)
  ) %>%
  # 替换目标行的字段值
  mutate(
    start = ifelse(to_replace, event_start, start),
    end = ifelse(to_replace, event_end, end),
    exon_identity = ifelse(to_replace, "Event", exon_identity)
  ) %>%
  # 过滤掉原Event行（仅保留替换后的Event行）
  filter(!(exon_identity == "Event" & !to_replace)) %>%
  # 移除临时辅助列并取消分组
  select(-up_pos, -down_pos, -to_replace, -event_start, -event_end) %>%
  ungroup()

代码逻辑说明 #

1. 分组定位：按gene_id分组后，找到每组中"Upstream"和"Downstream"的行位置，标记出两者之间需要替换的行
2. 提取替换值：获取同组内"Event"行的start和end数值
3. 字段替换：对标记的目标行，替换其start、end和exon_identity字段为Event行对应值
4. 过滤清理：移除原始的"Event"行，删除临时辅助列并取消分组

运行上述代码后，即可得到符合预期的输出结果。

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