You asked
Does anyone know any tools that can intersect two bigbed files
...but it appears you really wanted to subset a single very large bigbed file based on the contents of a much smaller bed file.
I disagree with your summation "there is not a very efficient way to work directly with the bigbed file". I see two options.
The first is to repeatedly call bigBedTobed once for each row in your .bed file. If you have multiple cores, this can be parallelized using GNU parallel. The following does that, emitting each overlapping result with additional column being the index into your .bed file from which the result stems.
cat your.bed | parallel --lb -j --colsep '\t' bigBedToBed -chrom={1} -start={2} -end={3} your.bb /dev/stdout '|' sed -e 's/$/\\t${PARALLEL_SEQ}/'
Depending upon your application, there is a second approach which uses R/BioConductor's rtracklayer::import.bb which can be parallelized across (chunks of) your .bed file.
I considered these approaches in what I expect is an almost identical situation:
My bigbed holds genome-wide transcription factor binding site predictions from JASPAR Genome Browser tracks as produced by wassermanlab / JASPAR-UCSC-tracks.
My bed file holds the loci of ~23K putative enhancers as determined by an analysis of regions of accessible chromatin from ATAC-Seq and various ChIP marks.
In my case I wished to produce a sparse Matrix with one row for each putative enhancer, and one column for each JASPAR motif (which will form the basis of downstream enrichment analyses).
bbf<-BigBedFile('path/to/your.bb')
gr<-import.bed('path/to/your.bed')
## I found I had to harmonize seqlevels (e.g. chromosome identifiers) to make this work. YMMV
gr<-renameSeqlevels(gr,paste0('chr',gsub('^(K.*)\\.(\\d+)$','Un_\\1v\\2',seqlevels(gr))))
mm<-motifMatrix(gr,bbf)
which results in a sparse Matrix, viz:
> rownames(mm)<-as.character(gr) # optional.
> mm[1:5,1:5]
5 x 5 sparse Matrix of class "dgCMatrix"
MA0002.2 MA0003.4 MA0004.1 MA0006.1 MA0007.3
chr1:11156-12180 14 . 2 28 2
chr1:17836-20023 25 . 27 63 2
chr1:22495-24367 22 1 22 137 2
chr1:27533-28038 2 2 4 7 .
chr1:36549-37079 6 . 2 41 .
This depends upon the following definition:
motifMatrix<-function(gr,bbf,motifID=c('name','TFName')[[1]],mc.jobs=100) {
## produce a sparse Matrix with one row for each region in
## GenomicRanges gr, and one column for each JASPAR motif indicated
## as overlapping <gr> accorinding to BigBedFile <bbf>. Default to
## using the 'name' column as the motif identifier, but allow for
## possibly using another, such as 'TFName' (as used by 2022 version
## of
## https://jaspar.genereg.net/genome-tracks/#ucsc_tracks). Parallelize
## into mc.jobs queries, defaulting to 100, which was found suitable
## for ~25000 loci in <gr> against JASPAR matches to zebrafish
## danRer11 genome running on a few dozen cores. In general,
## knowing how "big" each parallel job *should* be depends on
## available RAM and cores, and number of loci in <gr>.
##
## AUTHOR: [email protected]
##
## NB: you must install a version of `rtracklayer` which addresses
## issue [import\.bb returning incorrect ranges & results as
## compared with kent tools
## bigBedToBed](https://github.com/lawremi/rtracklayer/issues/59)
## (such as by `install_github('lawremi/rtracklayer')`)
si<-parallel::splitIndices(length(gr),mc.jobs)
mm<-mclapply(si,function(i) {
gri<-gr[i]
m<-import.bb(bbf,selection=BigBedSelection(gri, colnames = c(motifID)))
fo<-findOverlaps(gri,m)
xt<-
xtabs(~.,data=list(locus=i[queryHits(fo)],motif=mcols(m[subjectHits(fo)])[[motifID]]),sparse=TRUE)
xt
})
mm<-
## rBind.fill(mm) ## NOT! It does NOT preserve dimnames and
## returns matrix when given Matrix. Issue reported:
## [[https://github.com/cvarrichio/Matrix.utils/issues/5][rBind.fill
## on a list of Matrix arguably should return Matrix but returns
## matrix #5]]
Reduce(rBind.fill,mm) # which does preserve dimnames and return Matrix
mm
}
NB: If you choose this approach, you MUST install a version of rtracklayer which addresses issue the import.bb returning incorrect ranges & results as compared with kent tools bigBedToBed (such as by install_github('lawremi/rtracklayer'))
