2
$\begingroup$

I am now trying to find intersect genes of two files: one is .bigbed file (very large, 117GB), and the other is a .bed file.

At first, I tried to convert bigbed file to bed file and then applied bedtools. I didn't run the bigbedtobed program on the server because I am not allowed to install the bigbedtobed program on our server. Thus I planned to convert the .bb file on my own computer and then upload the converted .bed file to the server, and then apply bedtools. Firstly I tried to convert some smaller .bb data sets (< 10 GB) to .bed and it worked slowly but well on my computer. However, when I tried the 117GB data there was an error:

Invalid argument, Error reading 9188986873 bytes.

My code is :

bigBedToBed JASPAR2022_hg19.bb hg1 9.bed

Thank you in advance!


Now the bug is solved, and I got many other good suggestions. Below you can find many methods to work with the bigbed file, directly working with it or converting it to bed file.

$\endgroup$
0

4 Answers 4

1
$\begingroup$

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: malcolm_cook@stowers.org
  ##
  ## 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'))

$\endgroup$
3
  • $\begingroup$ Yes Malcook, this is exactly what I want at first! Thank you so much! $\endgroup$ Feb 21, 2022 at 9:28
  • $\begingroup$ glad to help - would you consider to change the accepted solution to mine? $\endgroup$
    – malcook
    Mar 29, 2022 at 7:26
  • $\begingroup$ This is the accepted answer for sure @BinhuanSun $\endgroup$
    – M__
    Mar 30, 2022 at 9:09
1
$\begingroup$

No matter how you approach this problem 110+ GB is big:

Your options are:

  1. Break up the file and perform the intersection on a chunk by chunk basis
  2. Ask your system administrator to install the program
  3. Load it into Docker, your systems administrator will understand this. Installing docker is a single line of code (which I forget). Docker is a container which protects the system from the program and no systems administrator in this day and age would refuse a request for Docker installation.

Your.systrm administrator might be concerned about the file sizes involved, but I think a better dialogue with them is the right way forward here.

$\endgroup$
8
  • $\begingroup$ the first option is maybe a good one, you can select individual chromosomes at a time with bigBedToBed $\endgroup$
    – Colin D
    Feb 13, 2022 at 7:43
  • $\begingroup$ Thank you for your suggestions! I will try the first one on my computer first :) $\endgroup$ Feb 13, 2022 at 7:49
  • $\begingroup$ @ColinD Yes, I will try. Thank you :) $\endgroup$ Feb 13, 2022 at 7:49
  • $\begingroup$ Hello M, because today is Sunday so I first tried the first method you suggested. I tried to break up the file using argument -Chrom of the bigBedtobed program: bigBedToBed -chrom=chr1 JASPAR2022_hg19.bb test.bed. However, I still got the same errorInvalid argument, Error reading 9188986873 bytes. I then tried a relatively smaller data (66GB) but got the same error. It seems that on my computer the bigBedtobed program only works on data that is smaller than 10 GB. I still have over 600 GB of drive memory, so I think this may due to the lack of RAM? $\endgroup$ Feb 13, 2022 at 15:42
  • 1
    $\begingroup$ Yes M. I think you are right, and today I asked my supervisor to help me apply for an account on another server that installed those UCSC programs, and I hope it can work on the server. Thank you again for the patient reply :) $\endgroup$ Feb 15, 2022 at 10:09
1
$\begingroup$

The question I asked a couple of days ago has been solved. Here I want to make a summary of this and hope it can help people who also encounter such problems later.

First, about the way to process the bigbed file. I didn't find a very good way to work directly with bigbed file at that time, so I still firstly converted the bigbed file to the bed file using bigBedTobed program and then apply bedtools. The bigBedTobed program can be downloaded through this link: UCSC programs, and it also contains many useful programs from UCSC genome browser. I finally converted my data (117 GB bigbed file) on the server, and it took ~2.5h to convert it to bed file. After conversion, the converted bed file is about 550 GB. And it took ~4h to intersect with a 4 Mb bed file using bedtools.

Second, the bug that I encountered when I tried to convert the bigbed file on my own laptop. The code I used is bigBedToBed JASPAR2022_hg19.bb hg1 9.bed and it has error Invalid argument, Error reading 9188986873 bytes. I reported this error to the UCSC genome browser group and they said it is a MacOS specific bug, reaching a default data limit. They have solved it in their development environment but it still needs some time to reach the public site. The current solution for this is to use a URL to the file instead of a local file allows the program to read appropriately sized pieces of data at a time(-udcDir= in bigBedToBed).

Finally, this is my first time asking questions here and I really appreciate those good suggestions from @M__ and @ColinD.


Renew: Methods from @malcook can work directly with the bigbed file! For a very large bigbed file, it is indeed a more efficient way.

$\endgroup$
2
  • $\begingroup$ Hi nice to hear from you @BinhuanSun. Cool glad its worked out and thanks for the feedback. $\endgroup$
    – M__
    Feb 19, 2022 at 21:05
  • $\begingroup$ MacOSX bug .. thanks thats interesting, I encountered something similar. They are nice machines, but it ain't a cool way to do things IMO. $\endgroup$
    – M__
    Feb 19, 2022 at 21:18
1
$\begingroup$

My 2p, a little late and maybe not too useful:

I didn't run the bigbedtobed program on the server because I am not allowed to install the bigbedtobed program on our server.

Presumably, you cannot install in /usr/local/bin or some other root directory because you don't have admin rights. This is the common setup on shared servers. However, many programs can be installed in the user's home directory without admin rights. Even better, install conda and mamba first (no need to be admin) and use that to install other programs.

The problem with docker is that some system administrators are concerned about security and refuse to install it (and docker does need admin rights). Whether such concerns are founded, I don't know...

I tried to convert bigbed file to bed file and then applied bedtools, but I found it was very slow and couldn't make it.

bedtools intersect supports the -sorted options that should speed things up quite a bit - Have you applied it? Your huge bigbed files should be already sorted after you convert to bed.

$\endgroup$
1
  • $\begingroup$ Hi Dariober, yes later I found I could install programs in my home directory. Sorry for this super easy problem because I am not familiar with the Linux server. I just checked my converted bed file, and it is indeed sorted, so I think it is just because the converted bed file is too big (560GB). $\endgroup$ Feb 21, 2022 at 9:35

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.