Repository navigation
Working with Table GT Format Databases
April 25, 2020
library(VariantAnnotation)
library(tidyverse)
library(dbplyr)con <- DBI::dbConnect(RSQLite::SQLite(), 'toy-data/table-genos/table-gt.db')filter() is primary tool for identifying variants of interest. Let’s
say we want to investigate variants in regulatory regions:
tbl(con, 'variant_impact') %>%
filter(consequence == "regulatory_region_variant") ## # Source: lazy query [?? x 53]
## # Database: sqlite 3.30.1
## # [/wdata/tkoomar/VCFdbR/toy-data/table-genos/table-gt.db]
## variant_id consequence impact symbol gene feature_type feature biotype exon
## <int> <chr> <chr> <chr> <chr> <chr> <chr> <chr> <chr>
## 1 10 regulatory… MODIF… "" "" RegulatoryF… ENSR00… "" ""
## 2 11 regulatory… MODIF… "" "" RegulatoryF… ENSR00… "" ""
## 3 12 regulatory… MODIF… "" "" RegulatoryF… ENSR00… "" ""
## 4 13 regulatory… MODIF… "" "" RegulatoryF… ENSR00… "" ""
## 5 88 regulatory… MODIF… "" "" RegulatoryF… ENSR00… "CTCF_… ""
## 6 89 regulatory… MODIF… "" "" RegulatoryF… ENSR00… "CTCF_… ""
## 7 90 regulatory… MODIF… "" "" RegulatoryF… ENSR00… "CTCF_… ""
## 8 91 regulatory… MODIF… "" "" RegulatoryF… ENSR00… "CTCF_… ""
## 9 92 regulatory… MODIF… "" "" RegulatoryF… ENSR00… "CTCF_… ""
## 10 93 regulatory… MODIF… "" "" RegulatoryF… ENSR00… "CTCF_… ""
## # … with more rows, and 44 more variables: intron <chr>, hgvsc <chr>,
## # hgvsp <chr>, cdna_position <chr>, cds_position <chr>,
## # protein_position <chr>, amino_acids <chr>, codons <chr>,
## # existing_variation <chr>, allele_num <chr>, distance <chr>, strand <chr>,
## # flags <chr>, variant_class <chr>, symbol_source <chr>, hgnc_id <chr>,
## # canonical <chr>, tsl <chr>, appris <chr>, ccds <chr>, ensp <chr>,
## # swissprot <chr>, trembl <chr>, uniparc <chr>, refseq_match <chr>,
## # source <chr>, given_ref <chr>, used_ref <chr>, bam_edit <chr>,
## # gene_pheno <chr>, sift <chr>, polyphen <chr>, domains <chr>,
## # hgvs_offset <chr>, af <chr>, max_af_pops <chr>, motif_name <chr>,
## # motif_pos <chr>, high_inf_pos <chr>, motif_score_change <chr>,
## # is_lof <int>, is_splicing <int>, is_exonic <int>, is_intronic <int>
We need to get a non-redundant list of the variant IDs
tbl(con, 'variant_impact') %>%
filter(consequence == "regulatory_region_variant") %>%
distinct(variant_id)## # Source: lazy query [?? x 1]
## # Database: sqlite 3.30.1
## # [/wdata/tkoomar/VCFdbR/toy-data/table-genos/table-gt.db]
## variant_id
## <int>
## 1 10
## 2 11
## 3 12
## 4 13
## 5 88
## 6 89
## 7 90
## 8 91
## 9 92
## 10 93
## # … with more rows
Now, let’s say we want to make sure the variants are relatively common
in our cohort (minor allele is more than 1%). We can do this by
cross-referencing these variant IDs with the variant_info table via an
inner_join():
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)## # Source: lazy query [?? x 34]
## # Database: sqlite 3.30.1
## # [/wdata/tkoomar/VCFdbR/toy-data/table-genos/table-gt.db]
## variant_id chr start end ref alt qual filter ciend cipos cs
## <int> <chr> <int> <int> <chr> <chr> <dbl> <chr> <chr> <chr> <chr>
## 1 11 22 1.61e7 1.61e7 A <CN2> 100 PASS c(NA… c(NA… DUP_…
## 2 12 22 1.61e7 1.61e7 A <CN3> 100 PASS c(NA… c(NA… DUP_…
## 3 713 22 1.61e7 1.61e7 C G 100 PASS c(NA… c(NA… <NA>
## # … with 23 more variables: imprecise <int>, mc <chr>, meinfo <chr>,
## # mend <int>, mlen <int>, mstart <int>, svlen <chr>, svtype <chr>, tsd <chr>,
## # ac <int>, af <dbl>, ns <int>, an <int>, eas_af <dbl>, eur_af <dbl>,
## # afr_af <dbl>, amr_af <dbl>, sas_af <dbl>, dp <int>, aa <chr>, vt <chr>,
## # ex_target <int>, multi_allelic <int>
In this case, from population-specific alternate allele frequencies:
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(chr, start, ref, alt, ends_with("_af")) %>%
collect() %>%
mutate(id = str_c(chr, start, ref, alt, sep = ":")) %>%
pivot_longer(cols = ends_with("_af"), names_to = "population") %>%
ggplot(aes(x = population, y = value, fill = population)) +
geom_col() +
facet_wrap(~id) +
theme_minimal() +
scale_fill_brewer(palette = 'Set2', guide = 'none') +
scale_x_discrete(guide = guide_axis(n.dodge = 2))
This is again achieved with the help of inner_join().
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) %>%
inner_join(tbl(con, 'variant_geno'), by = 'variant_id')## # Source: lazy query [?? x 4]
## # Database: sqlite 3.30.1
## # [/wdata/tkoomar/VCFdbR/toy-data/table-genos/table-gt.db]
## 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 more rows
Note that is important to first select(variant_id) so that you do
not create a massive table including all of the variant_info columns.
If you want to merge genotypes with annotations, it is recommended to
pivot_wider() first:
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) %>%
inner_join(tbl(con, 'variant_geno'), by = 'variant_id') %>%
collect() %>%
pivot_wider(id_cols = variant_id,
names_from = sample,
values_from = gt)## # A tibble: 3 x 2,505
## variant_id HG00096 HG00097 HG00099 HG00100 HG00101 HG00102 HG00103 HG00105
## <int> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 11 0 0 0 0 0 0 0 0
## 2 12 1 0 0 0 0 0 0 0
## 3 713 0 0 0 0 0 0 0 0
## # … with 2,496 more variables: HG00106 <dbl>, HG00107 <dbl>, HG00108 <dbl>,
## # HG00109 <dbl>, HG00110 <dbl>, HG00111 <dbl>, HG00112 <dbl>, HG00113 <dbl>,
## # HG00114 <dbl>, HG00115 <dbl>, HG00116 <dbl>, HG00117 <dbl>, HG00118 <dbl>,
## # HG00119 <dbl>, HG00120 <dbl>, HG00121 <dbl>, HG00122 <dbl>, HG00123 <dbl>,
## # HG00125 <dbl>, HG00126 <dbl>, HG00127 <dbl>, HG00128 <dbl>, HG00129 <dbl>,
## # HG00130 <dbl>, HG00131 <dbl>, HG00132 <dbl>, HG00133 <dbl>, HG00136 <dbl>,
## # HG00137 <dbl>, HG00138 <dbl>, HG00139 <dbl>, HG00140 <dbl>, HG00141 <dbl>,
## # HG00142 <dbl>, HG00143 <dbl>, HG00145 <dbl>, HG00146 <dbl>, HG00148 <dbl>,
## # HG00149 <dbl>, HG00150 <dbl>, HG00151 <dbl>, HG00154 <dbl>, HG00155 <dbl>,
## # HG00157 <dbl>, HG00158 <dbl>, HG00159 <dbl>, HG00160 <dbl>, HG00171 <dbl>,
## # HG00173 <dbl>, HG00174 <dbl>, HG00176 <dbl>, HG00177 <dbl>, HG00178 <dbl>,
## # HG00179 <dbl>, HG00180 <dbl>, HG00181 <dbl>, HG00182 <dbl>, HG00183 <dbl>,
## # HG00185 <dbl>, HG00186 <dbl>, HG00187 <dbl>, HG00188 <dbl>, HG00189 <dbl>,
## # HG00190 <dbl>, HG00231 <dbl>, HG00232 <dbl>, HG00233 <dbl>, HG00234 <dbl>,
## # HG00235 <dbl>, HG00236 <dbl>, HG00237 <dbl>, HG00238 <dbl>, HG00239 <dbl>,
## # HG00240 <dbl>, HG00242 <dbl>, HG00243 <dbl>, HG00244 <dbl>, HG00245 <dbl>,
## # HG00246 <dbl>, HG00250 <dbl>, HG00251 <dbl>, HG00252 <dbl>, HG00253 <dbl>,
## # HG00254 <dbl>, HG00255 <dbl>, HG00256 <dbl>, HG00257 <dbl>, HG00258 <dbl>,
## # HG00259 <dbl>, HG00260 <dbl>, HG00261 <dbl>, HG00262 <dbl>, HG00263 <dbl>,
## # HG00264 <dbl>, HG00265 <dbl>, HG00266 <dbl>, HG00267 <dbl>, HG00268 <dbl>,
## # HG00269 <dbl>, HG00271 <dbl>, …
To do this, it is helpful to first identify all of the variants in each
gene on the variant_impact table:
genes <- tbl(con, 'variant_impact') %>%
filter(gene != "") %>%
select(gene, variant_id) %>%
distinct() %>%
collect()
head(genes)## # A tibble: 6 x 2
## gene variant_id
## <chr> <int>
## 1 ENSG00000215270 819
## 2 ENSG00000215270 820
## 3 ENSG00000215270 821
## 4 ENSG00000215270 822
## 5 ENSG00000215270 823
## 6 ENSG00000215270 824
Let’s generate some random phenotypes to do a simple test against. . .
set.seed(1099015)
phenos <- tbl(con, 'samples') %>%
collect() %>%
transmute(sample = value,
phenotype = sample(0:1, n(), replace = TRUE))
head(phenos)## # A tibble: 6 x 2
## sample phenotype
## <chr> <int>
## 1 HG00096 1
## 2 HG00097 0
## 3 HG00099 1
## 4 HG00100 1
## 5 HG00101 1
## 6 HG00102 1
We can then pull in the genotypes for the variants in each table
gene_genos <- genes %>%
nest(variants = variant_id) %>%
mutate(
genos = map(variants, function(.x){
v <- .x$variant_id
tbl(con, 'variant_geno') %>%
filter(variant_id %in% v) %>%
collect() %>%
pivot_wider(id_cols = sample, names_from = variant_id, values_from = gt)
})
)
gene_genos## # A tibble: 5 x 3
## gene variants genos
## <chr> <list> <list>
## 1 ENSG00000215270 <tibble [84 × 1]> <tibble [2,504 × 85]>
## 2 ENSG00000223875 <tibble [206 × 1]> <tibble [2,504 × 207]>
## 3 ENSG00000229286 <tibble [56 × 1]> <tibble [2,504 × 57]>
## 4 ENSG00000233866 <tibble [281 × 1]> <tibble [2,504 × 282]>
## 5 ENSG00000235265 <tibble [40 × 1]> <tibble [2,504 × 41]>
gene_genos$genos[[1]][1:3, 1:5]## # A tibble: 3 x 5
## sample `819` `820` `821` `822`
## <chr> <dbl> <dbl> <dbl> <dbl>
## 1 HG00096 0 0 0 0
## 2 HG00097 0 0 0 0
## 3 HG00099 0 0 0 0
Pulling the genotypes in from the database (with collect()) could
quickly fill up memory, however. It is more efficient to run whatever
operation you want on the genotypes as part of the map() function and
just capture the output.
In this case, let’s do a simple logistic regression on our random phenotype.
genes_glm <- genes %>%
group_by(gene) %>%
transmute(v = variant_id) %>%
summarize(
model = list(
tbl(con, 'variant_geno') %>%
filter(variant_id %in% v) %>%
collect() %>%
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)
)## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
genes_glm %>%
unnest( model_summary)## # A tibble: 5 x 9
## gene model null.deviance df.null logLik AIC BIC deviance df.residual
## <chr> <lis> <dbl> <int> <dbl> <dbl> <dbl> <dbl> <int>
## 1 ENSG00000… <glm> 3471. 2503 -1681. 3527. 4011. 3361. 2421
## 2 ENSG00000… <glm> 3471. 2503 -1613. 3619. 4767. 3225. 2307
## 3 ENSG00000… <glm> 3471. 2503 -1702. 3515. 3841. 3403. 2448
## 4 ENSG00000… <glm> 3471. 2503 -1578. 3688. 5237. 3156. 2238
## 5 ENSG00000… <glm> 3471. 2503 -1710. 3501. 3734. 3421. 2464
You may have noticed in the above code we had to change the name of the
variant_id column to v:
transmute(v = variant_id)This was largely to remove ambiguity, for when we later filter the genotypes to only include those in the gene of interest with:
tbl(con, 'variant_geno') %>%
filter(variant_id %in% v)However, why couldn’t we use the trick from Example 1 and use
inner_join(), rather than filter(), like so:
inner_join(., tbl(con, 'variant_geno') )One important thing to keep in mind with the above – and all operations
– on SQLite databases is if the data you are working with is local (in
memory) or remote (in the database). A lot of operations automatically
occur without actually pulling data into memory. This this particular
case, the genes object was in memory
Here’s a simple example of how this can go wrong:
common_alt_remote <- tbl(con, 'variant_info') %>%
filter(af > 0.5)
common_alt_local <- tbl(con, 'variant_info') %>%
filter(af > 0.5) %>%
collect()By using, collect() we have pulled common_alt_local into memory. You
can see that it is slightly larger than common_alt_remote. If it
selected more rows, this difference would be substantial.
object.size(common_alt_local)## 11448 bytes
object.size(common_alt_remote)## 9272 bytes
Also note that, we get an error if we try to join common_alt_local
with variant_impact:
common_alt_local %>%
inner_join(tbl(con, 'variant_impact'))## Joining, by = c("variant_id", "af")
## Error: `x` and `y` must share the same src, set `copy` = TRUE (may be slow)
But not if we join with common_alt_remote
common_alt_remote %>%
inner_join(tbl(con, 'variant_impact'), by = 'variant_id')## # Source: lazy query [?? x 86]
## # Database: sqlite 3.30.1
## # [/wdata/tkoomar/VCFdbR/toy-data/table-genos/table-gt.db]
## variant_id chr start end ref alt qual filter ciend cipos cs
## <int> <chr> <int> <int> <chr> <chr> <dbl> <chr> <chr> <chr> <chr>
## 1 390 22 1.61e7 1.61e7 T C 100 PASS c(NA… c(NA… <NA>
## 2 391 22 1.61e7 1.61e7 C T 100 PASS c(NA… c(NA… <NA>
## 3 392 22 1.61e7 1.61e7 T C 100 PASS c(NA… c(NA… <NA>
## 4 394 22 1.61e7 1.61e7 T G 100 PASS c(NA… c(NA… <NA>
## 5 736 22 1.61e7 1.61e7 C G 100 PASS c(NA… c(NA… <NA>
## 6 736 22 1.61e7 1.61e7 C G 100 PASS c(NA… c(NA… <NA>
## 7 338 22 1.61e7 1.61e7 A G 100 PASS c(NA… c(NA… <NA>
## 8 554 22 1.61e7 1.61e7 C T 100 PASS c(NA… c(NA… <NA>
## 9 547 22 1.61e7 1.61e7 T C 100 PASS c(NA… c(NA… <NA>
## 10 958 22 1.61e7 1.61e7 A G 100 PASS c(NA… c(NA… <NA>
## # … with more rows, and 75 more variables: imprecise <int>, mc <chr>,
## # meinfo <chr>, mend <int>, mlen <int>, mstart <int>, svlen <chr>,
## # svtype <chr>, tsd <chr>, ac <int>, af.x <dbl>, ns <int>, an <int>,
## # eas_af <dbl>, eur_af <dbl>, afr_af <dbl>, amr_af <dbl>, sas_af <dbl>,
## # dp <int>, aa <chr>, vt <chr>, ex_target <int>, multi_allelic <int>,
## # consequence <chr>, impact <chr>, symbol <chr>, gene <chr>,
## # feature_type <chr>, feature <chr>, biotype <chr>, exon <chr>, intron <chr>,
## # hgvsc <chr>, hgvsp <chr>, cdna_position <chr>, cds_position <chr>,
## # protein_position <chr>, amino_acids <chr>, codons <chr>,
## # existing_variation <chr>, allele_num <chr>, distance <chr>, strand <chr>,
## # flags <chr>, variant_class <chr>, symbol_source <chr>, hgnc_id <chr>,
## # canonical <chr>, tsl <chr>, appris <chr>, ccds <chr>, ensp <chr>,
## # swissprot <chr>, trembl <chr>, uniparc <chr>, refseq_match <chr>,
## # source <chr>, given_ref <chr>, used_ref <chr>, bam_edit <chr>,
## # gene_pheno <chr>, sift <chr>, polyphen <chr>, domains <chr>,
## # hgvs_offset <chr>, af.y <chr>, max_af_pops <chr>, motif_name <chr>,
## # motif_pos <chr>, high_inf_pos <chr>, motif_score_change <chr>,
## # is_lof <int>, is_splicing <int>, is_exonic <int>, is_intronic <int>
If you’re frequently having trouble with such errors, read over the documentation of dbplyr in detail.
How might we do Example 2 keeping data out of memory longer?
genes_remote <- tbl(con, 'variant_impact') %>%
filter(gene != "") %>%
select(gene, variant_id) %>%
distinct()
genes_remote_glm <- genes_remote %>%
select(gene) %>%
distinct() %>%
collect() %>%
mutate(
model = map(gene, function(.x){
genes_remote %>%
filter(gene %in% .x) %>%
inner_join(tbl(con, 'variant_geno'), by = 'variant_id') %>%
collect() %>%
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)
)## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
genes_remote_glm %>%
unnest(model_summary)## # A tibble: 5 x 9
## gene model null.deviance df.null logLik AIC BIC deviance df.residual
## <chr> <lis> <dbl> <int> <dbl> <dbl> <dbl> <dbl> <int>
## 1 ENSG00000… <glm> 3471. 2503 -1681. 3527. 4011. 3361. 2421
## 2 ENSG00000… <glm> 3471. 2503 -1613. 3619. 4767. 3225. 2307
## 3 ENSG00000… <glm> 3471. 2503 -1702. 3515. 3841. 3403. 2448
## 4 ENSG00000… <glm> 3471. 2503 -1578. 3688. 5237. 3156. 2238
## 5 ENSG00000… <glm> 3471. 2503 -1710. 3501. 3734. 3421. 2464
And now for the ultimate question: which way is faster?
library(microbenchmark)
collect_sooner <- microbenchmark({
tbl(con, 'variant_impact') %>%
filter(gene != "") %>%
select(gene, variant_id) %>%
distinct() %>%
collect() %>%
group_by(gene) %>%
transmute(v = variant_id) %>%
summarize(
model = list(
tbl(con, 'variant_geno') %>%
filter(variant_id %in% v) %>%
collect() %>%
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)
collect_later <- microbenchmark({
genes_remote <- tbl(con, 'variant_impact') %>%
filter(gene != "") %>%
select(gene, variant_id) %>%
distinct()
genes_remote_glm <- genes_remote %>%
select(gene) %>%
distinct() %>%
collect() %>%
mutate(
model = map(gene, function(.x){
genes_remote %>%
filter(gene %in% .x) %>%
inner_join(tbl(con, 'variant_geno'), by = 'variant_id') %>%
collect() %>%
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)Collecting Sooner
summary(collect_sooner$time/1e9)## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 5.905 6.933 7.916 7.787 8.689 9.235
Collecting Latter
summary(collect_later$time/1e9)## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 5.955 6.376 6.981 7.220 7.943 8.805
Clearly collecting later is faster, and with more data the effect would
likely be larger. In general, if an operation seems like it is taking
longer than it should, the best solution (other than checking that you
have an index for each column that is used in a filter()) is to find a
way to postpone reading data into memory for as long as possible.
DBI::dbDisconnect(con)