如何用tidyverse或等效工具对带末端Ns的测序数据分组计数？

我不确定tidyverse是否支持模糊搜索或分组功能，希望将测序数据合并分组。原因是测序仪可能在序列末端生成人工Ns，但中间核心序列完全一致，需找到高效方法筛选这些序列并分组计数。

例如，序列ATCGACACACACACACACACACAATC、NTCGACACACACACACACACACAATC和NNCGACACACACACACACACACAATN应归为同一组，共用核心序列CGACACACACACACACACACAAT作为组名。

约束条件 #

• Ns仅允许出现在序列两端（一端或两端），若Ns位于序列中间，则需归为单独分组。
• 仅涉及标准DNA（或氨基酸）编码：DNA使用ATCG，N代表随机碱基；氨基酸使用ARNDBCEQZGHILKMFPSTWYV，X代表随机残基。

示例数据集及预期结果 #

输入数据集 #

data <- data.frame(DNA = c("ATCGACACACACACACACACACAATC", "NTCGACACACACACACACACACAATC", "NNCGACACACACACACACACACAATN", "NNCGACACACACACACACACACAATN", "NNCGACACACACACACACACACAATN", "NNCGACACACACACACACACACAATN", "ATCGACACACACACACACACACAATC", "NTCGACACACACACACACACACAATC", "ATCGACACACACACACACACACAATC", "NTCGACACACACACACACACACAATC", "NTCGACACACACACACACACACAATC", "AGAGTCTCGATCG", "TCTCTGATGCTAAA", "TAGCTAGACTAGCATCGACTACGACT", "TAGCTAGACNNGCATCGACTACGACT", "TAGCTAGACNAGCATCGACTACGACT", "TAGCTAGACTAGCATCGACTACGACT", "TAGCTAGACTAGCATCGACTACGACT", "NAGCTAGACTAGCATCGACTACGACT", "NNCTAGCATCGACTACGACT", "NNCTAGCATCGACTACNACT", "GTCGACACACACACACACACACAATC"))

预期输出结果 #

grouped_data <- data.frame(DNA_group = c("CGACACACACACACACACACAAT", # 如描述中的核心序列
                                         "AGAGTCTCGATCG", # 序列长度可能不同，理想情况下代码允许定义最小匹配长度，当前设为10，实际会超过100
                                         "TCTCTGATGCTAAA", "GTCGACACACACACACACACACAATC",
                                         "CTAGCATCGACTACGACT", # 同第一个示例
                                         "TAGCTAGACNNGCATCGACTACGACT", # N位于中间，单独成组
                                         "TAGCTAGACNAGCATCGACTACGACT", # N位于中间，单独成组
                                         "CTAGCATCGACTACNACT"), # N位于中间，即使两端有N也单独成组，理想情况可修剪末端Ns（可选）
                           count = c(11, 1, 1, 1, 5, 1, 1, 1))

已尝试方法 #

曾尝试用grep或gsub创建临时列并计数，但不确定能否实现动态grep（因为组名不同，且不清楚每组最长的代表性序列）。偏好使用R语言处理，也可考虑其他高效的开源工具。

核心思路 #

1. 先判断序列是否存在中间的N/X：通过正则表达式匹配[ATCG]+N[ATCG]+（DNA）或[ARNDBCEQZGHILKMFPSTWYV]+X[ARNDBCEQZGHILKMFPSTWYV]+（氨基酸），存在则保留原序列作为组名。
2. 对于无中间N/X的序列，修剪两端的N/X，得到核心序列作为组名。
3. 最后按组名分组计数。

R语言实现代码（基于tidyverse） #

library(tidyverse)

# 定义序列处理函数
process_sequence <- function(seq, type = "DNA") {
  # 定义对应类型的有效字符和模糊字符
  if (type == "DNA") {
    valid_chars <- "ATCG"
    fuzzy_char <- "N"
  } else if (type == "AA") {
    valid_chars <- "ARNDBCEQZGHILKMFPSTWYV"
    fuzzy_char <- "X"
  } else {
    stop("Type must be 'DNA' or 'AA'")
  }
  
  # 匹配中间存在模糊字符的正则表达式
  middle_fuzzy_pattern <- str_glue("[{valid_chars}]+{fuzzy_char}[{valid_chars}]+")
  
  # 检查是否有中间模糊字符
  if (str_detect(seq, middle_fuzzy_pattern)) {
    # 可选：修剪两端的模糊字符（按需启用）
    # seq <- str_remove_all(seq, str_glue("^{fuzzy_char}+|{fuzzy_char}+$"))
    return(seq)
  } else {
    # 修剪两端模糊字符得到核心序列
    core_seq <- str_remove_all(seq, str_glue("^{fuzzy_char}+|{fuzzy_char}+$"))
    # 可选：添加最小核心序列长度过滤（按需启用）
    # if (nchar(core_seq) >= 10) return(core_seq) else return(seq)
    return(core_seq)
  }
}

# 处理数据并分组计数
grouped_result <- data %>%
  mutate(DNA_group = map_chr(DNA, process_sequence, type = "DNA")) %>%
  group_by(DNA_group) %>%
  summarise(count = n(), .groups = "drop") %>%
  arrange(desc(count)) # 可选：按计数降序排列

# 查看结果
print(grouped_result)

代码说明 #

• process_sequence函数：根据序列类型自动适配规则，先识别中间模糊字符，再决定是保留原序列还是提取核心序列。
• 可选功能：取消注释即可启用末端模糊字符修剪、核心序列最小长度过滤。
• 效率：基于tidyverse的向量化操作，处理大规模测序数据性能稳定；若需进一步提速，可替换为data.table实现。

其他工具推荐 #

针对超大规模数据（千万级以上序列），可考虑：

• seqkit：开源命令行工具，支持快速修剪、分组计数，适合批量处理FASTA/FASTQ文件。
• vsearch：专业序列聚类工具，可通过调整相似度阈值适配末端N的分组场景。

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