Hi Duncan,

Here are a couple of things that could explain the differences you got between the method used in the generegulation workflow and your own code:

1) By default promoters() will extract the regions located 2000 bp upstream and 200 bp downstream of the TSS (quickly check this with args(promoters)). Since you said you wanted 1000 bp upstream and 500 bp downstream of the TSS, make sure you specify this in your call to promoters():

human_promoters <- promoters(TxDb.Hsapiens.UCSC.hg19.knownGene,
                             upstream=1000, downstream=500)

2) I guess the purpose of your last step is to keep only the hits returned by matchPWM() that are in the promoter regions. Instead of using intersect() for this (whose exact semantic can be a little bit tricky to grasp), I would suggest you do:

promoter_hit <- genome_hit[genome_hit %within% human_promoters]

Erratum: need to use %within% or %over% instead of %in% in the above code. Edited it to use %within%.

Unlike interesect(), this won't truncate TF binding sites that are not fully contained inside a promoter region, or merge TF binding sites that overlap. These are rare events but they can happen. In other words, by subsetting genome_hit, the ranges in promoter_hit are guaranteed to be a subset of the ranges in genome_hit, and with no alteration of the ranges.

Anyway, since your ultimate goal is to get the list of transcripts/genes for your TF of interest, some of these steps can be avoided. Also, since the string matching step is the expensive one, it's worth trying to minimize that cost by calling matchPWM()/countPWM() on the promoter regions only (instead of the entire genome). The getTFtranscripts() function below does that:

library(BSgenome)

## Returns a data.frame with 1 row per transcript and 2 columns (tx_name
## and gene_id).
getTFtranscripts <- function(pwm, txdb, genome,
                             upstream=2000, downstream=200, min.score="80%")
{
    prom <- suppressWarnings(
        promoters(txdb, upstream=upstream, downstream=downstream,
                        columns=c("tx_name", "gene_id")))
    prom <- trim(prom)
    genome <- getBSgenome(genome)
    prom_seqs <- getSeq(genome, prom)
    hits_per_prom <- suppressWarnings(sapply(
        seq_along(prom_seqs),
        function(i) countPWM(pwm, prom_seqs[[i]], min.score=min.score)
    ))
    ans <- mcols(prom)[hits_per_prom != 0L, ]
    ans[ , "gene_id"] <- as.character(ans[ , "gene_id"])
    as.data.frame(ans)
}

## Then:
library(TxDb.Hsapiens.UCSC.hg19.knownGene)
txdb <- TxDb.Hsapiens.UCSC.hg19.knownGene
pcm.rest.transcripts <- getTFtranscripts(pcm.rest, txdb, "hg19",
                            upstream=1000, downstream=500, min.score="85%")

Finally, also please note that what you pass to getTFtranscripts() (or to matchPWM(), or to countPWM()) should be a Position Weight Matrix (PWM) and not a Position Frequency Matrix (PFM). matchPWM() and countPWM() are really supposed to be used with a PWM. Note that a PFM can easily be turned into a PWM by just calling the PWM() function on it. See ?PWM for more information about this. The PFM -> PWM transformation performed by PWM() is based on the method described in the Wasserman & Sandelin paper (see ?PWM for the reference to the paper). It's probably not a big deal since the results you will get when using a PWM instead of a PFM will generally be very close or even the same. But sometimes they won't.

Cheers,

H.