用 R 提取基因座的基因 ID

Question

我有几个 100 个具有开始和结束位置的基因组区域。我想从参考文件中提取每个区域的基因 ID（参考文件包含数千个基因及其开始和结束位置）。我提取了基因的区域信息，例如，输出文件告诉我基因是否存在于区域一或区域二等。然而，每个区域都包含许多基因，并且有多个基因位于多个区域。我想要一个输出文件，将所有基因的 ID 放在该区域旁边的单元格中（每个基因 ID 用逗号分隔）。如果一个基因存在于多个区域，它将与该区域的其他基因一起出现在与这些区域相对应的多个细胞中。这可以用 R 代码完成吗？请帮忙。

示例输入。

区域信息文件

Region  Chr Start   End
1   1A  1   12345
2   1A  23456   34567
3   2A  1234    23456
***
1830 7D 123     45678

基因信息文件

Gene    Chr Start   End
GeneID1 1A  831 1437
GeneID2 1A  1487    2436
GeneID3 1B  2665    5455
***
GeneID10101 7D  13456  56789

RequiredOutPutFile

Region  Chr Start   End   Gene
1       1A  1       12345 GeneID1, GeneID2, GeneID5, GeneID6 ***
***
1830    7D  123     45689 GeneID7, GeneID100, GeneID200 ***

Answer 1

您正在寻找的通常是通过 Bioconductor 上出色的 GenomicRanges 包来实现的。您可以从区域创建 GenomicRanges 对象，并使用 findOverlaps() 函数查找它们重叠的位置。要控制什么算作重叠，请参阅 findOverlaps(..., type="") 中的 "type" 参数。这是 R 中的一个例子。我相信有更优雅的方法可以做到这一点，但这就是我在午餐时想到的：

# install.packages('readr')
library(readr)
# install.packages('BiocManager')
# BiocManager::install('GenomicRanges')
library(GenomicRanges)

# Read data into R. Needs to have column names "chr", "start" and "end"
x <- read_csv('~/Downloads/regions.csv')
y <- read_csv('~/Downloads/genes.csv')

# Set up two GenomicRanges objects
gr_x <- makeGRangesFromDataFrame(x, keep.extra.columns = TRUE)
gr_y <- makeGRangesFromDataFrame(y, keep.extra.columns = TRUE)

# Overlap the regions containing genes with the trait data
ovl <- findOverlaps(gr_y, gr_x, type = "any", select = "all", ignore.strand = TRUE)
# Group hits by trait regions
hits_by_region <- split(from(ovl), to(ovl))

# Create a data.frame that matches a trait region index to a string of genes
hits_df <- do.call(rbind, lapply(seq_along(hits_by_region), function(f) {
  idx <- as.integer(names(hits_by_region[f]))
  genes <- gr_y$Gene[hits_by_region[[f]]]
  data.frame(Index = idx, Genes = paste(genes, collapse=','), stringsAsFactors = FALSE)
}))

# Add the string of genes as metadata to the GRanges object
gr_x$Genes <- ""
gr_x$Genes[hits_df$Index] <- hits_df$Genes

# Convert back to data.frame (if needed)
genes_by_regio_df <- as.data.frame(gr_x)

用 R 提取基因座的基因 ID

Extracting gene IDs for loci with R

annotations

r

示例输入。

区域信息文件

基因信息文件

RequiredOutPutFile