Skip to content
Testing pleiotropy vs. separate QTL in multiparental populations
Branch: master
Clone or download
Fetching latest commit…
Cannot retrieve the latest commit at this time.
Permalink
Type Name Latest commit message Commit time
Failed to load latest commit information.
R
inst
man
src
tests
vignettes
.Rbuildignore
.gitignore
.travis.yml
CONDUCT.md
DESCRIPTION
LICENSE
NAMESPACE
NEWS.md
README.Rmd
README.md
_pkgdown.yml
codecov.yml
codemeta.json
install.R
qtl2pleio.Rproj
runtime.txt

README.md

qtl2pleio

Binder Travis-CI Build Status Coverage Status Project Status: WIP – Initial development is in progress, but there has not yet been a stable, usable release suitable for the public.

Overview

qtl2pleio is a software package for use with the R statistical computing environment. qtl2pleio is freely available for download and use. I share it under the MIT license. The user will also want to download and install the qtl2 R package.

Click here to explore qtl2pleio within a live Rstudio session in “the cloud”.

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

Installation

To install qtl2pleio, use install_github() from the devtools package.

install.packages("devtools")
devtools::install_github("fboehm/qtl2pleio")

You may also wish to install R/qtl2 and the qtl2convert package. We will use both below.

install.packages(c("qtl2", "qtl2convert"), repos="http://rqtl.org/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, qtl2convert, and qtl2pleio packages.

library(qtl2)
library(qtl2convert)
library(qtl2pleio)

Reading data from qtl2data repository on github

We’ll consider the DOex data in the qtl2data repository. We’ll download the DOex.zip file before calculating founder allele dosages.

file <- paste0("https://raw.githubusercontent.com/rqtl/",
               "qtl2data/master/DOex/DOex.zip")
DOex <- read_cross2(file)
probs <- calc_genoprob(DOex)

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

pr <- genoprob_to_alleleprob(probs)

We now have an allele probabilities object stored in pr.

names(pr)
#> [1] "2" "3" "X"
dim(pr$`2`)
#> [1] 261   8 127

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.

kinship <- calc_kinship(probs = pr, type = "loco")
str(kinship)
#> List of 3
#>  $ 2: num [1:261, 1:261] 0.6934 0.0705 0.2356 0.0558 0.0513 ...
#>   ..- attr(*, "n_pos")= int 195
#>   ..- attr(*, "dimnames")=List of 2
#>   .. ..$ : chr [1:261] "1" "4" "5" "6" ...
#>   .. ..$ : chr [1:261] "1" "4" "5" "6" ...
#>  $ 3: num [1:261, 1:261] 0.6662 0.0647 0.2024 0.1129 0.0772 ...
#>   ..- attr(*, "n_pos")= int 220
#>   ..- attr(*, "dimnames")=List of 2
#>   .. ..$ : chr [1:261] "1" "4" "5" "6" ...
#>   .. ..$ : chr [1:261] "1" "4" "5" "6" ...
#>  $ X: num [1:261, 1:261] 0.4871 0.0831 0.1953 0.1043 0.1125 ...
#>   ..- attr(*, "n_pos")= int 229
#>   ..- attr(*, "dimnames")=List of 2
#>   .. ..$ : chr [1:261] "1" "4" "5" "6" ...
#>   .. ..$ : chr [1:261] "1" "4" "5" "6" ...

Statistical model specification

We use the multivariate linear mixed effects model:

[vec(Y) = X vec(B) + vec(G) + 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.

# set up the design matrix, X
pp <- pr[[2]] #we'll work with Chr 3's genotype data
#Next, we prepare a design matrix X
X <- gemma2::stagger_mats(pp[ , , 50], pp[ , , 50])
# 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)
# call to sim1
Ypre <- sim1(X = X, B = B, Vg = diag(2), Ve = diag(2), kinship = kinship[[2]])
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.

s1 <- scan1(genoprobs = pr, pheno = Y, kinship = kinship)

Here is a plot of the results.

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

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

find_peaks(s1, map = DOex$pmap, threshold=8)
#>   lodindex lodcolumn chr       pos      lod
#> 1        1       tr1   3  82.77806 16.19704
#> 2        2       tr2   3  82.77806 18.26406
#> 3        2       tr2   X 103.79061 16.19708

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

out <- scan_pvl(probs = pp,
                pheno = Y,
                kinship = kinship[[2]], # 2nd entry in kinship list is Chr 3
                start_snp = 38,
                n_snp = 25, n_cores = 1
                )

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.

library(dplyr)
out %>%
  tidy_scan_pvl(DOex$pmap[[2]]) %>% # pmap[[2]] is physical map for Chr 3
  add_intercepts(intercepts_univariate = c(82.8, 82.8)) %>%
  plot_pvl(phenames = c("tr1", "tr2"))

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.

(lrt <- calc_lrt_tib(out))
#> [1] 0

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.

(pleio_index <- find_pleio_peak_tib(out, start_snp = 38))
#> loglik13 
#>       50
set.seed(2018-11-25) # set for reproducibility purposes.
system.time(b_out <- boot_pvl(probs = pp,
         pheno = Y,
         pleio_peak_index = pleio_index,
         kinship = kinship[[2]], # 2nd element in kinship list is Chr 3
         nboot_per_job = 10,
         start_snp = 38,
         n_snp = 25
         ))
#>    user  system elapsed 
#>  77.357   1.996  79.458
(pvalue <- mean(b_out >= lrt))
#> [1] 1

Code of Conduct

Please note that this project is released with a Contributor Code of Conduct. By participating in this project you agree to abide by its terms.

Citation information

citation("qtl2pleio")
#> 
#> To cite qtl2pleio in publications use:
#> 
#>   Frederick Boehm (2019). qtl2pleio: Testing pleiotropy vs.
#>   separate QTL in multiparental populations. R package version
#>   0.1.2.9001. URL https://github.com/fboehm/qtl2pleio
#> 
#> A BibTeX entry for LaTeX users is
#> 
#>   @Manual{,
#>     title = {qtl2pleio: Testing pleiotropy vs. separate QTL in multiparental populations},
#>     author = {Frederick Boehm},
#>     year = {2019},
#>     note = {R package version 0.1.2.9001},
#>     url = {https://github.com/fboehm/qtl2pleio},
#>   }

Icon image of mice: https://www.jax.org/~/media/JaxWeb/images/jax-mice-and-services/mice/datasheets/002468

Icon image of Arabidopsis: https://sites.cns.utexas.edu/sites/default/files/styles/os_files_large/public/juenger_lab/files/pp1244cvr001-1.gif

Icon image of Drosophila: https://www.yourgenome.org/sites/default/files/styles/banner/public/banners/stories/fruit-flies-in-the-laboratory/single-fruit-fly-drosophila-melanogaster-on-white-background-cropped.jpg

You can’t perform that action at this time.