求可统计多序列比对中各位置氨基酸频率的工具

R语言实现 #

直接用基础R函数即可完成，无需额外安装包，代码示例如下：

# 输入氨基酸序列（仅保留序列部分）
seqs <- c("MMRSA", "MMLSA", "MMRTA")

# 将每个序列拆分为单个氨基酸的列表
split_seqs <- strsplit(seqs, "")

# 转换为矩阵，每行对应一条序列，每列对应一个位置
seq_matrix <- do.call(rbind, split_seqs)

# 逐列统计氨基酸频率并格式化输出
for (pos in 1:ncol(seq_matrix)) {
  aa_counts <- table(seq_matrix[, pos])
  count_str <- paste(names(aa_counts), "=", aa_counts, collapse = ", ")
  cat(sprintf("Position %d - %s;\n", pos, count_str))
}

运行后会直接输出你期望的位置频率结果。

Linux命令行实现（awk） #

如果习惯用Linux工具，可通过awk脚本处理，假设你的序列存储在align_seqs.txt文件中（内容为你提供的物种+序列格式），执行以下命令：

awk '{
    seq = $NF  # 提取每行最后一段的氨基酸序列
    for (i=1; i<=length(seq); i++) {
        aa = substr(seq, i, 1)
        counts[i][aa]++
    }
}
END {
    for (pos in counts) {
        printf "Position %d - ", pos
        first = 1
        for (aa in counts[pos]) {
            if (!first) printf ", "
            printf "%s = %d", aa, counts[pos][aa]
            first = 0
        }
        print ";"
    }
}' align_seqs.txt

执行后会按位置输出各氨基酸的出现次数。

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