R中如何将seqinr SeqFastadna对象转为Biostrings DNAStringSet用于多序列比对

实现原理 #

seqinr::read.fasta()返回的是列表格式，每个元素为单个碱基拆分的字符向量，msa要求输入为Biostrings::DNAStringSet类对象，只需在内存中完成格式转换即可，无需导出中间文件。

转换代码 #

在你筛选得到seq_sub之后，插入以下代码即可完成转换并运行多序列比对：

# 转换seqinr列表为DNAStringSet
library(Biostrings)
seq_sub_dna <- DNAStringSet(
  x = sapply(seq_sub, function(seq) paste(toupper(seq), collapse = "")),
  use.names = TRUE
)

# 正常运行msa
seq_alignment <- msa(seq_sub_dna, method="ClustalOmega")

说明 #

• 代码中toupper()是为了适配seqinr默认读入小写碱基的特性，如果你在read.fasta时设置了as.string = TRUE或者已手动处理大小写，可以去掉该函数
• use.names = TRUE会完整保留你在seqinr处理阶段修改过的序列表头，不需要额外赋值

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