Hi, 

I am using nearest to annotate my genomic intervals (breakpoints from whole genome sequencing). Those are not ChIP-seq peaks and the strandness for each interval is critical to me.
6 column bed file format:

chr1  10000 100200  breakpoint_1  .   +
chr2  20000 200100 breakpoint_2  .    -

I want to find the closest gene that is upstream of the breakpoint.

library(GenomicFeatures)
hg19.refseq.db <- makeTxDbFromUCSC(genome="hg19", table="refGene")
refseq.genes<- genes(hg19.refseq.db)

for breakpoint_1, I only want to search genes with start smaller than 10000.

for breakpoint_2, I only want to search genes with end greater than 200100.

Note that, I am not saying the strandness of the genes (as mentioned here http://ygc.name/2014/01/14/bug-of-r-package-chippeakanno/), but rather the strandness of the intervals.

It looks like `distanceToNearest()` always returns [an obsolute value](distance to TSS upstream/downstream) for distance no matter it is upstream or downstream of the subject hits. We can [get around with it](https://github.com/genomicsclass/labs/blob/master/bioc/grangesERBSExample.Rmd#finding-nearest-gene-for-each-binding-event). Instead, one can use follow and precede.

The bumphunter bioconductor package by Rafael A. Irizarry et.al has a function called [`annoateNearest`](http://finzi.psych.upenn.edu/library/bumphunter/html/annotateNearest.html) can better annotate the nearest distances considering the strandness, but it DID NOT restrict nearest to find only upstream or downstream relative to the query.

I was trying to avoid using bedtools and put everything in R for the sake of reproducibility, although [bedtools closest](http://bedtools.readthedocs.org/en/latest/content/tools/closest.html) documentation has a much better explannation of all the cases I want to use. 

In a word, for GRanges, is it possible to find the upstream nearest or downstream nearest as bedtools closest -D -id

and bedtools closest -D -iu

Thanks very much!

Ming

I think you can just do something like:

upstreamGenes <- follow(breakpoints, unstrand(refseq.genes))

If you also want distance, then do:

distance(breakpoints, refseq.genes[upstreamGenes])

But beware that there may be a breakpoint that has no upstream gene.

Hi,

You could also have a try with ChIPpeakAnno like this,

library(ChIPpeakAnno)
gr <- GRanges(c("chr1", "chr2"), IRanges(c(10000, 20000), c(100200, 200100), names=c("breakpoint_1", "breakpoint_2")), strand=c("+", "-")) ## toGRanges("bed.filename", format="BED")
library(GenomicFeatures)
hg19.refseq.db <- makeTxDbFromUCSC(genome="hg19", table="refGene")
genes <- annoGR(hg19.refseq.db, feature="gene")
 gr.anno <- annotatePeakInBatch(gr, AnnotationData=genes, output="nearestLocation", ignore.strand=FALSE)

------------------------------------------------------------------------------------------------------------------------------------------

About the post http://ygc.name/2014/01/14/bug-of-r-package-chippeakanno/,  I want to explain a little more.

In TSS.human.GRCh37, the annotations for NIPAL3(ENSG00000001461) and STPG1(ENSG00000001460) are,

GRanges object with 2 ranges and 1 metadata column:
                  seqnames               ranges strand |                                                        description
                     <Rle>            <IRanges>  <Rle> |                                                        <character>
  ENSG00000001461        1 [24742285, 24799466]      + |       NIPA-like domain containing 3 [Source:HGNC Symbol;Acc:25233]
  ENSG00000001460        1 [24683490, 24743424]      - | UPF0490 protein C1orf201  [Source:UniProtKB/Swiss-Prot;Acc:Q5TH74]
  -------
  seqinfo: 72 sequences from an unspecified genome; no seqlengths

At the time we generate the annotation file, STPG1 is still not named yet. And it's range is from 24683490 to 24743424. If you calculate the distance of TSS to IRanges(24736757, 24737528), you will get:

24742285-24737528= 4757 for NIPAL3
24743424-24737528= 5896 for C1orf201 (named as STPG1 later).

And Guangchuan ignored this fact and changed the annotation file to TxDb.Hsapiens.UCSC.hg19.knownGene and want to say our result is not true.
However, if you try to change the annotations, you will get,
-----------------------
> annoData <- annoGR(TxDb.Hsapiens.UCSC.hg19.knownGene)
> annotatePeakInBatch(peak, AnnotationData=annoData)
GRanges object with 1 range and 9 metadata columns:
           seqnames               ranges strand |        peak     feature start_position end_position feature_strand insideFeature distancetoFeature
              <Rle>            <IRanges>  <Rle> | <character> <character>      <integer>    <integer>    <character>      <factor>         <numeric>
  X1.90529     chr1 [24736757, 24737528]      * |           1       90529       24683489     24741587              -        inside              4830
           shortestDistance fromOverlappingOrNearest
                  <integer>              <character>
  X1.90529             4059          NearestLocation
  -------
  seqinfo: 1 sequence from an unspecified genome; no seqlengths
> get("90529", org.Hs.egSYMBOL)
[1] "STPG1"
---------------------------
You will see how Guangchuan play with the data. The range of STPG1 is very different from the data we saved at the release time (we do not change this just for repeatability of our package).

Please understand that his post spend a lot of our time to explain why this happens.