In the course of a genome-wide association study, the situation often arises that some phenotypes are known with greater precision than others. It could be that some individuals are known to harbor more micro-environmental variance than others. In the case of inbred strains of model organisms, it could be the case that more organisms were observed from some strains than others, so the strains with more organisms have better-estimated means.
wISAM handles this situation by allowing for weighting of each observation accoding to residual variance. Specifically, the
weight parameter to the function
conduct_scan takes the precision of each observation (one over the variance).
You can install wISAM from github with:
or from CRAN with:
Load the package and the example data will be available automatically. You may want to have a look at the data format and shape with
str(). Create a standard genome scan (
gs, below) or one that uses weights (
wgs, below). You can see that this package uses a reference class to implement the GenomeScan object.
library(wISAM) gs <- GenomeScan$new(y = phenotype, X = covariate_mat, G = locus_list, K = kinship_mat) wgs <- GenomeScan$new(y = phenotype, X = covariate_mat, G = locus_list, K = kinship_mat, w = 1/se_mean_per_strain)
You can then call
prep_scan, which computes the quantities that will only need to be computed once in the course of the genome scan. Depending on the size of the population, this can be quite slow. The methods of the GenomeScan reference class return the instance, so you can chain them together. Unfortunately, this means you will need to wrap the calls in invisible, or save them to an object to avoid getting a big print to screen.
invisible(gs$prep_scan()) #> Preparing GenomeScan... invisible(wgs$prep_scan()) #> Preparing GenomeScan...
You may want to investigate what changed in
wgs when you ran
prep_scan(). The next step is to conduct the genome scan.
invisible(gs$conduct_scan()) #> Conducting GenomeScan... invisible(wgs$conduct_scan()) #> Conducting GenomeScan...
After conducting a scan, you can have a look at the output.
plot(x = gs$Results$LR, col = c('red', rep('black', 99)))
plot(x = wgs$Results$LR, col = c('red', rep('black', 99)))
plot(x = gs$Results$LR, y = wgs$Results$LR, col = c('red', rep('black', 99))) abline(a = 0, b = 1)
You can see that the weighted and unweighted version give similar results, but
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.