如何从FASTA文件计算距离矩阵？R语言dist.aa报错求助

错误原因分析 #

运行dist.aa()时报错的核心原因是输入格式不匹配：dist.aa()属于ape包，要求输入为AAbin格式对象，但你用Biostrings的readAAStringSet()读取的是AAStringSet类型，直接传入会触发计算逻辑错误；当两条序列完全一致时，格式不兼容的问题会被放大，导致numeric()函数接收到无效长度参数。

解决方案代码 #

1. 正确读取并转换序列格式 #

先将AAStringSet转换为ape包支持的AAbin格式，再计算氨基酸距离：

library(tidyverse)
library(ape)
library(Biostrings)

# 读取FASTA文件
naseq1 <- readAAStringSet("File1.fasta")
naseq2 <- readAAStringSet("File2.fasta")

# 转换为AAbin格式（ape包要求的格式）
aabin1 <- as.AAbin(naseq1)
aabin2 <- as.AAbin(naseq2)

# 重新计算氨基酸距离
dist.aa(aabin1)
dist.aa(aabin2)

运行后，File1会返回两条完全相同序列的距离（0），File2会返回两条序列的差异距离。

2. 计算欧氏距离矩阵 #

先将氨基酸序列编码为数值矩阵，再计算欧氏距离：

# 定义氨基酸转数值的映射（示例：按字母顺序编码，可按需替换为BLOSUM等生化矩阵）
aa_map <- setNames(1:20, c("A","R","N","D","C","Q","E","G","H","I","L","K","M","F","P","S","T","W","Y","V"))

# 把AAStringSet转为数值矩阵的工具函数
seq_to_matrix <- function(aa_set) {
  aa_set %>% 
    as.character() %>% 
    str_split("", simplify = TRUE) %>% 
    apply(2, function(x) aa_map[x])
}

mat1 <- seq_to_matrix(naseq1)
mat2 <- seq_to_matrix(naseq2)

# 计算欧氏距离
euclidean_dist1 <- dist(t(mat1), method = "euclidean")
euclidean_dist2 <- dist(t(mat2), method = "euclidean")

# 转为矩阵格式（可选）
as.matrix(euclidean_dist1)
as.matrix(euclidean_dist2)

3. 计算皮尔逊相关系数矩阵 #

基于上述数值矩阵，直接计算皮尔逊相关系数：

# 计算皮尔逊相关系数
pearson_cor1 <- cor(t(mat1), method = "pearson")
pearson_cor2 <- cor(t(mat2), method = "pearson")

# 输出结果
pearson_cor1
pearson_cor2

关键说明 #

• AAbin是ape包专门用于存储氨基酸序列的格式，和Biostrings的AAStringSet不兼容，必须转换后才能使用dist.aa()。
• 氨基酸转数值的映射可根据需求调整，比如使用BLOSUM矩阵得分，能更精准反映氨基酸的生化特性差异。

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