如何基于digest和gene分组移除被包含短序列（tidyverse方案）

现有如下R语言数据框w：

> w
   digest    gene          seq
1     InS  AB0583          AAB
2     InS  AB0583        AABKR
3     InS  AB0583      GFHGHGG
4     PAC PU83022          EUT
5     PAC PU83022      HSFSFJF
6     PAC PU83022        EUTCK
7     PAC PU83022       EUTCKJ
8     InS PO93853         HDGJ
9     InS PO93853        HDGJU
10    InS PO93853       YTYEYD
11    InS PO93853 YTYEYDJHSGSG
12    InS PO93853   SALGHAGGEE

数据构造代码：

w <- data.frame(
  digest = c(rep("InS", 3), rep("PAC", 4), rep("InS", 5)),
  gene = c(rep("AB0583", 3), rep("PU83022", 4), rep("PO93853", 5)),
  seq = c("AAB", "AABKR", "GFHGHGG",
          "EUT", "HSFSFJF", "EUTCK", "EUTCKJ",
          "HDGJ", "HDGJU", "YTYEYD", "YTYEYDJHSGSG", "SALGHAGGEE")
)

需求说明 #

按digest和gene分组，移除同组内被其他seq包含的短序列，仅保留无法被其他序列包含的最长序列，使用tidyverse工具包实现，具体步骤：

• 新增count字段统计seq字符数，按digest、gene、count排序
• 分组后过滤掉被包含的短序列

期望输出 #

> w
   digest    gene          seq count
1     InS  AB0583        AABKR     5 
2     InS  AB0583      GFHGHGG     7
3     InS PO93853        HDGJU     5 
4     InS PO93853   SALGHAGGEE    10
5     InS PO93853 YTYEYDJHSGSG    12 
6     PAC PU83022       EUTCKJ     6 
7     PAC PU83022      HSFSFJF     7

解决方案代码 #

library(tidyverse)

w_processed <- w %>%
  # 新增序列长度字段并排序
  mutate(count = str_length(seq)) %>%
  arrange(digest, gene, count) %>%
  # 分组过滤被包含的短序列
  group_by(digest, gene) %>%
  filter(!map_lgl(seq, ~ any(str_detect(seq[count > str_length(.)], .)))) %>%
  ungroup() %>%
  # 按期望格式排序输出
  arrange(digest, gene, count)

print(w_processed)

代码解释 #

1. mutate(count = str_length(seq))：用stringr工具的str_length函数生成序列长度字段count
2. arrange(digest, gene, count)：按分组字段和序列长度升序排序，便于后续判断包含关系
3. group_by(digest, gene)：按指定字段分组处理
4. filter(!map_lgl(seq, ~ any(str_detect(seq[count > str_length(.)], .))))：对每个序列，检查同组内所有更长的序列是否包含它，若被包含则过滤该短序列
5. ungroup()：取消分组，恢复普通数据框结构
6. 最后再次排序，保证输出格式与期望一致

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