从具有 R 中许多物种序列的 fasta 中获取每个物种的共识序列

Question

在 R 中，我有一个 fasta（例如：Andrenidae.FASTA），其中包含来自十几个物种的数百个核苷酸序列。我想为每个物种获得 1 个共识序列。每条序列由A/C/T/G/N组成，N代表未知核苷酸，序列名称为物种名称。所有序列的长度完全相同，我已经对齐了它们。

问题是我希望在每个位置选择最常见的核苷酸作为共识序列，并且我希望 A/C/T/G 中的任何一个覆盖 N。因此，如果每个序列都应该只选择 N 作为共识因为该物种在该位置包含 N。如果即使有 1 个序列有 A/C/T/G，无论哪一个最常见都会进入共识而不是 N。这是因为这些比对序列来自该基因不同部分的几个较长序列，因此通常是基因的某些部分该基因仅对这些序列中的一小部分进行了测序。

感谢您的帮助！

目前，这个 R 代码平等对待 A/C/T/G/N，所以我最终在我的共识序列中有太多的 N：

library(Biostrings)
rm(list=ls())

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')

Answer 1

我认为您需要从函数中删除

ambiguityMap

和

thresholdd

。 设置

ambiguityMap = N

和

thresholdd = .0001

意味着包含在输入集中出现频率低于 0.01% 的核苷酸的序列将在共识计算中替换为

N

。

这里有一个小例子：

# 3 Ns and one A at 3rd pos
consensusString(DNAStringSet(c("ACNG","ACNR", "ACNG", "ACAG")))
# output:
# "ACAG"

# all Ns 
consensusString(DNAStringSet(c("ACNG","ACNR", "ACNG", "ACNG")))
# output:
#"ACNG"

# with those parameters set
consensusString(DNAStringSet(c("ACNG","ACNR", "ACNG", "ACAG")), ambiguityMap = "N", threshold = 0.001)
# output
#  "ACNN"