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

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

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.

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

```
<- system.time(fit <-
runtime 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))
<- which(fit$pip > 0.5)
j <- gwscan.gemma[[trait]]
r is.na(r)] <- 0
r[<- levels(gwscan.bvsr$chr)
chromosomes <- rep(0,length(chromosomes))
xticks names(xticks) <- chromosomes
<- 0
pos for (i in chromosomes) {
<- sum(gwscan.bvsr$chr == i)
n <- pos + n/2
xticks[i] <- pos + n + 25
pos
}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()`