BeviMed with VCFs

Daniel Greene

2021-01-31

The BeviMed package comes with a script containing functions to simplify reading allele count matrices from VCF files. The functions depend on tabix, but have the advantage of allowing variants in local regions to be read in, reducing the amount of memory consumed at any one time. However, if you want to analyse many regions, it may be more efficient to read in larger parts of the file - in which case, a package such as vcfR might be more appropriate.

In order to make the functions available, we must source the script:

library(BeviMed)
source(paste0(system.file(package="BeviMed", "/scripts/vcf.R")))

The script creates the function vcf2matrix, which depends on the external program tabix (available from http://www.htslib.org/download/) for reading allele count matrices from VCF files. It uses arguments:

You can invoke the function simply to obtain the allele count matrix and pass straight to bevimed, along with phenotype label:

ac_matrix <- vcf2matrix("my-vcf.vcf.gz", chr="2", from=1, to=1e4)
pheno <- read.table(file="my-phenotype-data.txt", header=TRUE)

bevimed(y=pheno$disease_status, G=ac_matrix)