识别两条等长序列突变位置及类型的优雅实现方案（R/Python/bash）

R 实现（适配Biostrings使用场景） #

无需显式循环，采用R原生矢量化操作实现，适配Biostrings读取的序列对象，性能可适配长序列场景：

library(Biostrings)
# 读取fasta序列文件
seqs <- readAAStringSet("PSE-1_Round20.fas")
# 拆分参考序列与待比较序列为单字符向量
ref_char <- strsplit(as.character(seqs[1]), "")[[1]]
query_char <- strsplit(as.character(seqs[2]), "")[[1]]
# 筛选差异位置并拼接为要求格式
diff_pos <- which(ref_char != query_char)
mutation_res <- paste0(ref_char[diff_pos], diff_pos, query_char[diff_pos], collapse = ", ")

如果需要处理多序列比对场景，可直接调用Biostrings内置的compareStrings或pairwiseAlignment函数提取差异位点，无需手动拆分字符。

Python 实现 #

采用列表推导式实现，逻辑简洁易集成到自定义分析流程：

def fetch_mutation_info(ref_seq: str, query_seq: str, start_idx: int = 1) -> str:
    """
    识别两条等长序列的突变信息
    start_idx: 序列位置计数起始值，生物序列通常从1开始计数
    """
    mut_list = [
        f"{r}{i + start_idx}{q}" 
        for i, (r, q) in enumerate(zip(ref_seq, query_seq)) 
        if r != q
    ]
    return ", ".join(mut_list)

# 示例调用
print(fetch_mutation_info("AGGCTAC", "AGCCTTC"))
# 输出结果：G3C, A6T

Bash 实现 #

适合命令行批量处理场景，无需依赖额外编程环境：

# 假设参考序列存放在ref.fasta，待比较序列存放在query.fasta，自动跳过fasta注释行
awk '
BEGIN{
    # 读取参考序列
    while((getline < "ref.fasta") > 0) {
        if(substr($0,1,1) != ">") ref = ref $0
    }
    # 读取待比较序列
    while((getline < "query.fasta") > 0) {
        if(substr($0,1,1) != ">") query = query $0
    }
    len = length(ref)
}
END{
    for(i=1; i<=len; i++) {
        r = substr(ref, i, 1)
        q = substr(query, i, 1)
        if(r != q) res = res (res?", ":"") r i q
    }
    print res
}'

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