如何用R语言Biostrings提取比对中长度≥10nt的完全匹配片段及位置

需求说明 #

使用Biostrings完成两个DNA序列的全局比对后，需要提取所有长度至少为10nt的完全匹配序列片段，同时获取这些片段在原始pattern和subject序列中的起始、结束位置。实际处理的序列长度超过15000nt，需要高效解决方案，偏好使用tidyverse工具。

可复现代码示例 #

library('Biostrings')
set.seed(123)
dna_alphabet <- DNA_ALPHABET[1:4]
seq1_random <- sample(dna_alphabet, 100, replace = TRUE)
seq2_random <- sample(dna_alphabet, 100, replace = TRUE)

match_12 <- sample(dna_alphabet, 12, replace = TRUE)
match_15 <- sample(dna_alphabet, 15, replace = TRUE)
match_23 <- sample(dna_alphabet, 23, replace = TRUE)

seq1_random[15:(15+11)] <- match_12
seq2_random[17:(17+11)] <- match_12 
seq1_random[41:(41+14)] <- match_15
seq2_random[40:(40+14)] <- match_15
seq1_random[70:(70+22)] <- match_23
seq2_random[73:(73+22)] <- match_23

seq1 <- DNAString(paste(seq1_random, collapse = ""))
seq2 <- DNAString(paste(seq2_random, collapse = ""))

alignment <- pairwiseAlignment(pattern = seq1, subject = seq2, type = "global")

alignment

运行后比对输出：

Global PairwiseAlignmentsSingleSubject (1 of 1)
pattern: GGGCGCCCGATCCA--CCTATTCGTTTGTCGCACGTCAGGATATTTCACGTGAGTCATCACAAT----TGACAAGACGACCATTTCTAACATGAATTCACTAACGG
subject: ATCATCAGGTTCTGACCCTATTCGTTTGA---ATGCCCTCCCATTTCACGTGAGTCAAGGGAGGCCCGTGGCACTACGACCATTTCTAACATGAATTCGGTAT---
score: -138.1462 

预期结果 #

> df_results
            perfect_match start_pos_pattern end_pos_pattern start_pos_subject end_pos_subject
1            CCTATTCGTTTG                15              26                17              28
2         ATTTCACGTGAGTCA                41              55                40              54
3 ACGACCATTTCTAACATGAATTC                70              92                73              95

当前尝试的问题 #

已提取比对后的序列，但用滑动窗口方法存在以下不足：

• 仅得到比对后pattern序列的位置，未考虑缺口（-）影响，与原始序列位置不符；
• 未获取subject序列中的匹配位置；
• 结果以列表存储，未整理为数据框。

高效解决方案（tidyverse风格） #

使用向量操作和tidyverse的分组统计，高效处理长序列，步骤如下：

library(tidyverse)

# 提取比对后的pattern和subject序列
aligned_pat <- as.character(pattern(alignment))
aligned_subj <- as.character(subject(alignment))

# 生成比对位置到原始序列位置的映射（跳过缺口）
pat_pos_map <- cumsum(str_split(aligned_pat, "", simplify = TRUE) != "-")
subj_pos_map <- cumsum(str_split(aligned_subj, "", simplify = TRUE) != "-")

# 构建包含所有比对位置信息的数据框
alignment_df <- tibble(
  align_pos = seq_along(pat_pos_map),
  pat_char = str_split(aligned_pat, "", simplify = TRUE),
  subj_char = str_split(aligned_subj, "", simplify = TRUE),
  pat_original_pos = pat_pos_map,
  subj_original_pos = subj_pos_map
)

# 筛选完全匹配的位置（无缺口且字符相等），并分组连续匹配片段
matches_df <- alignment_df %>%
  filter(pat_char != "-", subj_char != "-", pat_char == subj_char) %>%
  mutate(group = cumsum(c(TRUE, diff(align_pos) != 1)))

# 统计每个匹配片段的信息，筛选长度≥10的片段
result_df <- matches_df %>%
  group_by(group) %>%
  summarise(
    perfect_match = str_c(pat_char, collapse = ""),
    start_pos_pattern = first(pat_original_pos),
    end_pos_pattern = last(pat_original_pos),
    start_pos_subject = first(subj_original_pos),
    end_pos_subject = last(subj_original_pos),
    .groups = "drop"
  ) %>%
  filter(str_length(perfect_match) >= 10)

# 查看最终结果
result_df

方法说明 #

1. 位置映射：通过cumsum跳过缺口字符，计算每个比对位置对应的原始序列坐标，解决缺口导致的位置偏移问题；
2. 连续片段识别：利用diff(align_pos) != 1判断连续匹配的位置，通过cumsum给连续片段分组；
3. 高效统计：使用tidyverse的分组汇总，一次性提取每个片段的序列内容和原始位置，适合处理长序列。

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