Repository navigation
Schema of File 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/file-genos/file-gt.db')src_dbi(con)## src: sqlite 3.30.1 [/wdata/tkoomar/VCFdbR/toy-data/file-genos/file-gt.db]
## tbls: ALT, contig, fileDate, fileformat, FILTER, FORMAT, gene_map, INFO,
## reference, samples, source, 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 35]
## # Database: sqlite 3.30.1
## # [/wdata/tkoomar/VCFdbR/toy-data/file-genos/file-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 24 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>, geno <chr>
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.
When a database is built in File-GT mode, the variant_info table will
contain a geno column that has the absolute path to the genotype file
corresponding to that varint. This makes it easy to lapply/map over
this column to read in genotypes (see [Working with File GT Format
Databases]).
Additionally, it means that you can move the database, but leave the
folder containing the genotypes in place (e.g. place the database on an
SSD and leave the genotypes on a larger but slower hard disk). If you
move the [prefix]-geno directory, the geno column will no longer
point to the correct location. You will need to rewite this column, or
remove it entirely and use paste0() or similar to genrate the path to
the geno files.
tbl(con, 'variant_info') %>% select(variant_id, geno) %>% head()## # Source: lazy query [?? x 2]
## # Database: sqlite 3.30.1
## # [/wdata/tkoomar/VCFdbR/toy-data/file-genos/file-gt.db]
## variant_id geno
## <int> <chr>
## 1 1 /wdata/tkoomar/VCFdbR/toy-data/file-genos/file-gt-genos/1.rds
## 2 2 /wdata/tkoomar/VCFdbR/toy-data/file-genos/file-gt-genos/2.rds
## 3 3 /wdata/tkoomar/VCFdbR/toy-data/file-genos/file-gt-genos/3.rds
## 4 4 /wdata/tkoomar/VCFdbR/toy-data/file-genos/file-gt-genos/4.rds
## 5 5 /wdata/tkoomar/VCFdbR/toy-data/file-genos/file-gt-genos/5.rds
## 6 6 /wdata/tkoomar/VCFdbR/toy-data/file-genos/file-gt-genos/6.rds
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/file-genos/file-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 genotypes of a File GT database are stored in a special,
write-protected directory. They are .rds files named after the variant
to which they correspond. Because 1 file is generated per variant, this
can result in a directory with millions of files. Make sure your
filesystem is capable of handling this scenario before you create a File
GT format
database.
list.files("toy-data/file-genos/file-gt-genos") %>% head()## [1] "1000.rds" "100.rds" "101.rds" "102.rds" "103.rds" "104.rds"
These files can be read in with read_rds()/read.rds():
read_rds("toy-data/file-genos/file-gt-genos/1.rds") %>% head()## variant_id sample gt gt_raw
## 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).
The more samples present in your cohort, the larger these files will
become. Removing FORMAT
fields from
your VCF that you do not plan to use is one of the most effective ways
to save disk space and reduce database creation time.
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/file-genos/file-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/file-genos/file-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/file-genos/file-gt-variant-ranges.rds') ## GRanges object with 1000 ranges and 4 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
## geno
## <character>
## [1] /wdata/tkoomar/VCFdbR/toy-data/file-genos/file-gt-genos/1.rds
## [2] /wdata/tkoomar/VCFdbR/toy-data/file-genos/file-gt-genos/2.rds
## [3] /wdata/tkoomar/VCFdbR/toy-data/file-genos/file-gt-genos/3.rds
## [4] /wdata/tkoomar/VCFdbR/toy-data/file-genos/file-gt-genos/4.rds
## [5] /wdata/tkoomar/VCFdbR/toy-data/file-genos/file-gt-genos/5.rds
## ... ...
## [996] /wdata/tkoomar/VCFdbR/toy-data/file-genos/file-gt-genos/996.rds
## [997] /wdata/tkoomar/VCFdbR/toy-data/file-genos/file-gt-genos/997.rds
## [998] /wdata/tkoomar/VCFdbR/toy-data/file-genos/file-gt-genos/998.rds
## [999] /wdata/tkoomar/VCFdbR/toy-data/file-genos/file-gt-genos/999.rds
## [1000] /wdata/tkoomar/VCFdbR/toy-data/file-genos/file-gt-genos/1000.rds
## -------
## 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)