Question: Problems when exporting GRanges file as bigwig after hg38 to hg19 lift over using rtracklayer and annotation hub
0
Entering edit mode

I'm trying to use rtracklayer and AnnotationHub in R to lift over a .bw file I have that is aligned to hg38. I need the data to be aligned to hg19.

By following this vingette, I think I have successfully done this using the R packages mentioned above with the code below.

library(rtracklayer)
library(AnnotationHub)

mphg_bw_hg38 <- import(con = 'CEPBb_macrophage.bw', format = 'bigWig')

ahub <- AnnotationHub()
ahub.chain <- subset(ahub, rdataclass == "ChainFile" & species == "Homo sapiens")
query(ahub.chain, c("hg38", "hg19"))
chain <- ahub.chain[ahub.chain$title == "hg38ToHg19.over.chain.gz"]
chain <- chain[[1]]

mphg_bw_hg19 <- liftOver(mphg_bw_hg38, chain)

My issue is that I'm having difficulty converting the file created during the liftover back to .bw for visualisation in IGV.

The output of the liftover is a compressed GRanges list. I tried to export to bigwig this using:

 export(mphg_bw_hg19, con = 'CEPBb_macrophage_hg19.bw', format = 'bigWig')

But I got the following: Error in .local(object, con, format, ...) : 'object' must have names

So I tried using unlist to convert the file to GRanges then export ...

mphg_bw_hg19_a <- unlist(mphg_bw_hg19)
export(mphg_bw_hg19_a, con = 'CEPBb_macrophage_hg19.bw', format = 'bigWig')

But this threw: Error in .local(object, con, format, ...) : Unable to determine seqlengths; either specify 'seqlengths' or specify a genome on 'object' that is known to BSgenome or UCSC

I assume I need to assign hg19 to this GRanges file somehow and/or specify the seqlengths of all the chromosomes in the file - is there a simple way to do this? Not an expert on GRanges!

ADD COMMENTlink 9 months ago camerond • 70 • updated 9 months ago benformatics • 870
1
Entering edit mode

You need to get the chromosome information for the human hg19 build

My suggestion:

libary(TxDb.Hsapiens.UCSC.hg19.knownGene)
txdb <- TxDb.Hsapiens.UCSC.hg19.knownGene
seqinfo(bw_hg19_a) <- seqinfo(txdb)[seqnames(seqinfo(bw_hg19_a))]
export(mphg_bw_hg19_a, con = 'CEPBb_macrophage_hg19.bw', format = 'bigWig')

However if you have overlapping ranges (following liftOver) you will get a new error. That's because you probably shouldn't be liftOver'ing bigWig files...

ADD COMMENTlink 9 months ago benformatics • 870
Entering edit mode
0

Thanks for the suggestion, but your right that threw an overlapping error. Error in FUN(extractROWS(unlisted_X, IRanges(X_elt_start[i], X_elt_end[i])), : BigWig ranges cannot overlap

The vignettes for the packages I used suggest that it is possible to lift over using a bigwig file ...

ADD REPLYlink 9 months ago
camerond
• 70

Login before adding your answer.

Powered by the version 1.8