Mapping QTLs in outbred mice using varbvs

The hardware and bandwidth for this mirror is donated by dogado GmbH, the Webhosting and Full Service-Cloud Provider. Check out our Wordpress Tutorial.
If you wish to report a bug, or if you are interested in having us mirror your free-software or open-source project, please feel free to contact us at mirror[@]dogado.de.

2023-05-31

In this vignette, we use varbvs to map QTLs for phenotypes measured in CFW (Carworth Farms White) outbred mice. Phenotypes include muscle weights—EDL and soleus muscle—and testis weight measured at sacrifice. Running this script with trait = "testis" reproduces the results and figures shown in the second example of a forthcoming paper (Carbonetto et al, 2016).

Vignette parameters

Begin by loading packages into the R environment.

library(curl)
library(lattice)
library(varbvs)

These script parameters specify the candidate prior log-odds settings, the prior variance of the coefficients, and which trait to analyze. Set trait to “edl”, “soleus” or “testis”.

trait      <- "testis"
covariates <- "sacwt"
logodds    <- seq(-5,-3,0.25)
sa         <- 0.05

Set the random number generator seed.

set.seed(1)

Load the genotype and phenotype data

Retrieve the data from the Zenodo repository.

load(curl("https://zenodo.org/record/546142/files/cfw.RData"))

Only analyze samples for which the phenotype and all the covariates are observed.

rows <- which(apply(pheno[,c(trait,covariates)],1,
                    function (x) sum(is.na(x)) == 0))
pheno <- pheno[rows,]
geno  <- geno[rows,]

Fit variational approximation to posterior

runtime <- system.time(fit <-
  varbvs(geno,as.matrix(pheno[,covariates]),pheno[,trait],
         sa = sa,logodds = logodds,verbose = FALSE))
cat(sprintf("Model fitting took %0.2f minutes.\n",runtime["elapsed"]/60))

Summarize the results of model fitting

print(summary(fit))

Show three genome-wide scans: (1) one using the posterior inclusion probabilities (PIPs) computed in the BVS analysis of all SNPs; (2) one using the p-values computed using GEMMA; and (3) one using the PIPs computed from the BVSR model in GEMMA.

trellis.par.set(axis.text     = list(cex = 0.7),
                par.ylab.text = list(cex = 0.7),
                par.main.text = list(cex = 0.7,font = 1))
j <- which(fit$pip > 0.5)
r <- gwscan.gemma[[trait]]
r[is.na(r)] <- 0
chromosomes   <- levels(gwscan.bvsr$chr)
xticks        <- rep(0,length(chromosomes))
names(xticks) <- chromosomes
pos           <- 0
for (i in chromosomes) {
  n         <- sum(gwscan.bvsr$chr == i)
  xticks[i] <- pos + n/2
  pos       <- pos + n + 25
}
print(plot(fit,groups = map$chr,vars = j,gap = 1500,cex = 0.6,
           ylab = "probability",main = "a. multi-marker (varbvs)",
           scales = list(y = list(limits = c(-0.1,1.2),at = c(0,0.5,1))),
           vars.xyplot.args = list(cex = 0.6)),
      split = c(1,1,1,3),more = TRUE)
print(plot(fit,groups = map$chr,score = r,vars = j,cex = 0.6,gap = 1500,
           draw.threshold = 5.71,ylab = "-log10 p-value",
           main = "b. single-marker (GEMMA -lm 2)",
           scales = list(y = list(limits = c(-1,20),at = seq(0,20,5))),
           vars.xyplot.args = list(cex = 0.6)),
     split = c(1,2,1,3),more = TRUE)
print(xyplot(p1 ~ plot.x,gwscan.bvsr,pch = 20,col = "midnightblue",
             scales = list(x = list(at = xticks,labels = chromosomes),
                           y = list(limits = c(-0.1,1.2),at = c(0,0.5,1))),
             xlab = "",ylab = "probability",main = "c. multi-marker (BVSR)"),
      split = c(1,3,1,3),more = FALSE)

Session information

This is the version of R and the packages that were used to generate these results.

sessionInfo()

These binaries (installable software) and packages are in development.
They may not be fully stable and should be used with caution. We make no claims about them.
Health stats visible at Monitor.