Извлечение идентификаторов генов для локусов с помощью R

Question

У меня есть несколько 100 геномных областей с начальными конечными положениями. Я хочу извлечь идентификаторы генов в каждом регионе из эталонного файла (эталонный файл содержит тысячи генов с их начальной конечной позицией). Я извлек информацию о регионе генов, например, выходной файл сообщает мне, присутствует ли ген в области один или два, и так далее. Однако каждый регион содержит много генов, и есть несколько генов, которые лежат в нескольких регионах. Мне нужен выходной файл, который бы помещал идентификаторы всех генов в ячейку рядом с этой областью (каждый идентификатор гена отделялся запятой). Если ген присутствует в нескольких областях, он появится в нескольких клетках, соответствующих этим областям, наряду с другими генами этой области. Можно ли это сделать с помощью кода R? Пожалуйста, помогите.

Пример ввода.

RegionInfoFile

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

GeneInfoFile

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

То, что вы ищете, обычно достигается с помощью превосходного пакета GenomicRanges на Bioconductor. Вы можете создавать GenomicRanges объекты из регионов и находить, где они перекрываются, используя функцию findOverlaps(). Чтобы контролировать то, что считается перекрытием, смотрите аргумент «type» в findOverlaps(..., 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)