In kullrich/vGWAS: Variance Heterogeneity Genome-wide Association Study - Reimplementation

knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>"
)

Vignette - vGWAS basic tutorial

A Brief Tutorial of the R Package vGWAS

Setup Notes

To use the vGWAS package, of course, an R environment is required. Visit:

http://www.r-project.org

and install R for the operating system.

Start R and in the R console, type the following command to install the package dependencies first:

#install.packages("devtools")
#install.packages("knitr")
#install.packages("dglm")
#install.packages("doParallel")
#install.packages("foreach")
#install.packages("genio")
#install.packages("hglm")
#install.packages("Matrix")
#install.packages("onewaytests")
#library(devtools)
#devtools::install_github("kullrich/vGWAS", build_vignettes = TRUE, dependencies = FALSE)

Now the package is installed in the R library.

Example

Load the library in the R console, and open this vignette with the command:

library(vGWAS)

Nine main functions in the package are ready to use:

1. bfmedian.test
2. brown.forsythe.tes
3. getMAF
4. vGWAS
5. vGWASparallel
6. vGWAS.gc
7. vGWAS.variance
8. plot (S3 method for vGWAS object)
9. summary (S3 method for vGWAS object)

Run the following commands to load the example data:

data(pheno)
data(geno)
data(chr)
data(map)

pheno is a numeric vector of the simulated phenotypic values. By running:

hist(pheno, breaks = 30, density = 15, col = "darkred")

The command:

table(chr)

shows exactly the number of markers on each of the five simulated chromosomes.

To pre-filter the geno.matrix data for minor-allele-frequency, use the getMAF function as follows:

data(geno)
dim(geno)
maf <- getMAF(geno.matrix = geno, geno.snp = "col")
geno.maf <- geno[, maf > 0.05]
dim(geno.maf)
chr.maf <- chr[maf > 0.05]
map.maf <- map[maf > 0.05]

Now, the objects loaded in R are ready for a vGWA scan, which can be done using the single command:

vgwa <- vGWAS(
  phenotype = pheno,
  geno.matrix = geno,
  marker.map = map,
  chr.index = chr,
  pB = FALSE)

A progress bar will indicate the progress of the scan if pB is set to TRUE.

The same can be run using multiple cores.

Note: using the vGWASparallel function one need to specify, if the SNPs in the geno.matrix are at columns (geno.snp = "col") or at rows (geno.snp = "row")

#vgwa <- vGWASparallel(
#  phenotype = pheno,
#  geno.matrix = geno,
#  geno.snp = "col",
#  marker.map = map,
#  chr.index = chr,
#  pB = FALSE,
#  ncores = 2)

When the scan is finished, all the output statistics will be returned as a list into the object vgwa, which belongs to the class ‘vGWAS’. Any object that has a structure belonging to class ‘vGWAS’ can be directly passed into S3 method function plot. For instance, simply run the following command, we can plot the results in vgwa:

plot(vgwa)

which produces Figure 2. There is a clear peak above the Bonferroni corrected threshold (dashed orange line).

Regarding the marker that gave the highest score, the heritability explained by the mean and variance can be split and calculated via:

# get marker with lowest p.value
marker.lowest <- vgwa[["p.value"]] == min(vgwa[["p.value"]])
vGWAS.variance(
  phenotype = pheno,
  marker.genotype = geno[, marker.lowest])

The output can also be stored if assigning the function call to an object.

Remarks

The package source and further development information are on the R-Forge project page: https://r-forge.r-project.org/projects/vgwas/

A reimplementation of vGWAS introducing the parallel function vGWASparallel as a devtool package can be found here: https://github.com/kullrich/vGWAS

kullrich/vGWAS documentation built on June 10, 2025, 3:56 a.m.