README.Rmd

---
output: github_document
---

<!-- README.md is generated from README.Rmd. Please edit README.Rmd -->



```{r setup, include = FALSE}
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>",
  fig.path = "man/figures/README-",
  out.width = "100%",
  message = FALSE,
  warning = FALSE
)
```

# qtl2pleio

[![R-CMD-check](https://github.com/fboehm/qtl2pleio/workflows/R-CMD-check/badge.svg)](https://github.com/fboehm/qtl2pleio/actions)
[![CRAN_Status_Badge](https://www.r-pkg.org/badges/version/qtl2pleio)](https://cran.r-project.org/package=qtl2pleio)
[![Coverage Status](https://img.shields.io/codecov/c/github/fboehm/qtl2pleio/master.svg)](https://codecov.io/github/fboehm/qtl2pleio?branch=master)
[![Project Status: Active – The project has reached a stable, usable state and is being actively developed.](https://www.repostatus.org/badges/latest/active.svg)](https://www.repostatus.org/#active)
[![status](https://joss.theoj.org/papers/66bca5dc3d2e72b6259159bad07aafaf/status.svg)](https://joss.theoj.org/papers/10.21105/joss.01435)
[![DOI](https://zenodo.org/badge/104493705.svg)](https://zenodo.org/badge/latestdoi/104493705)




## Overview

`qtl2pleio` is a software package for use with the [R statistical computing environment](https://cran.r-project.org/). `qtl2pleio` is freely available for download and use. I share it under the [MIT license](https://opensource.org/licenses/mit-license.php). The user will also want to download and install the [`qtl2` R package](https://kbroman.org/qtl2/).

Click [here](https://mybinder.org/v2/gh/fboehm/qtl2pleio/master?urlpath=rstudio) to explore `qtl2pleio` within a live [Rstudio](https://rstudio.com/) session in "the cloud".

## Contributor guidelines

We eagerly welcome contributions to `qtl2pleio`. All pull requests will be considered. Features requests and bug reports may be filed as Github issues. All contributors must abide by the [code of conduct](https://github.com/fboehm/qtl2pleio/blob/master/CONDUCT.md). 

## Technical support 

For technical support, please open a Github issue. If you're just getting started with `qtl2pleio`, please examine the [vignettes](https://fboehm.us/software/qtl2pleio) on the [package's web site](https://fboehm.us/software/qtl2pleio). You can also email <frederick.boehm@gmail.com> for assistance. 




## Goals

The goal of `qtl2pleio` is, for a pair of traits that show evidence for
a QTL in a common region, to distinguish between pleiotropy (the null
hypothesis, that they are affected by a common QTL) and the
alternative that they are affected by separate QTL. It extends the
likelihood ratio test of [Jiang and Zeng
(1995)](https://www.genetics.org/content/140/3/1111.long) for
multiparental populations, such as Diversity Outbred mice, including
the use of multivariate polygenic random effects to account for population structure.
`qtl2pleio` data structures are those used in the
[`rqtl/qtl2`](https://kbroman.org/qtl2/) package.


## Installation


To install qtl2pleio, use `install_github()` from the
[devtools](https://devtools.r-lib.org) package.


```{r install-qtl2pleio, eval = FALSE}
install.packages("qtl2pleio")
```

You may also wish to install the [R/qtl2](https://kbroman.org/qtl2/). We
will use it below.

```{r install-qtl2, eval = FALSE}
install.packages("qtl2")
```


## Example

Below, we walk through an example analysis with Diversity Outbred
mouse data. We need a number of preliminary steps before we can
perform our test of pleiotropy vs. separate QTL. Many procedures rely
on the R package `qtl2`. We first load the `qtl2`
and `qtl2pleio` packages.

```{r pkgs}
library(qtl2)
library(qtl2pleio)
library(ggplot2)
```


### Reading data from `qtl2data` repository on github

We'll consider the
[`DOex`](https://github.com/rqtl/qtl2data/tree/master/DOex/) data in
the [`qtl2data`](https://github.com/rqtl/qtl2data/) repository.
We'll download the DOex.zip file before calculating founder allele dosages.

```{r download-DOex}
file <- paste0("https://raw.githubusercontent.com/rqtl/",
               "qtl2data/master/DOex/DOex.zip")
DOex <- read_cross2(file)
```


```{r calc-genoprobs}
probs <- calc_genoprob(DOex)
```


Let's calculate the founder allele dosages from the 36-state genotype probabilities.

```{r calc-allele-probs}
pr <- genoprob_to_alleleprob(probs)
```


We now have an allele probabilities object stored in `pr`.

```{r check-pr}
names(pr)
dim(pr$`2`)
```

We see that `pr` is a list of 3 three-dimensional arrays - one array for each of 3 chromosomes.

### Kinship calculations

For our statistical model, we need a kinship matrix. We get one with the `calc_kinship` function in the `rqtl/qtl2` package.


```{r calc-kinship}
kinship <- calc_kinship(probs = pr, type = "loco")
```

### Statistical model specification

We use the multivariate linear mixed effects model:

$$
\text{vec}(Y) = X \text{vec}(B) + \text{vec}(G) + \text{vec}(E)
$$

where $Y$ contains phenotypes, X contains founder allele probabilities and covariates, and B contains founder allele effects. G is the polygenic random effects, while E is the random errors. We provide more details in the vignette.



### Simulating phenotypes with `qtl2pleio::sim1`


The function to simulate phenotypes in `qtl2pleio` is `sim1`.

```{r pp-def}
# set up the design matrix, X
pp <- pr[[2]] #we'll work with Chr 3's genotype data
```


```{r X-def}
#Next, we prepare a design matrix X
X <- gemma2::stagger_mats(pp[ , , 50], pp[ , , 50])
```

```{r B-def}
# assemble B matrix of allele effects
B <- matrix(data = c(-1, -1, -1, -1, 1, 1, 1, 1, -1, -1, -1, -1, 1, 1, 1, 1), nrow = 8, ncol = 2, byrow = FALSE)
# set.seed to ensure reproducibility
set.seed(2018-01-30)
Sig <- calc_Sigma(Vg = diag(2), Ve = diag(2), kinship = kinship[[2]])
# call to sim1
Ypre <- sim1(X = X, B = B, Sigma = Sig)
Y <- matrix(Ypre, nrow = 261, ncol = 2, byrow = FALSE)
rownames(Y) <- rownames(pp)
colnames(Y) <- c("tr1", "tr2")
```

Let's perform univariate QTL mapping for each of the two traits in the Y matrix.

```{r 1d-scans}
s1 <- scan1(genoprobs = pr, pheno = Y, kinship = kinship)
```

Here is a plot of the results.

```{r 1d-lod-plots, include=FALSE}
plot(s1, DOex$pmap)
plot(s1, DOex$pmap, lod=2, col="violetred", add=TRUE)
legend("topleft", colnames(s1), lwd=2, col=c("darkslateblue", "violetred"), bg="gray92")
```

```{r, echo = FALSE}
knitr::include_graphics("https://raw.githubusercontent.com/fboehm/qtl2pleio/master/man/figures/README-1d-lod-plots-1.png")
```

And here are the observed QTL peaks with LOD > 8.

```{r find-peaks}
find_peaks(s1, map = DOex$pmap, threshold=8)
```





### Perform two-dimensional scan as first step in pleiotropy vs. separate QTL hypothesis test


We now have the inputs that we need to do a pleiotropy vs. separate QTL test. We have the founder allele dosages for one chromosome, *i.e.*, Chr 3, in the R object `pp`, the matrix of two trait measurements in `Y`, and a LOCO-derived kinship matrix, `kinship[[2]]`.


```{r 2d-scan}
out <- suppressMessages(scan_pvl(probs = pp,
                pheno = Y,
                kinship = kinship[[2]], # 2nd entry in kinship list is Chr 3
                start_snp = 38,
                n_snp = 25
                ))
```


### Create a profile LOD plot to visualize results of two-dimensional scan

To visualize results from our two-dimensional scan, we calculate profile LOD for each trait. The code below makes use of the R package `ggplot2` to plot profile LODs over the scan region.


```{r profile-plot, include = TRUE}
library(dplyr)
out %>%
  calc_profile_lods() %>%
  add_pmap(pmap = DOex$pmap$`3`) %>%
  ggplot() + geom_line(aes(x = marker_position, y = profile_lod, colour = trait))
```




### Calculate the likelihood ratio test statistic for pleiotropy v separate QTL

We use the function `calc_lrt_tib` to calculate the likelihood ratio test statistic value for the specified traits and specified genomic region.

```{r lrt-calc}
(lrt <- calc_lrt_tib(out))
```

### Bootstrap analysis to get p-values


Before we call `boot_pvl`, we need to identify the index (on the chromosome under study) of the marker that maximizes the likelihood under the pleiotropy constraint. To do this, we use the `qtl2pleio` function `find_pleio_peak_tib`.

```{r get-pleio-index}
(pleio_index <- find_pleio_peak_tib(out, start_snp = 38))
```



```{r boot}
set.seed(2018-11-25) # set for reproducibility purposes.
b_out <- suppressMessages(boot_pvl(probs = pp,
         pheno = Y,
         pleio_peak_index = pleio_index,
         kinship = kinship[[2]], # 2nd element in kinship list is Chr 3
         nboot = 10,
         start_snp = 38,
         n_snp = 25
         ))
```


```{r pval}
(pvalue <- mean(b_out >= lrt))
```



## Citation information

```{r cite}
citation("qtl2pleio")
```