Compatibility between `EnsDb` objects and `BSgenome` objects

[mschubert]

I'm trying to use the VariantAnnotation package to annotate a VCF file from nf-core/sarek (UCSC style) using AnnotationHub EnsDb objects (Ensembl style), where I encountered the following issues:

1. The VariantAnnotation package does not accept EnsDb objects. I raised this (Allow 'ANY' txdb, eg. EnsDb objects Bioconductor/VariantAnnotation#74), but it will probably not get fixed; I worked around this by providing my own S4 method in my code
2. The seqlevelsStyles() do not match

For point (2), I can change the EnsDb style to UCSC:

ens106 = AnnotationHub::AnnotationHub()[["AH100643"]]
seqlevelsStyle(ens106) = "UCSC"

and then either load UCSC-style genome or change the Ensembl-style genome to UCSC:

asm = BSgenome.Hsapiens.NCBI.GRCh38::BSgenome.Hsapiens.NCBI.GRCh38
seqlevelsStyle(asm) = "UCSC"

asm = BSgenome.Hsapiens.UCSC.hg38::BSgenome.Hsapiens.UCSC.hg38

However, this does not work, because the genome() of the EnsDb object is GRCh38, while the one of the assembly is hg38, raising an assertion error in VariantAnnotation. So this needs an additional line changing the internal state of the S4 genome object (we can't change the EnsDb object):

genome(ens106)[] = "hg38"
# Error: unable to find an inherited method for function 'seqinfo<-' for signature 'x = "EnsDb"'

asm@seqinfo@genome[] = "GRCh38" # works, but is messing with internals
# this was previously the default, but they explicitly changed it

I'm not sure what a good solution is here. It seems to be that a check if they genomes are identical (as performed by VariantAnnotation) is reasonable. I'm raising this issue more to document it rather than suggesting a change in ensembldb.

[jorainer]

Thanks for reporting. Is there no way to create or load a BSgenome or other supported genome sequence object from AnnotationHub? There should also be the twobit files from Ensembl accessible through AnnotationHub - maybe these can be used? The advantage would be that you could ensure to use matching EnsDb and genome sequence from the same Ensembl release. I'm just a bit worried that hg38 is not exactly identical to the GRCh38 version used by Ensembl...

[mschubert]

This is still an unresolved issue for us, so I checked the twobit files for the release now.

The problem with those is that they do not provide methods to interact with Bioconductor core infrastructure:

asm = AnnotationHub::AnnotationHub()[["AH103844"]] # 2bit
seqlevelsStyle(asm) = "UCSC"
# Error: unable to find an inherited method for function 'seqinfo<-' for signature 'x = "TwoBitFile"'
seqlevelsInUse(asm)
# Error: unable to find an inherited method for function 'seqlevelsInUse' for signature 'x = "TwoBitFile"'

Regarding whether this is the exact same genome, the BSgenome packages switch between GRCh38 and hg38 depending on the seqlevelsStyle parameter.

[jorainer]

I added now a genome<- method to manually overwrite the Genome build info. Can you please install the latest devel version BiocManager::install("jorainer/ensembldb", ref = "devel") to see if this works for you @mschubert ?

Basically,

ens106 = AnnotationHub::AnnotationHub()[["AH100643"]]
seqlevelsStyle(ens106) = "UCSC"
genome(ens106) <- "hg38"

should do the trick (note that genome<- expects a character of length 1!).