如何在R中将CLUSTAL格式蛋白比对结果转为指定表格？

步骤1：清洗比对序列 #

首先清理clustal.res中序列包含的数字、制表符等干扰字符，得到纯净的比对序列（包含空位-）：

# 加载seqinr包
library(seqinr)

# 清洗序列：移除数字、制表符和换行符
clean_seq <- lapply(clustal.res$seq, function(x) {
  gsub("[0-9\t\n]", "", x)
})

# 校验两个序列长度一致（比对后序列长度必须相同）
stopifnot(nchar(clean_seq[[1]]) == nchar(clean_seq[[2]]))

步骤2：生成残基位置向量 #

为两个序列分别计算有效残基的位置（跳过空位-，仅对实际氨基酸计数）：

# 生成SMARCE1的位置向量，空位位置设为NA
smarce1_pos <- cumsum(clean_seq[[1]] != "-")
smarce1_pos[clean_seq[[1]] == "-"] <- NA

# 生成swsn-3的位置向量，空位位置设为NA
swsn3_pos <- cumsum(clean_seq[[2]] != "-")
swsn3_pos[clean_seq[[2]] == "-"] <- NA

步骤3：提取比对共识符号 #

直接从原始CLUSTAL文件提取共识符号行（避免read.alignment丢失信息）：

# 读取CLUSTAL文件文本
clustal_text <- readLines("sample.aln")
# 筛选出仅包含.:*的共识行
consensus_line <- grep("^\\s*[.:*]+", clustal_text, value = TRUE)
# 清理无关字符，得到纯净的符号序列
consensus_symbols <- gsub("[^.:*]", "", consensus_line)

# 校验符号序列长度与比对序列一致
stopifnot(nchar(consensus_symbols) == nchar(clean_seq[[1]]))

# 将无符号的位置设为"N/A"
consensus_vec <- strsplit(consensus_symbols, "")[[1]]
consensus_vec[consensus_vec == ""] <- "N/A"

步骤4：合并生成目标表格 #

拆分残基、替换空位为N/A，最终合并为数据框：

# 拆分序列为单个残基向量
smarce1_residues <- strsplit(clean_seq[[1]], "")[[1]]
swsn3_residues <- strsplit(clean_seq[[2]], "")[[1]]

# 将空位替换为"N/A"
smarce1_residues[smarce1_residues == "-"] <- "N/A"
swsn3_residues[swsn3_residues == "-"] <- "N/A"

# 生成最终表格
result_table <- data.frame(
  "Row" = 1:length(smarce1_residues),
  "SMARCE1位置" = smarce1_pos,
  "SMARCE1残基" = smarce1_residues,
  "swsn-3_protein位置" = swsn3_pos,
  "swsn-3_protein残基" = swsn3_residues,
  "比对符号" = consensus_vec,
  stringsAsFactors = FALSE
)

# 查看前几行结果
head(result_table)

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