如何使用R筛选与TIMEPOINT_1相似度≥95%的TIMEPOINT_2 DNA序列？

核心思路 #

通过计算序列间的莱文斯坦编辑距离（衡量两个字符串的差异程度），将其转换为相似度指标，筛选出TIMEPOINT_2中与TIMEPOINT_1任意序列相似度≥95%的条目。

实现步骤与代码 #

1. 安装并加载依赖包

# 安装（首次运行时）并加载stringdist包
if (!require(stringdist)) {
  install.packages("stringdist")
  library(stringdist)
}

2. 定义相似度计算函数
   该函数以两条DNA序列为输入，基于最长序列长度计算相似度：

calc_similarity <- function(seq_a, seq_b) {
  max_length <- max(nchar(seq_a), nchar(seq_b))
  edit_distance <- stringdist(seq_a, seq_b, method = "lv")
  # 相似度 = (最长长度 - 编辑距离) / 最长长度
  (max_length - edit_distance) / max_length
}

3. 筛选符合条件的晚期序列
   遍历TIMEPOINT_2的每条序列，检查是否存在TIMEPOINT_1中满足相似度阈值的序列：

# 筛选TIMEPOINT_2中符合条件的序列
filtered_time2 <- TIMEPOINT_2[sapply(TIMEPOINT_2$sequence, function(target_seq) {
  # 检查与TIMEPOINT_1中任意序列的相似度是否≥95%
  any(sapply(TIMEPOINT_1$sequence, function(ref_seq) {
    calc_similarity(target_seq, ref_seq) >= 0.95
  }))
}), ]

# 查看筛选结果
print(filtered_time2)

示例结果解释 #

基于你提供的测试数据，筛选后的结果会包含：

• time2_seq1：与time1_seq1的相似度约95.45%
• time2_seq3：与time1_seq3完全匹配（相似度100%）
• time2_seq4：与time1_seq4的相似度约95.83%

常见问题说明 #

若之前使用stringdistmatrix未得到预期结果，大概率是未将编辑距离转换为相似度，或未处理序列长度不一致的情况。上述方法通过基于最长序列长度计算相似度，解决了不同长度序列的比较问题。

内容的提问来源于stack exchange，提问作者Gerald Vasquez Aleman