如何使用R从GenBank文件提取指定或全部基因核苷酸序列

你已经完成了GenBank文件读取和基因注释信息提取，后续序列提取可直接基于现有gb和GENES对象实现，无需额外读取参考基因组文件，具体方法如下：

R语言实现方案 #

依赖加载与基础序列获取 #

首先加载序列处理包，从已读入的GenBank对象中提取全基因组序列：

library(Biostrings)
# 从genbankr读取的对象中直接获取全基因组参考序列
genome_seq <- getSeq(gb)

按locus_tag提取单个基因序列 #

注意负链基因需要做反向互补转换，以下函数自动处理链方向问题：

# 定义单基因提取函数
get_gene_by_locus <- function(target_locus, gene_set = GENES, ref_seq = genome_seq) {
  target_interval <- gene_set[gene_set$locus_tag == target_locus]
  if (length(target_interval) == 0) {
    stop("输入的locus_tag无匹配结果，请检查拼写")
  }
  # 自动根据链信息提取正确方向的序列
  res_seq <- getSeq(ref_seq, target_interval)
  names(res_seq) <- target_locus
  return(res_seq)
}

# 示例：提取BJE04_RS00275的序列
single_gene <- get_gene_by_locus("BJE04_RS00275")
# 查看序列内容
single_gene
# 转换为纯字符串格式使用 as.character(single_gene)

批量导出所有基因序列为FASTA格式 #

直接批量提取所有基因区间序列，导出为标准FASTA文件：

# 可选：先过滤假基因，不需要可跳过
GENES_filtered <- GENES[!GENES$pseudo]

# 批量提取所有基因序列，序列名设置为locus_tag
all_gene_seq <- getSeq(genome_seq, GENES_filtered)
names(all_gene_seq) <- GENES_filtered$locus_tag

# 导出到本地文件
writeXStringSet(all_gene_seq, filepath = "YJ016_all_genes.fasta", format = "fasta")

导出的FASTA文件完全符合你给出的示例格式，每条序列的标题为对应locus_tag，下行为核苷酸序列。

备选命令行方案（无R环境可使用） #

如果不想配置R环境，可使用本地安装的seqkit工具，单条命令即可完成提取：

• 批量提取所有基因序列导出FASTA：

seqkit gengb -g gene -n locus_tag YJ016_I.gb > YJ016_all_genes.fasta

• 按locus_tag调取单个基因序列：

seqkit grep -p "BJE04_RS00275" YJ016_all_genes.fasta

注意：所有操作均基于本地文件完成，不需要上传数据到外部服务器

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