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).

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)`

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,]
```

```
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))
```

`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)
```

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

`sessionInfo()`