.libPaths("~/work/libs")
knitr::opts_chunk$set(dev="CairoPNG")
knitr::opts_chunk$set(echo = TRUE)

Read the csv files

cluster_traits <- readr::read_csv("../data/pheno_eigenscores_from_WGCNA_all60_win.csv")
load("../data/reduced_map_qtl2_mapping_objects.Rdata")
kinship <- qtl2::calc_kinship(probs = genoprobs.1, type = "loco")
phe <- as.matrix(cluster_traits[, c("greenyellow", "magenta")])
rownames(phe) <- cluster_traits$X1
start <- max(which(map.1$`10` < 20)) 
end <- min(which(map.1$`10` > 45))
library(magrittr)
library(ggplot2)

WCGNA traits

library(future.batchtools)
library(future)
options(future.globals.maxSize= 4 * 1024 ^ 3)
plan("batch_multicore")
sp_out <- qtl2pleio::scan_pvl(probs = genoprobs.1$`10`, 
                                pheno = phe,
                                kinship = kinship$`10`,
                                start_snp = start,
                                n_snp = end - start + 1
                                )
(lrt <- qtl2pleio::calc_lrt_tib(sp_out))
sp_out %>%
  qtl2pleio::calc_profile_lods() %>%
  qtl2pleio::add_pmap(pmap = map.1$`10`) %>%
  ggplot2::ggplot() + ggplot2::geom_line(ggplot2::aes(x = marker_position, y = profile_lod, colour = trait))
bfn <- "../data/bivariate-greenyellow-magenta-Chr10-boot.rds"

if (!file.exists(bfn)){
  (pp_index <- qtl2pleio::find_pleio_peak_tib(tib = sp_out, 
                                             start_snp = start
                                             ))
  set.seed(3411192)
  b_out <- qtl2pleio::boot_pvl(probs = genoprobs.1$`10`, 
                      pheno = phe,
                      kinship = kinship$`10`,
                      start_snp = start,
                      n_snp = end - start + 1, 
                      pleio_peak_index = pp_index,
                      nboot = 400
                     )
  (pvalue <- mean(b_out >= lrt))
} else {b_out <- readRDS(bfn)}
(pvalue <- mean(b_out >= lrt))

Colophon

This report was generated on r Sys.time() using the following computational environment and dependencies:

# which R packages and versions?
if ("devtools" %in% installed.packages()) devtools::session_info()

The current Git commit details are:

# what commit is this file at? 
if ("git2r" %in% installed.packages() & git2r::in_repository(path = ".")) git2r::repository(here::here())  


fboehm/qtl2effects documentation built on Sept. 20, 2021, 7:34 p.m.