I have the following problem: I want to find the next 5 gene in (chr1_gene) on both side of a interval give from the GRanges (gr3) and generate a dataframe containing chr start end (from gr3) of the interval and i column for each ensembl_transcript_id with eventually NA.

Thanks a lot for suggestions!

df <- data.frame(chrom=c("chr1","chr1"), start=c(5087459, 9995206 ),   end=c(5097899, 10015020 ))
gr3 <- as(df, "GRanges")


library(biomaRt)
ensembl = useEnsembl(biomart="ensembl", dataset="hsapiens_gene_ensembl")
chr1_genes <- getBM(attributes=c('ensembl_gene_id',
'ensembl_transcript_id','hgnc_symbol','chromosome_name','start_position','end_position'), filters =
'chromosome_name', values ="1", mart = ensembl)

head(chr1_gene)

ensembl_gene_id ensembl_transcript_id hgnc_symbol chromosome_name start_position end_position
1 ENSG00000231510       ENST00000443270                           1        5086459      5090899
2 ENSG00000162444       ENST00000315901        RBP7               1        9997206     10016020
3 ENSG00000162444       ENST00000294435        RBP7               1        9997206     10016020
4 ENSG00000270171       ENST00000602640                           1        7693124      7694844
5 ENSG00000225643       ENST00000412797                           1       25581478     25590356
6 ENSG00000116497       ENST00000530710     S100PBP               1       32816767     32858879

I have fond this  but is quite slow:

enh25gens <- function(gr, GR){

idxf.1 <- follow(gr, subject = GR)
idxf.2 <- follow(GR[idxf.1],subject = GR)
idxf.3 <- follow(GR[idxf.2],subject = GR)
idxf.4 <- follow(GR[idxf.3],subject = GR)
idxf.5 <- follow(GR[idxf.4],subject = GR)

fol.gene.list <- list()
fol.gene.list <- append( fol.gene.list, list(mcols(GR[idxf.1])$ensembl_gene_id))
fol.gene.list <- append( fol.gene.list, list(mcols(GR[idxf.2])$ensembl_gene_id))
fol.gene.list <- append( fol.gene.list, list(mcols(GR[idxf.3])$ensembl_gene_id))
fol.gene.list <- append( fol.gene.list, list(mcols(GR[idxf.4])$ensembl_gene_id))
fol.gene.list <- append( fol.gene.list, list(mcols(GR[idxf.5])$ensembl_gene_id))

fol.gene <- unlist(fol.gene.list)

idxf.1 <- precede(gr, subject = GR)
idxf.2 <- precede(GR[idxf.1],subject = GR)
idxf.3 <- precede(GR[idxf.2],subject = GR)
idxf.4 <- precede(GR[idxf.3],subject = GR)
idxf.5 <- precede(GR[idxf.4],subject = GR)

pre.gene.list <- list()
pre.gene.list <- append( pre.gene.list, list(mcols(GR[idxf.1])$ensembl_gene_id))
pre.gene.list <- append( pre.gene.list, list(mcols(GR[idxf.2])$ensembl_gene_id))
pre.gene.list <- append( pre.gene.list, list(mcols(GR[idxf.3])$ensembl_gene_id))
pre.gene.list <- append( pre.gene.list, list(mcols(GR[idxf.4])$ensembl_gene_id))
pre.gene.list <- append( pre.gene.list, list(mcols(GR[idxf.5])$ensembl_gene_id))

pre.gene <- unlist(pre.gene.list)

list.enh2gene <- unlist(c(as.data.frame(gr[1]), pre.gene,fol.gene ))
return(list.enh2gene)
}

df <-  do.call(rbind.data.frame, lapply(gr3, function(x) enh25gens((x),chr1_genes.GR))  )
colnames(df) <- c("chr", "start", "end", "width", "strand", 
                  "pre1", "pre2", "pre3", "pre4", "pre5", "fol1", "fol2", "fol3", "fol4", "fol5")

Here is the function:

findKNN <- function(query, subject, k=5L) {
    seqlevels(subject) <- seqlevels(query)
    
    starts <- with(subject, GRanges(seqnames, IRanges(start, width=1L)))
    ends <- with(subject, GRanges(seqnames, IRanges(end, width=1L)))

    start_ord <- order(starts)
    end_ord <- order(ends)

    starts <- starts[start_ord]
    ends <- ends[end_ord]

    phits <- precede(query, starts)
    fhits <- follow(query, ends)
    
    seqends <- end(seqnames(starts))[findRun(phits, seqnames(starts))]
    phits[is.na(phits)] <- 1L
    seqends[is.na(seqends)] <- 0L
    pwindows <- restrict(IRanges(phits, width = k), end=seqends)

    seqstarts <- start(seqnames(starts))[findRun(fhits, seqnames(starts))]
    seqstarts[is.na(seqstarts)] <- 1L
    fhits[is.na(fhits)] <- 0L
    fwindows <- restrict(IRanges(end=fhits, width = k), seqstarts)
    
    w <- pc(extractList(start(starts), pwindows) - end(query),
            end(query) - extractList(end(ends), fwindows))

    ans <- pc(extractList(start_ord, pwindows), extractList(end_ord, fwindows))
    ans[phead(order(w), k)]
}

The function can be slow in the release branch, because of the suboptimal order() method for the list at the end. There is an optimized order() method in devel.

First, you should get your genes into a GRanges via a TxDb from makeTxDbFromBiomart(). I just used UCSC and generated 10,000 query ranges across the genome.

library(BSgenome.Hsapiens.UCSC.hg19)
tiles <- tileGenome(seqinfo(BSgenome.Hsapiens.UCSC.hg19), ntile=10000)
query <- resize(unlist(tiles), 10)

library(TxDb.Hsapiens.UCSC.hg19.knownGene)
subject <- genes(TxDb.Hsapiens.UCSC.hg19.knownGene)

system.time(findKNN(query, subject))

   user  system elapsed
  0.153   0.010   0.163

So at least in devel it takes 0.15 seconds for that query. Not too bad, but that's only 10k queries. Caveat: strand is ignored and genes overlapping the query will not be detected (just as with precede and follow). The return value is a List of subscripts into subject. We can for example extract the gene IDs:

ids <- extractList(subject$gene_id, knn)

You wanted a matrix, and we can make one very quickly thanks to partitionings:

as.matrix.AtomicList <- function(x, ...) {
    p <- PartitioningByEnd(x)
    m <- matrix(nrow=length(x), ncol=max(width(p)))
    ind <- cbind(togroup(p), as.integer(IRanges(1, width(p))))
    m[ind] <- unlist(x, use.names=FALSE)
    m
}
setMethod("as.matrix", "AtomicList", function(x) as.matrix.AtomicList(x))

knn_matrix <- as.matrix(knn)

I will put that as.matrix method in IRanges. The findKNN one might need some more love.