在R中从多物种FASTA文件按物种生成指定规则的共识序列

问题背景 #

处理一个已比对完成的FASTA文件（示例：Andrenidae.FASTA），文件包含十几种物种的数百条长度一致的核苷酸序列，序列由A/C/T/G/N组成（N代表未知核苷酸），序列名称以物种名开头。

需要为每个物种生成符合以下规则的共识序列：

• 每个位点选取出现频率最高的核苷酸
• A/C/T/G优先级高于N：仅当该物种所有序列的该位点均为N时，共识序列该位才取N；若存在至少一条序列该位为A/C/T/G，则从这些碱基中选择频率最高的

当前使用的代码将A/C/T/G/N同等处理，导致共识序列中出现大量N，原代码如下：

library(Biostrings) 

seqs <- readDNAStringSet("Andrenidae.FASTA")
species_names <- sapply(names(seqs), function(x) strsplit(x, " ")[[1]][1])
species_sequences <- split(seqs, species_names)

#get the species names in a vector
sp<- unique(species_names)
#create empty list to fill with consensus sequences
con_seq <- list()

#calculate consensus sequences
for(i in 1:10){
  con_seq[[i]] <- consensusString(species_sequences[[i]], ambiguityMap="N", 
                                  threshold=0.0001)
}

#unlist con seq into one large vector
cs_all<-unlist(con_seq)
#create a dataframe with the species names and their corresponding con seq
cs_all_df <- as.data.frame(cbind(sp, cs_all))
#write out df
write.csv(cs_all_df, file='Andrenidae_con.csv')

解决方案 #

通过Biostrings包的consensusMatrix统计每个位点的碱基计数，先过滤N的影响，再按规则生成共识序列，修改后的代码如下：

library(Biostrings)

# 读取FASTA文件
seqs <- readDNAStringSet("Andrenidae.FASTA")

# 提取物种名（假设序列名以物种名开头，空格分隔）
species_names <- sapply(names(seqs), function(x) strsplit(x, " ")[[1]][1])
species_sequences <- split(seqs, species_names)

# 获取所有物种名
sp <- unique(species_names)
con_seq <- list()

# 为每个物种生成共识序列
for (i in seq_along(sp)) {
  # 获取当前物种的所有序列
  current_seqs <- species_sequences[[sp[i]]]
  # 生成碱基计数矩阵：行=碱基(A/C/G/T/N)，列=位点
  count_matrix <- consensusMatrix(current_seqs)
  
  # 遍历每个位点确定共识碱基
  consensus_bases <- character(ncol(count_matrix))
  for (pos in 1:ncol(count_matrix)) {
    # 提取A/C/G/T的计数，排除N
    acgt_counts <- count_matrix[c("A", "C", "G", "T"), pos]
    total_acgt <- sum(acgt_counts)
    
    if (total_acgt == 0) {
      # 所有序列该位点都是N，共识位取N
      consensus_bases[pos] <- "N"
    } else {
      # 找到A/C/G/T中计数最高的碱基（若多个碱基计数相同，取第一个出现的）
      max_count <- max(acgt_counts)
      consensus_bases[pos] <- names(acgt_counts)[acgt_counts == max_count][1]
    }
  }
  
  # 拼接成完整的共识序列
  con_seq[[i]] <- paste(consensus_bases, collapse = "")
}

# 整理成数据框并导出
cs_all_df <- data.frame(
  species = sp, 
  consensus_sequence = unlist(con_seq), 
  stringsAsFactors = FALSE
)
write.csv(cs_all_df, file = 'Andrenidae_con.csv', row.names = FALSE)

关键说明 #

1. 使用consensusMatrix生成每个位点的碱基计数，清晰统计每个碱基的出现次数
2. 对每个位点优先判断是否存在非N碱基：
   • 若全为N，则共识位保留N
   • 若存在非N碱基，仅在A/C/G/T中选取频率最高的碱基，彻底排除N的干扰
3. 遍历所有物种（原代码只处理了前10个物种，修改后覆盖全部），确保每个物种都生成对应的共识序列

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