在R语言中基于局部比对识别氨基酸替换的技术问询

当然能搞定！结合你已经在用的DECIPHER包，再加个Bioconductor里的Biostrings，就能精准提取目标区域的氨基酸差异，生成你要的结构化表格。我拿你给的示例序列做演示，一步一步来：

1. 准备工作：加载所需包 #

先确保你装了必要的工具包，没装的话先运行这段代码安装：

if (!requireNamespace("BiocManager", quietly = TRUE))
    install.packages("BiocManager")
BiocManager::install(c("DECIPHER", "Biostrings"))

然后加载包：

library(DECIPHER)
library(Biostrings)

2. 读取FASTA序列并完成比对 #

把你给的示例序列存成FASTA文件（比如命名为seqs.fasta），然后读取并执行比对——你已经用过AlignSeqs()，这里直接沿用：

# 读取氨基酸序列
seqs <- readAAStringSet("seqs.fasta")

# 执行序列比对
aligned_seqs <- AlignSeqs(seqs)

3. 映射参考序列的位置（从68开始） #

你的参考序列首氨基酸位置是68，我们先提取参考序列，再构建位置映射规则：

# 提取参考序列（假设FASTA中第一条是ref）
ref_seq <- aligned_seqs["ref"]
# 找出参考序列中非空位的原始位置
ref_valid_pos <- which(strsplit(as.character(ref_seq), "")[[1]] != "-")
# 映射为从68开始的编号（第一个有效位置对应68，所以用67+位置索引）
mapped_pos <- 67 + ref_valid_pos

4. 提取差异并生成目标表格 #

接下来遍历所有查询序列，对比参考序列找出差异，整理成你要的表格结构：

# 初始化结果列表
result_list <- list()

# 遍历所有非参考序列
for (seq_id in names(aligned_seqs)[names(aligned_seqs) != "ref"]) {
    query_seq <- aligned_seqs[seq_id]
    # 拆分序列为单个氨基酸字符
    ref_aa_vec <- strsplit(as.character(ref_seq), "")[[1]]
    query_aa_vec <- strsplit(as.character(query_seq), "")[[1]]
    
    # 筛选出：参考非空位、查询非空位、且氨基酸不同的位置
    diff_idx <- which(ref_aa_vec != "-" & query_aa_vec != "-" & ref_aa_vec != query_aa_vec)
    
    # 对应到映射后的参考位置编号
    diff_pos <- mapped_pos[match(diff_idx, ref_valid_pos)]
    
    # 整理成数据框行
    result_row <- data.frame(
        ID = rep(seq_id, length(diff_pos)),
        Reference_AA = ref_aa_vec[diff_idx],
        Sample_AA = query_aa_vec[diff_idx],
        Pos = diff_pos,
        stringsAsFactors = FALSE
    )
    
    result_list[[seq_id]] <- result_row
}

# 合并所有结果为最终表格
final_table <- do.call(rbind, result_list)

5. 查看结果 #

运行完上面的代码后，打印final_table就能得到你预期的结果：

print(final_table)

输出如下：

ID Reference_AA Sample_AA Pos
query1 query1            R         T  70
query2 query2            A         L  72

关键说明 #

• 这个方法会自动忽略比对中的空位，只聚焦参考序列有氨基酸的目标区域
• 位置映射严格按照你要求的“参考首氨基酸位置为68”设置，确保Pos列完全对应参考序列的编号
• 真实数据中有更多查询序列时，代码可以直接批量处理，不用修改逻辑

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