locateVariants failing on intronic variants

[jayoung]

hi there,

I'm experiencing a bug that's the same as an issue that was closed back in November, but it still seems to be a problem.

I am using the release version (1.50.0) - apologies if this is fixed in devel. I cannot right now wrap my head around running the devel version - I used to know how to switch back and forth between versions more easily....

Anyway:

locateVariants() fails to annotate some intronic variants, unless I create and use a cached version of the txdb where I've removed empty elements. It's very similar to the older report, but here's the code I use to illustrate it:

I'm using example data as seen in the vignette

thanks!

Janet

library(TxDb.Hsapiens.UCSC.hg19.knownGene)
library(VariantAnnotation)

## read in test vcf file and sort out seqlevels
fl <- system.file("extdata", "chr22.vcf.gz", package="VariantAnnotation")
vcf <- readVcf(fl, "hg19")
seqlevels(vcf) <- "chr22"

## define txdb
txdb <- TxDb.Hsapiens.UCSC.hg19.knownGene

## get rowRanges, as in the vignette
rd <- rowRanges(vcf)

### use just the first 10 snps
test_snps <- rd[1:10]

## show that they're in an intron:
txdb_introns <- intronsByTranscript(txdb, use.names=TRUE)
findOverlaps(test_snps, txdb_introns)
# Hits object with 10 hits and 0 metadata columns:
#     queryHits subjectHits
# <integer>   <integer>
#     [1]         1       75253
# [2]         2       75253
# [3]         3       75253
# [4]         4       75253
# [5]         5       75253
# [6]         6       75253
# [7]         7       75253
# [8]         8       75253
# [9]         9       75253
# [10]        10       75253

#### locateVariants returns no annotation for these SNPs, whether we use AllVariants() or IntronVariants()
test_loc <- locateVariants(test_snps, txdb, AllVariants())
test_loc
# 
# GRanges object with 0 ranges and 9 metadata columns:
#     seqnames    ranges strand | LOCATION  LOCSTART    LOCEND   QUERYID      TXID
# <Rle> <IRanges>  <Rle> | <factor> <integer> <integer> <integer> <integer>
#     CDSID      GENEID       PRECEDEID        FOLLOWID
# <IntegerList> <character> <CharacterList> <CharacterList>
#     -------
#     seqinfo: no sequences

test_loc_2 <- locateVariants(test_snps, txdb, IntronVariants())
test_loc_2
# GRanges object with 0 ranges and 9 metadata columns:
#     seqnames    ranges strand | LOCATION  LOCSTART    LOCEND   QUERYID      TXID
# <Rle> <IRanges>  <Rle> | <factor> <integer> <integer> <integer> <integer>
#     CDSID      GENEID       PRECEDEID        FOLLOWID
# <IntegerList> <character> <CharacterList> <CharacterList>
#     -------
#     seqinfo: no sequences


#### if we make a cached txdb where we remove empty elements it DOES find the intronic variants
## https://github.com/Bioconductor/VariantAnnotation/issues/76

cache <- new.env(parent = emptyenv())

# CodingVariants
txdb <- TxDb.Hsapiens.UCSC.hg19.knownGene::TxDb.Hsapiens.UCSC.hg19.knownGene
z <- cdsBy(txdb)
z <- z[lengths(z) > 0L]
cache[["cdsbytx"]] <- z

# IntronVariants
z <- intronsByTranscript(txdb)
z <- z[lengths(z) > 0L]
cache[["intbytx"]] <- z

# ThreeUTRVariants 
z <- threeUTRsByTranscript(txdb)
z <- z[lengths(z) > 0L]
cache[["threeUTRbytx"]] <- z

# FiveUTRVariants 
z <- fiveUTRsByTranscript(txdb)
z <- z[lengths(z) > 0L]
cache[["fiveUTRbytx"]] <- z

# IntergenicVariants 
z <- transcriptsBy(txdb, "gene")
z <- z[lengths(z) > 0L]
cache[["txbygene"]]  <- z

# SpliceSiteVariants 
# no need to populate cache further ?

# PromoterVariants 
z <- transcripts(txdb)
z <- z[lengths(z) > 0L]
cache[["tx"]] <- splitAsList(z, seq_len(length(z))) 
names(cache[["tx"]]) <- z$tx_id

test_loc_3 <- locateVariants(test_snps, txdb, AllVariants(), cache = cache)
test_loc_3
# GRanges object with 10 ranges and 9 metadata columns:
#     seqnames    ranges strand | LOCATION  LOCSTART    LOCEND   QUERYID        TXID
# <Rle> <IRanges>  <Rle> | <factor> <integer> <integer> <integer> <character>
#     rs7410291    chr22  50300078      - |   intron     10763     10763         1       75253
# rs147922003    chr22  50300086      - |   intron     10755     10755         2       75253
# rs114143073    chr22  50300101      - |   intron     10740     10740         3       75253
# rs141778433    chr22  50300113      - |   intron     10728     10728         4       75253
# rs182170314    chr22  50300166      - |   intron     10675     10675         5       75253
# rs115145310    chr22  50300187      - |   intron     10654     10654         6       75253
# rs186769856    chr22  50300268      - |   intron     10573     10573         7       75253
# rs77627744    chr22  50300346      - |   intron     10495     10495         8       75253
# rs193230365    chr22  50300423      - |   intron     10418     10418         9       75253
# rs9627788    chr22  50300438      - |   intron     10403     10403        10       75253
# CDSID      GENEID       PRECEDEID        FOLLOWID
# <IntegerList> <character> <CharacterList> <CharacterList>
#     rs7410291                     79087                                
# rs147922003                     79087                                
# rs114143073                     79087                                
# rs141778433                     79087                                
# rs182170314                     79087                                
# rs115145310                     79087                                
# rs186769856                     79087                                
# rs77627744                     79087                                
# rs193230365                     79087                                
# rs9627788                     79087                                
# -------
#     seqinfo: 1 sequence from an unspecified genome; no seqlengths

[jayoung]

sessionInfo():

R version 4.4.1 (2024-06-14)
Platform: x86_64-apple-darwin20
Running under: macOS Sonoma 14.5

Matrix products: default
BLAS:   /System/Library/Frameworks/Accelerate.framework/Versions/A/Frameworks/vecLib.framework/Versions/A/libBLAS.dylib 
LAPACK: /Library/Frameworks/R.framework/Versions/4.4-x86_64/Resources/lib/libRlapack.dylib;  LAPACK version 3.12.0

locale:
[1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8

time zone: America/Los_Angeles
tzcode source: internal

attached base packages:
[1] stats4    stats     graphics  grDevices utils     datasets  methods   base     

other attached packages:
 [1] VariantAnnotation_1.50.0                Rsamtools_2.20.0                       
 [3] Biostrings_2.72.1                       XVector_0.44.0                         
 [5] SummarizedExperiment_1.34.0             MatrixGenerics_1.16.0                  
 [7] matrixStats_1.3.0                       TxDb.Hsapiens.UCSC.hg19.knownGene_3.2.2
 [9] GenomicFeatures_1.56.0                  AnnotationDbi_1.66.0                   
[11] Biobase_2.64.0                          GenomicRanges_1.56.1                   
[13] GenomeInfoDb_1.40.1                     IRanges_2.38.1                         
[15] S4Vectors_0.42.1                        BiocGenerics_0.50.0                    

loaded via a namespace (and not attached):
 [1] SparseArray_1.4.8        bitops_1.0-8             RSQLite_2.3.7           
 [4] lattice_0.22-6           grid_4.4.1               fastmap_1.2.0           
 [7] blob_1.2.4               jsonlite_1.8.8           Matrix_1.7-0            
[10] restfulr_0.0.15          DBI_1.2.3                httr_1.4.7              
[13] BSgenome_1.72.0          UCSC.utils_1.0.0         XML_3.99-0.17           
[16] codetools_0.2-20         abind_1.4-5              cli_3.6.3               
[19] rlang_1.1.4              crayon_1.5.3             bit64_4.0.5             
[22] cachem_1.1.0             DelayedArray_0.30.1      yaml_2.3.10             
[25] S4Arrays_1.4.1           tools_4.4.1              parallel_4.4.1          
[28] BiocParallel_1.38.0      memoise_2.0.1            GenomeInfoDbData_1.2.12 
[31] curl_5.2.1               vctrs_0.6.5              R6_2.5.1                
[34] png_0.1-8                BiocIO_1.14.0            rtracklayer_1.64.0      
[37] zlibbioc_1.50.0          KEGGREST_1.44.1          bit_4.0.5               
[40] pkgconfig_2.0.3          GenomicAlignments_1.40.0 rstudioapi_0.16.0       
[43] rjson_0.2.21             compiler_4.4.1           RCurl_1.98-1.16