its me (yet again)...
I am currently using AnnotationHub to get the latest ens_db object for GRCh38. I noticed, that for some loci, I get the message "No genes to plot". However, I know for a fact that there are genes in the region. I am also able to manually search the ens_db object and list the genes that should be plotted.
This is code to reproduce the issue:
ah <- AnnotationHub()
# this is the latest ensembl annotation
# query(ah, c("EnsDb", "v109"))
ensDB_v109 <- ah[["AH109606"]]
x <- fread("example.tsv")
chr=1
start = 247416156 - 110e3 # these are different from the coloc_susie script because we use hg38 here
stop = 247449108 + 110e3
tmp_loc <- locus(
data =x,
seqname = "chr1",
chrom = "Chromosome_phe",
xrange = c(start, stop),
flank = 0,
p = "Nominal_pvalue",
pos = "Pos_variant",
ens_db = ensDB_v109,
labs = "gene_snp")
locus_plot(tmp_loc) # -> No genes to plot
# Manually look up the ensDB object:
genes(ensDB_v109) %>% as.data.frame() %>%
filter(start > 247416156 - 110e3) %>%
filter(start < 247449108 + 110e3) %>% View() # shows the genes that I expect to be in the region.
The issue is the "chr" prefix in your chromosome name. I have added a line of code to locus() which removes "chr" from seqname. Try now, it works for me.
Hey Myles,
its me (yet again)... I am currently using
AnnotationHub
to get the latestens_db
object for GRCh38. I noticed, that for some loci, I get the message "No genes to plot". However, I know for a fact that there are genes in the region. I am also able to manually search theens_db
object and list the genes that should be plotted. This is code to reproduce the issue:The example data can be downloaded here: https://we.tl/t-wgfpRz9aHd
As always, any insights are much appreciated!
Cheers, Nick