Multivariate, one-QTL scans

Source: vignettes/multivariate-one-QTL-scans.Rmd

Our goal is to perform a multivariate one-QTL scan across the genome. We’ll also do an approximate permutation test to get p-values for any QTL that we find.

The analysis is multivariate in that we’re analyzing three traits simultaneously. The “one-QTL” phrase denotes the fact that we’re looking for a single QTL that affects the three traits simultaneously.

Load Recla data from qtl2data github repository

We illustrate functions from qtl2pleio by analyzing data from 261 Diversity Outbred mice [@recla2014precise,@logan2013high].

file <- paste0("https://raw.githubusercontent.com/rqtl/",
               "qtl2data/master/DO_Recla/recla.zip")
recla <- read_cross2(file)
# make sex a covariate for use in scans
recla[[6]][ , 1, drop = FALSE] -> sex
# insert pseudomarkers
insert_pseudomarkers(recla, step = 0.10) -> pseudomap

We use the hidden Markov model from @broman2012genotype and @broman2012haplotype (as implemented in @broman2019rqtl2) to calculate 36-state genotype probabilities for autosomal markers.

probs <- calc_genoprob(recla, map = pseudomap, cores = 1)

We then convert the genotype probabilities to haplotype dosages.

aprobs <- genoprob_to_alleleprob(probs)

We calculate kinship matrices, by the “leave one chromosome out (loco)” method.

kinship <- calc_kinship(aprobs, "loco")

Before performing our QTL mapping, we transform the phenotypes.

recla$pheno -> ph
log(ph) -> lph
apply(FUN = broman::winsorize, X = lph, MARGIN = 2) -> wlph
tibble::as_tibble(wlph) -> wlph_tib

We next perform the bivariate, one-QTL scan for two correlated traits.

sex2 <- matrix(as.numeric(sex == "female"), ncol = 1)
colnames(sex2) <- "female"
rownames(sex2) <- rownames(aprobs$`1`)
system.time(out <- qtl2pleio::scan_multi_onechr(probs = aprobs$`1`, 
             pheno = wlph[, c(7, 10)], 
             kinship = NULL, 
             addcovar = sex2
             )
)
#> Warning: replacing previous import 'vctrs::data_frame' by 'tibble::data_frame'
#> when loading 'dplyr'
#>    user  system elapsed 
#> 634.608  74.049 667.948
out2 <- qtl2pleio::scan_multi_oneqtl(probs_list = list(aprobs$`1`), pheno = wlph[, c(7, 10)], kinship_list = NULL, addcovar = sex2)