Hi everyone,

I'm learning the course ph525.5x, which is a fantastic course to learn Bioconductor I think. When I repeated the code "res <- precede(erbs, ghs)" as at 5:10 of this video, I could not get the same result with the instructor. Then I looked though precede() and found that it does not perform exactly as what is described in its help page. There is a sentence in "Details" of ?`precede,GenomicRanges,GenomicRanges-method`:

"x on * strand can match to ranges on any of +, - or * strands. In the case of a tie the first range by order is chosen."

First look at the source code.

> getMethod(precede, signature(x="GenomicRanges", subject="GenomicRanges"))
Method Definition:

function (x, subject = x, ...)
{
    .local <- function (x, subject, select = c("arbitrary", "all"),
        ignore.strand = FALSE)
    {
        select <- match.arg(select)
        .GenomicRanges_findPrecedeFollow(x, subject, select,
            ignore.strand, "precede")
    }
    .local(x, subject, ...)
}
<environment: namespace:GenomicRanges>

Signatures:
        x               subject       
target  "GenomicRanges" "GenomicRanges"
defined "GenomicRanges" "GenomicRanges"

> GenomicRanges:::.GenomicRanges_findPrecedeFollow
function (query, subject, select, ignore.strand, where = c("precede",
    "follow"))
{
    if (!length(query) || !length(subject))
        return(Hits(nLnode = length(query), nRnode = length(subject),
            sort.by.query = TRUE))
    leftOf <- "precede" == match.arg(where)
    if (ignore.strand)
        strand(query) <- strand(subject) <- "+"
    if (leftOf) {
        plusfun <- function(xstart, xend, ystart, yend, sentinel) .GenomicRanges_findNearest0(xend,
            ystart, sentinel, leftOf)
        minusfun <- function(xstart, xend, ystart, yend, sentinel) .GenomicRanges_findNearest0(xstart,
            yend, sentinel, !leftOf)
    }
    else {
        plusfun <- function(xstart, xend, ystart, yend, sentinel) .GenomicRanges_findNearest0(xstart,
            yend, sentinel, leftOf)
        minusfun <- function(xstart, xend, ystart, yend, sentinel) .GenomicRanges_findNearest0(xend,
            ystart, sentinel, !leftOf)
    }
    endq <- start(query) + width(query)
    ends <- start(subject) + width(subject)
    maxend <- max(max(endq), max(ends)) + 1
    lvls <- union(seqlevels(query), seqlevels(subject))
    offset <- setNames((seq_along(lvls) - 1) * maxend, lvls)
    stopifnot(typeof(offset) == "double")
    sentinel <- c(0, seq_along(lvls) * maxend)
    queryOff <- unname(offset[as.character(seqnames(query))])
    queryStart <- start(query) + queryOff
    queryEnd <- end(query) + queryOff
    qid <- seq_along(query)
    subjectOff <- unname(offset[as.character(seqnames(subject))])
    subjectStart <- start(subject) + subjectOff
    subjectEnd <- end(subject) + subjectOff
    spid <- which(strand(subject) != "-")
    smid <- which(strand(subject) != "+")
    idx <- which(strand(query) == "+")
    phit <- plusfun(queryStart[idx], queryEnd[idx], subjectStart[spid],
        subjectEnd[spid], sentinel)
    phit <- .Hits(qid[idx][queryHits(phit)], spid[subjectHits(phit)])
    idx <- which(strand(query) == "-")
    mhit <- minusfun(queryStart[idx], queryEnd[idx], subjectStart[smid],
        subjectEnd[smid], sentinel)
    mhit <- .Hits(qid[idx][queryHits(mhit)], smid[subjectHits(mhit)])
    idx <- which(strand(query) == "*")
    bhit <- local({
        qid <- qid[idx]
        phit <- plusfun(queryStart[idx], queryEnd[idx], subjectStart,
            subjectEnd, sentinel)
        .Hits(qid[queryHits(phit)], subjectHits(phit), length(query),
            length(subject))
    })
    qryHits <- c(queryHits(phit), queryHits(mhit), queryHits(bhit))
    subjHits <- c(subjectHits(phit), subjectHits(mhit), subjectHits(bhit))
    if ("arbitrary" == select) {
        hits <- integer()
        hits[length(query)] <- NA_integer_
        idx <- !duplicated(qryHits)
        hits[qryHits[idx]] <- subjHits[idx]
    }
    else {
        hits <- .Hits(qryHits, subjHits, length(query), length(subject))
    }
    hits
}
<environment: namespace:GenomicRanges>

In fact, that sentence is right but misleading. For the ranges of x on "*" strand, what precede() does is identical to precede(., ignore.strand = TRUE), say, just to treat all ranges x on * strand, and all ranges of subject as they are all on the "+" strand. However, what users desire (also what the instructor desire) is to treat the x as they were on "+" and do precede() with ranges of subject on "+" and "*", then treat x as they were on "-" and do precede() with ranges of subject on "-" and "*", then return the range for each x that is nearest to x among which x precede on either strand. We want it that way because it is very useful to find the nearest gene that each TFBS peak precede.

I write some code to realize that demand.

precedeGenes <- function(x, genes) {
  xPlus <- GRanges(x, strand = "+")
  xMinus <- GRanges(x, strand = "-")
  idxGenesPlus <- precede(xPlus, genes, ignore.strand = FALSE)
  idxGenesMinus <- precede(xMinus, genes, ignore.strand = FALSE)
  disPlus <- rep(Inf, length(x))
  disMinus <- rep(Inf, length(x))
  disPlus[!is.na(idxGenesPlus)] <- distance(xPlus[!is.na(idxGenesPlus)], genes[na.omit(idxGenesPlus)])
  disMinus[!is.na(idxGenesMinus)] <- distance(xMinus[!is.na(idxGenesMinus)], genes[na.omit(idxGenesMinus)])
  dis <- pmin(disPlus, disMinus)
  idxGenes <- ifelse(disPlus > disMinus, idxGenesMinus, idxGenesPlus)
  return(idxGenes)
}

precedeGenes(erbs, ghs)

 [1] 22817 16173 21772  7870 20199 22131 13257 20469 21943  8564
[11] 19061 15836 21282  8959   278 14620  1657 22774 21746 10777
[21] 14568 20910  6626 12830  9838 19991 21643 15740 12911 13326
[31]  8303 20834   934 10485 23027  4443  6613 21335 21125 14508
[41] 15979  2693  4122  9485 17665 18376 17932  2696 16990 16952
[51] 15735  9386  1642  3634  6037 20891 15062  1779  3316 11746
[61] 12168 21689 22406 19979 14099 13558  3653  6077 21280  9705
[71] 13095  3322  6343  1553  9552

It is right but very fragile in the wild. So I wish the precede() could be modified to do the right thing. Now I have some idea but not complete. In the bold Italic,

phit <- plusfun(queryStart[idx], queryEnd[idx], subjectStart, 
            subjectEnd, sentinel)

can be replaced by

phit <- plusfun(queryStart[idx], queryEnd[idx], subjectStart[spid], 
        subjectEnd[spid], sentinel)

mhit <- minusfun(queryStart[idx], queryEnd[idx], subjectStart[smid], 
        subjectEnd[smid], sentinel)

Then use some code to combine them appropriately.

You can reproduce results in the video as below.

# library(devtools)
# install_github("genomicsclass/ERBS")
library(ERBS)
library(Homo.sapiens)
data(HepG2)
data(GM12878)
res <- findOverlaps(HepG2, GM12878)
index <- queryHits(res)
erbs <- HepG2[index, ]
erbs <- granges(erbs)
ghs <- genes(Homo.sapiens)
res <- precede(erbs, ghs)
res
 [1] 16845 16173 21772 13617 18834 15172  4543 22522 21943 22493
[11]   846 15836  3903 21053 18112 17682 19650 22112 21746 22839
[21] 14568  9043 20477 12830  9841  5379 21643 14441 12911   553
[31]  7980  6028 22699  4480 11293 20120 11859 17016 14199  6719
[41] 22946  2693  6923 20570  6680 15953 17932 20311  8049 11969
[51]  3738  2914  1642  4173  6443  4446 15062 11268 16695  3496
[61]  8458   568 22406  3201 14099 13558  8545  7569 13718 16702
[71] 13095 17686 19348 12849  9552

precede(erbs, ghs, ignore.strand = TRUE)
 [1] 16845 16173 21772 13617 18834 15172  4543 22522 21943 22493
[11]   846 15836  3903 21053 18112 17682 19650 22112 21746 22839
[21] 14568  9043 20477 12830  9841  5379 21643 14441 12911   553
[31]  7980  6028 22699  4480 11293 20120 11859 17016 14199  6719
[41] 22946  2693  6923 20570  6680 15953 17932 20311  8049 11969
[51]  3738  2914  1642  4173  6443  4446 15062 11268 16695  3496
[61]  8458   568 22406  3201 14099 13558  8545  7569 13718 16702
[71] 13095 17686 19348 12849  9552

identical(precede(erbs, ghs, ignore.strand = TRUE), precede(erbs, ghs, ignore.strand = FALSE)) # should not be TRUE I think
[1] TRUE

precedeGenes <- function(x, genes) {
  xPlus <- GRanges(x, strand = "+")
  xMinus <- GRanges(x, strand = "-")
  idxGenesPlus <- precede(xPlus, genes, ignore.strand = FALSE)
  idxGenesMinus <- precede(xMinus, genes, ignore.strand = FALSE)
  disPlus <- rep(Inf, length(x))
  disMinus <- rep(Inf, length(x))
  disPlus[!is.na(idxGenesPlus)] <- distance(xPlus[!is.na(idxGenesPlus)], genes[na.omit(idxGenesPlus)])
  disMinus[!is.na(idxGenesMinus)] <- distance(xMinus[!is.na(idxGenesMinus)], genes[na.omit(idxGenesMinus)])
  dis <- pmin(disPlus, disMinus)
  idxGenes <- ifelse(disPlus > disMinus, idxGenesMinus, idxGenesPlus)
  return(idxGenes)
}
precedeGenes(erbs, ghs) # the result is identical to that in the video

 [1] 22817 16173 21772  7870 20199 22131 13257 20469 21943  8564
[11] 19061 15836 21282  8959   278 14620  1657 22774 21746 10777
[21] 14568 20910  6626 12830  9838 19991 21643 15740 12911 13326
[31]  8303 20834   934 10485 23027  4443  6613 21335 21125 14508
[41] 15979  2693  4122  9485 17665 18376 17932  2696 16990 16952
[51] 15735  9386  1642  3634  6037 20891 15062  1779  3316 11746
[61] 12168 21689 22406 19979 14099 13558  3653  6077 21280  9705
[71] 13095  3322  6343  1553  9552

sessionInfo()
R version 3.3.2 (2016-10-31)
Platform: x86_64-w64-mingw32/x64 (64-bit)
Running under: Windows 7 x64 (build 7601) Service Pack 1

locale:
[1] LC_COLLATE=Chinese (Simplified)_People's Republic of China.936 
[2] LC_CTYPE=Chinese (Simplified)_People's Republic of China.936   
[3] LC_MONETARY=Chinese (Simplified)_People's Republic of China.936
[4] LC_NUMERIC=C                                                   
[5] LC_TIME=Chinese (Simplified)_People's Republic of China.936    

attached base packages:
[1] parallel  stats4    stats     graphics  grDevices utils    
[7] datasets  methods   base     

other attached packages:
 [1] Homo.sapiens_1.3.1                     
 [2] TxDb.Hsapiens.UCSC.hg19.knownGene_3.2.2
 [3] org.Hs.eg.db_3.4.0                     
 [4] GO.db_3.4.0                            
 [5] OrganismDbi_1.16.0                     
 [6] GenomicFeatures_1.26.2                 
 [7] GenomicRanges_1.26.1                   
 [8] GenomeInfoDb_1.10.1                    
 [9] AnnotationDbi_1.36.0                   
[10] IRanges_2.8.1                          
[11] S4Vectors_0.12.1                       
[12] Biobase_2.34.0                         
[13] BiocGenerics_0.20.0                    
[14] ERBS_1.0                               

loaded via a namespace (and not attached):
 [1] graph_1.52.0               Rcpp_0.12.8               
 [3] XVector_0.14.0             zlibbioc_1.20.0           
 [5] GenomicAlignments_1.10.0   BiocParallel_1.8.1        
 [7] lattice_0.20-34            tools_3.3.2               
 [9] grid_3.3.2                 SummarizedExperiment_1.4.0
[11] DBI_0.5-1                  RBGL_1.50.0               
[13] digest_0.6.10              Matrix_1.2-7.1            
[15] rtracklayer_1.34.1         bitops_1.0-6              
[17] RCurl_1.95-4.8             biomaRt_2.30.0            
[19] memoise_1.0.0              RSQLite_1.1-1             
[21] BiocInstaller_1.24.0       Biostrings_2.42.1         
[23] Rsamtools_1.26.1           XML_3.98-1.5 

Thanks in advance.

Can

Hi Can,

As you can imagine, people have different opinions on 'right behavior' of a function. precede() and follow() handle strand the way they do because of user input at the time they were written.

It sounds like you want to make specific strand comparisons in a certain order. Instead of modifying the function to meet a specific use case it's better to make several calls to the function while imposing your strand restrictions as necessary.

This statement you made above is not true (or I've misinterpreted it):
"... For the ranges of x on "*" strand, what precede() does is identical to precede(., ignore.strand = TRUE), say, just to treat all ranges x on * strand, and all ranges of subject as they are all on the "+" strand. ..."

As stated in the man page, a '*' strand range can match to '+', '-' or '*'. In the case of a tie, the first range by order is chosen.

The '*' query matches to the "+' range subject because it's first by order.

query <- GRanges("A", IRanges(start=5, width=1), strand="*")
subject <- GRanges("A", IRanges(start=c(10, 10), width=1), strand=c("+", "-"))
> precede(query, subject)
[1] 1

Reversing the order of the subject now the '-' range is chosen because it's first by order.

> precede(query, rev(subject))
[1] 1

If you feel the documentation does not accurately reflect the behavior please cite the section and I'll try to make it more clear.

Thanks.
Valerie