求助：如何按Bed文件中单个区域逐一获取对应范围内的基因

我懂你遇到的问题——你想要每个Bed区域和它对应的基因一一绑定，而不是得到一堆杂乱无章的基因列表。下面分别针对bedtools和biomaMart给出具体的解决办法：

bedtools默认输出可能会让匹配的基因和区域信息脱节，但只要加上两个关键参数，就能轻松把原始区域和匹配基因的信息绑定在一起：

假设你的查询区域文件是query_regions.bed，基因注释Bed文件（比如从Ensembl/NCBI下载的标准基因坐标文件）是genes.bed，执行以下命令：

bedtools intersect -wa -wb -a query_regions.bed -b genes.bed > region_gene_mappings.bed

参数解释：

• -wa：完整保留你的查询区域（文件a）的每一行内容，包括你Bed里的MIRxxx标识
• -wb：完整保留匹配到的基因（文件b）的每一行内容

输出的每一行都会同时包含你的原始区域信息和对应的基因信息，一眼就能看出哪个区域对应哪些基因。如果一个区域匹配到多个基因，会生成多行，每行对应一个基因-区域对；如果某个区域没有匹配到基因，就不会出现在结果里。

要是想保留所有查询区域（哪怕没有匹配到基因的），可以加上-loj（左外连接）参数：

bedtools intersect -wa -wb -loj -a query_regions.bed -b genes.bed > region_gene_mappings_with_null.bed

这样无匹配的区域也会出现在结果中，基因部分会用.填充。

biomaMart批量查询默认会返回所有匹配基因，但我们可以通过两种方式把基因和对应的区域关联起来：

方式1：批量查询后绑定区域标识 #

先给你的每个Bed区域加上唯一标识（就是你Bed里的MIRxxx名称），然后在查询时把这个标识和区域坐标一起传入，最后在结果中保留关联：

library(biomaRt)

# 读取你的Bed文件，给列命名方便处理
bed_df <- read.table("query_regions.bed", sep = " ", col.names = c("chrom", "start", "end", "region_id"))

# 连接到对应物种的Ensembl数据库（这里以人类为例，按需替换）
ensembl <- useMart("ensembl", dataset = "hsapiens_gene_ensembl")

# 批量查询，同时保留区域ID的关联
results <- getBM(
  attributes = c("external_gene_name", "chromosome_name", "start_position", "end_position"),
  filters = c("chromosome_name", "start", "end"),
  values = list(
    chromosome_name = bed_df$chrom,
    start = bed_df$start,
    end = bed_df$end
  ),
  mart = ensembl
)

# 将区域ID与查询结果关联（确保顺序匹配）
results$region_id <- rep(bed_df$region_id, sapply(1:nrow(bed_df), function(x) {
  sum(results$chromosome_name == bed_df$chrom[x] & 
      results$start_position >= bed_df$start[x] & 
      results$end_position <= bed_df$end[x])
}))

方式2：循环单个区域查询（更直观稳妥） #

如果担心批量查询的关联出错，可以循环遍历每个区域单独查询，直接把区域ID和对应基因绑定成列表：

library(biomaRt)

bed_df <- read.table("query_regions.bed", sep = " ", col.names = c("chrom", "start", "end", "region_id"))
ensembl <- useMart("ensembl", dataset = "hsapiens_gene_ensembl")

# 创建空列表存储结果
region_gene_map <- list()

for (i in 1:nrow(bed_df)) {
  current_region <- bed_df[i, ]
  # 查询当前区域内的基因
  genes <- getBM(
    attributes = "external_gene_name",
    filters = c("chromosome_name", "start", "end"),
    values = list(
      chromosome_name = current_region$chrom,
      start = current_region$start,
      end = current_region$end
    ),
    mart = ensembl
  )
  # 把区域ID和基因列表关联，无匹配则标记为NA
  region_gene_map[[current_region$region_id]] <- ifelse(nrow(genes) == 0, NA, genes$external_gene_name)
}

# 转换成数据框方便查看和后续处理
region_gene_df <- stack(region_gene_map)
colnames(region_gene_df) <- c("gene_name", "region_id")

这个方法的好处是绝对不会出现区域和基因的关联错误，每个区域的结果都独立处理，一目了然。

内容的提问来源于stack exchange，提问作者Emilio Mármol Sánchez