Working with File GT Format Databases

April 25, 2020

Load Packages

library(VariantAnnotation)
library(tidyverse)
library(dbplyr)

Connect to the toy database

con <- DBI::dbConnect(RSQLite::SQLite(), 'toy-data/file-genos/file-gt.db')

The main complicating factor of a File GT database is how genotypes must be read in.

Please read through the Working with Table GT Format Databases tutorial first if you are not already familiar with dbplyr

Let’s look at Example 1 from Working with Table GT Format Databases in the context of a File GT database. This is what the variant_info table looks like in a File GT database.

tbl(con, 'variant_impact') %>%
  filter(consequence == "regulatory_region_variant")  %>%
  distinct(variant_id) %>%
  inner_join(tbl(con, 'variant_info'), by = 'variant_id') %>%
  filter(af > 0.01 & af < 0.99) %>%
  select(variant_id:alt, geno)

## # Source:   lazy query [?? x 7]
## # Database: sqlite 3.30.1
## #   [/wdata/tkoomar/VCFdbR/toy-data/file-genos/file-gt.db]
##   variant_id chr     start     end ref   alt   geno                             
##        <int> <chr>   <int>   <int> <chr> <chr> <chr>                            
## 1         11 22     1.61e7  1.61e7 A     <CN2> /wdata/tkoomar/VCFdbR/toy-data/f…
## 2         12 22     1.61e7  1.61e7 A     <CN3> /wdata/tkoomar/VCFdbR/toy-data/f…
## 3        713 22     1.61e7  1.61e7 C     G     /wdata/tkoomar/VCFdbR/toy-data/f…

The most efficient way to read in the genotypes for a set of variants of interest is to lapply()/map() over the geno column of the variant_impact table:

tbl(con, 'variant_impact') %>%
  filter(consequence == "regulatory_region_variant")  %>%
  distinct(variant_id) %>%
  inner_join(tbl(con, 'variant_info'), by = 'variant_id') %>%
  filter(af > 0.01 & af < 0.99) %>%
  select(geno) %>%
  collect() %>%
  mutate(geno = map(geno, read_rds)) %>%
  unnest(geno)

## # A tibble: 7,512 x 4
##    variant_id sample     gt gt_raw
##         <int> <chr>   <dbl> <chr> 
##  1         11 HG00096     0 0|0   
##  2         11 HG00097     0 0|0   
##  3         11 HG00099     0 0|0   
##  4         11 HG00100     0 0|0   
##  5         11 HG00101     0 0|0   
##  6         11 HG00102     0 0|0   
##  7         11 HG00103     0 0|0   
##  8         11 HG00105     0 0|0   
##  9         11 HG00106     0 0|0   
## 10         11 HG00107     0 0|0   
## # … with 7,502 more rows

Note that in the case of a File GT database, you will always need to collect() your genos/variant_ids before iterating over them. This can add memory overhead compared to a Table GT database.

set.seed(1099015)
phenos <- tbl(con, 'samples') %>%
  collect() %>%
  transmute(sample = value, 
            phenotype = sample(0:1, n(), replace = TRUE))

library(microbenchmark)

file_gt_table_glm <- microbenchmark({
  tbl(con, 'variant_impact') %>%
  filter(gene != "") %>%
  select(gene, variant_id) %>%
  distinct() %>%
  inner_join(tbl(con, 'variant_info')) %>%
  select(gene, geno) %>%
  collect() %>%
  group_by(gene) %>%
  summarize(
    model = list(
      map(geno, read_rds) %>%
        bind_rows() %>%
        pivot_wider(id_cols = sample, names_from = variant_id, values_from = gt) %>%
        inner_join(phenos, by = 'sample') %>%
        {glm(phenotype ~ . - sample, family = 'binomial', data = .)}
    ), 
    model_summary = map(model, broom::glance)
  )
}, times = 10)

summary(file_gt_table_glm$time/1e9)

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   6.575   6.820   6.850   7.071   7.076   8.543

Note that this is somewhere in-between the speed of the two examples in Working with Table GT Format Databases.

Always disconnect when done working with the database!

DBI::dbDisconnect(con)