Repository navigation
Schema of Table GT Format Databases
The only package that is totally necessary to load in order to work with
the database is dbplyr, but the other tidyverse packages provide
significant additional utility.
library(tidyverse)
library(dbplyr)con <- DBI::dbConnect(RSQLite::SQLite(), 'toy-data/table-genos/table-gt.db')src_dbi(con)## src: sqlite 3.30.1 [/wdata/tkoomar/VCFdbR/toy-data/table-genos/table-gt.db]
## tbls: ALT, contig, fileDate, fileformat, FILTER, FORMAT, gene_map, INFO,
## reference, samples, source, variant_geno, variant_impact, variant_info, VEP
The variant_info table contains the core information describing the
each variant (chromosome, position, ref, alt, etc), as well as the INFO
fields that describe each variant uniquely. Critically, this means that
each variant should be in variant_info exactly once.
By default, the variant_id, an, af and ac columns will be the
only one indexed (if present at all). To index additional columns for
quicker filtering (like the chr and start columns) see the
[Manually Indexing Columns] tutorial.
tbl(con, 'variant_info') %>% head()## # 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 1 22 1.61e7 1.61e7 A G 100 PASS c(NA… c(NA… <NA>
## 2 2 22 1.61e7 1.61e7 G A 100 PASS c(NA… c(NA… <NA>
## 3 3 22 1.61e7 1.61e7 C T 100 PASS c(NA… c(NA… <NA>
## 4 4 22 1.61e7 1.61e7 C T 100 PASS c(NA… c(NA… <NA>
## 5 5 22 1.61e7 1.61e7 C A 100 PASS c(NA… c(NA… <NA>
## 6 6 22 1.61e7 1.61e7 C A 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 <chr>, af <chr>, ns <int>, an <int>, eas_af <chr>, eur_af <chr>,
## # afr_af <chr>, amr_af <chr>, sas_af <chr>, dp <int>, aa <chr>, vt <chr>,
## # ex_target <int>, multi_allelic <int>
Notice that in this toy example, there are a number of columns that are
all NA (e.g. ciend). This occurs because these fields are declared
in the header, but do not actually appear in the VCF. Such ‘ghost’
variables would add significant bulk to the database. Removing specific
INFO
columns that
you don’t plan to use from the VCF before processing is always a good
idea.
Your database will only have a variant_impact table if it is annotated
with VEP using the default nomenclature (an INFO field named CSQ).
For specifics on what the various annotations mean, refer to the VEP
documentation.
As can be seen below, each variant may appear in the variant_impact
table multiple times, as it may effect multiple genes or transcripts.
This is why this impact information is on a separate table from the
variant_info table: everything on variant_info would be needlessly
duplicated multiple times if it were merged with the variant_impact
table.
tbl(con, 'variant_impact') %>% filter(variant_id %in% 840:841) %>% head()## # 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 840 upstream_g… MODIF… LA16c… ENSG… Transcript ENST00… proces… ""
## 2 840 intron_var… MODIF… NBEAP3 ENSG… Transcript ENST00… unproc… ""
## 3 840 non_coding… MODIF… NBEAP3 ENSG… Transcript ENST00… unproc… ""
## 4 841 upstream_g… MODIF… LA16c… ENSG… Transcript ENST00… proces… ""
## 5 841 intron_var… MODIF… NBEAP3 ENSG… Transcript ENST00… unproc… ""
## 6 841 non_coding… MODIF… NBEAP3 ENSG… Transcript ENST00… unproc… ""
## # … with 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>
By default, the variant_id, symbol, gene, consequence,
biotype, is_lof, is_exonic and is_splicing columns will be
indexed. To index additional columns for quicker filtering, see the
[Manually Indexing Columns] tutorial.
The variant_geno table describes the genotype of each sample for each
variant:
tbl(con, 'variant_geno') %>% head()## # 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 1 HG00096 0 0|0
## 2 1 HG00097 0 0|0
## 3 1 HG00099 0 0|0
## 4 1 HG00100 0 0|0
## 5 1 HG00101 0 0|0
## 6 1 HG00102 0 0|0
Notably, this data is in “long” format, compared to the VCF standard’s “wide” format. This is critical, as SQLite has a limit on the number of column in tables that would be violated quickly by large cohorts. Additionally, it allows variants to be quickly filtered based on genotype or other metrics (like genotype quality or read depth).
By default, the variant_id column will be indexed. To index additional
columns for quicker filtering, see the [Manually Indexing Columns]
tutorial.
You can see that there are quite a few other tables in the database. Most of these are summaries of various parts of the VCF header, and will look different depending on what is declared in the header of your VCF.
For example, INFO lists all of the various INFO columns from the VCF
and gives their description.
tbl(con, 'INFO') %>% collect() %>% slice(14:19)## # A tibble: 6 x 4
## name Number Type Description
## <chr> <chr> <chr> <chr>
## 1 AC A Integer Total number of alternate alleles in called genotypes
## 2 AF A Float Estimated allele frequency in the range (0,1)
## 3 NS 1 Integer Number of samples with data
## 4 AN 1 Integer Total number of alleles in called genotypes
## 5 EAS_AF A Float Allele frequency in the EAS populations calculated from…
## 6 EUR_AF A Float Allele frequency in the EUR populations calculated from…
None of these tables are indexed, as they are typically quite small.
If your VCF was annotated with VEP, a table called gene_map should
be placed in the database. It provides a summary of the various features
that actually appear in your variant_impact table. It can be very
useful, for example for iterating over all genes in the database.
tbl(con, 'gene_map') %>% head(7)## # Source: lazy query [?? x 10]
## # Database: sqlite 3.30.1
## # [/wdata/tkoomar/VCFdbR/toy-data/table-genos/table-gt.db]
## symbol symbol_source gene source feature canonical ensp ccds motif_name
## <chr> <chr> <chr> <chr> <chr> <chr> <chr> <chr> <chr>
## 1 <NA> <NA> <NA> <NA> <NA> <NA> <NA> <NA> <NA>
## 2 LA16c… Clone_based_… ENSG… Ensem… ENST00… YES <NA> <NA> <NA>
## 3 <NA> <NA> <NA> <NA> ENSR00… <NA> <NA> <NA> <NA>
## 4 LA16c… Clone_based_… ENSG… Ensem… ENST00… YES <NA> <NA> <NA>
## 5 LA16c… Clone_based_… ENSG… Ensem… ENST00… YES <NA> <NA> <NA>
## 6 <NA> <NA> <NA> <NA> ENSR00… <NA> <NA> <NA> <NA>
## 7 NBEAP3 HGNC ENSG… Ensem… ENST00… YES <NA> <NA> <NA>
## # … with 1 more variable: feature_type <chr>
You may have noticed that a couple .rds files are produced along with
the database itself. You can read them into R with the read_rds() or
read.rds() functions.
The [prefix]-chunk-ranges.rds file specifies GeomicRanges that each
contain ~1,000 variants in your input VCF. This is used when building
the VCFdb, but has little utility afterward.
read_rds('toy-data/table-genos/table-gt-chunk-ranges.rds') ## GRanges object with 2 ranges and 0 metadata columns:
## seqnames ranges strand
## <Rle> <IRanges> <Rle>
## 22.1 22 16050075-16064513 *
## 22.2 22 16064763-16139658 *
## -------
## seqinfo: 86 sequences from an unspecified genome; no seqlengths
More useful is the [prefix]-variant-ranges.rds file. It contains 1
GRange per variant. These can be used to quickly subset your variants by
position, or to generate new annotations you might want to add to the
variant_info table after the database is created. If you are not
familiar with GenomicRanges, they are the workhorse of the
Bioconductor genomics tools, an extremely powerful tool, and worth
learning more
about
read_rds('toy-data/table-genos/table-gt-variant-ranges.rds') ## GRanges object with 1000 ranges and 3 metadata columns:
## seqnames ranges strand | variant_id ref alt
## <Rle> <IRanges> <Rle> | <integer> <character> <character>
## [1] 22 16050075 * | 1 A G
## [2] 22 16050115 * | 2 G A
## [3] 22 16050213 * | 3 C T
## [4] 22 16050319 * | 4 C T
## [5] 22 16050527 * | 5 C A
## ... ... ... ... . ... ... ...
## [996] 22 16139438 * | 996 T C
## [997] 22 16139473 * | 997 C A
## [998] 22 16139477 * | 998 A G
## [999] 22 16139539 * | 999 T C
## [1000] 22 16139658 * | 1000 C A
## -------
## seqinfo: 1 sequence from an unspecified genome; no seqlengths
Note Many Bioconductor packages (particularly GenomicRanges and
VariantAnnotation) have lots of namespace conflicts with
dplyr/dbplyr. Most notably, this includes the select() function,
but also functions like reduce(), expand() and slice(). Be careful
about the order you load packages, or prepare to append package names to
lots of your function calls (e.g. dplyr::select()).
The least tedious option tends to be to source VariantAnnotation
first (it loads many other useful Bioconductor packages), and then
dbplyr/tidyverse:
library(VariantAnnotation)
library(tidyverse)
library(dbplyr)DBI::dbDisconnect(con)