如何在R中匹配字符串模式并提取GenBank(gbk)文件的CDS信息

我们可以使用Bioconductor生态下的genbankr包完成该需求，它原生支持GenBank格式解析，无需手动编写正则匹配规则，可高效处理数千条级别的CDS条目。 #

1. 环境准备 #

首先安装依赖包：

# 若未安装BiocManager先执行安装
if (!require("BiocManager", quietly = TRUE))
    install.packages("BiocManager")
# 安装genbankr
BiocManager::install("genbankr")
# 加载依赖包
library(genbankr)
library(stringr)

2. 单gbk文件处理逻辑 #

核心处理步骤： #

• 读取gbk文件
• 提取所有CDS特征
• 提取目标字段，处理跨行列、缺失值问题
• 按要求格式化输出

# 替换为你的gbk文件路径
gbk_path <- "your_file.gbk"
# 读取gbk文件，ret.seqs设为FALSE可跳过核酸序列读取提升速度
gb <- readGenBank(gbk_path, ret.seqs = FALSE)
# 提取所有CDS的特征数据框
cds_df <- as.data.frame(cds(gb))
# 预处理目标字段
processed <- lapply(1:nrow(cds_df), function(i){
  # 提取各字段，缺失值替换为空字符串
  gene <- ifelse(is.null(cds_df$gene[i]) || is.na(cds_df$gene[i]), "", cds_df$gene[i])
  product <- ifelse(is.null(cds_df$product[i]) || is.na(cds_df$product[i]), "", cds_df$product[i])
  locus_tag <- ifelse(is.null(cds_df$locus_tag[i]) || is.na(cds_df$locus_tag[i]), "", cds_df$locus_tag[i])
  old_locus_tag <- ifelse(is.null(cds_df$old_locus_tag[i]) || is.na(cds_df$old_locus_tag[i]), "", cds_df$old_locus_tag[i])
  # 从inference字段提取RefSeq号
  refseq <- str_extract(cds_df$inference[i], "RefSeq:([A-Z0-9_.]+)") %>% str_remove("RefSeq:")
  refseq <- ifelse(is.na(refseq), "", refseq)
  protein_id <- ifelse(is.null(cds_df$protein_id[i]) || is.na(cds_df$protein_id[i]), "", cds_df$protein_id[i])
  # 提取坐标，把..替换为:
  coord <- str_extract(cds_df$location[i], "\\d+\\.\\.\\d+") %>% str_replace("\\.\\.", ":")
  coord <- ifelse(is.na(coord), "", coord)
  # 拼接翻译序列，去掉换行和空格
  translation <- str_remove_all(cds_df$translation[i], "\\s+")
  
  # 按要求拼接标题行
  header <- paste0(">", paste(gene, product, locus_tag, old_locus_tag, refseq, protein_id, coord, sep = "|"))
  # 返回标题和翻译序列
  return(c(header, translation))
})
# 展开为向量，准备写入文件
output <- unlist(processed)
# 写入输出文件，替换为你的输出路径
writeLines(output, "cds_output.fasta")

3. 多gbk文件批量处理 #

如果需要批量处理文件夹下所有gbk文件，使用以下逻辑：

# 替换为你的gbk文件所在文件夹路径
gbk_dir <- "your_gbk_dir/"
# 提取文件夹下所有.gbk后缀的文件
gbk_files <- list.files(gbk_dir, pattern = "\\.gbk$", full.names = TRUE)
# 循环处理所有文件
all_output <- lapply(gbk_files, function(f){
  gb <- readGenBank(f, ret.seqs = FALSE)
  cds_df <- as.data.frame(cds(gb))
  # 同上单文件的字段处理逻辑
  processed <- lapply(1:nrow(cds_df), function(i){
    gene <- ifelse(is.null(cds_df$gene[i]) || is.na(cds_df$gene[i]), "", cds_df$gene[i])
    product <- ifelse(is.null(cds_df$product[i]) || is.na(cds_df$product[i]), "", cds_df$product[i])
    locus_tag <- ifelse(is.null(cds_df$locus_tag[i]) || is.na(cds_df$locus_tag[i]), "", cds_df$locus_tag[i])
    old_locus_tag <- ifelse(is.null(cds_df$old_locus_tag[i]) || is.na(cds_df$old_locus_tag[i]), "", cds_df$old_locus_tag[i])
    refseq <- str_extract(cds_df$inference[i], "RefSeq:([A-Z0-9_.]+)") %>% str_remove("RefSeq:")
    refseq <- ifelse(is.na(refseq), "", refseq)
    protein_id <- ifelse(is.null(cds_df$protein_id[i]) || is.na(cds_df$protein_id[i]), "", cds_df$protein_id[i])
    coord <- str_extract(cds_df$location[i], "\\d+\\.\\.\\d+") %>% str_replace("\\.\\.", ":")
    coord <- ifelse(is.na(coord), "", coord)
    translation <- str_remove_all(cds_df$translation[i], "\\s+")
    header <- paste0(">", paste(gene, product, locus_tag, old_locus_tag, refseq, protein_id, coord, sep = "|"))
    return(c(header, translation))
  })
  return(unlist(processed))
})
# 合并所有结果并写入
writeLines(unlist(all_output), "all_cds_output.fasta")

注意事项 #

• 若部分gbk文件的字段命名存在差异，可先打印colnames(cds_df)查看所有可用字段，调整代码中的字段名即可
• 处理超大gbk文件时，可分批读取避免内存占用过高
• 缺省字段会自动填充为空字符串，不会中断处理流程

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