Code to generate Figure 3 of “Genetic architecture modulates diet induced hepatic mRNA and miRNA expression profiles”
Figure 3: eQTL and mirQTL differ in heritability and effect size in Diversity Outbred mice. (A, B) Narrow-sense heritability scores and R2 differences between full and null Haley-Knott regression models (phenotypic variance explained) were calculated for every significant eQTL and mirQTL in the additive genome scan models and colored by their LOD score. Polynomial and simple regression models were fit to the mRNA and miRNA data, respectively, to understand the association between heritability (x) and phenotypic variance explained by the eQTL (y). (C, D) Wilcoxon rank sum tests reveal a significant difference in the phenotypic variance explained between cis and trans-eQTL and a similar but insignificant (p ≥ 0.05) difference in the phenotypic variance explained between the cis and trans-mirQTL.
Load libraries and set working directories
library(qtl2)
library(ggplot2)
library(ggrepel)
library(ggpubr)
library(reshape)
library(dunn.test)
library(devtools)
library(dplyr)
library(epiR)
library(tidyr)
library(plyr)
library(scales)
library(OmicCircos)
getwd() # paths set such that wd is in markdown folder
data.dir <- "data/"
results.dir <- "results/"
fig.dir <- "figures/"
script.dir <- "scripts"
Version of R and packages used:
sessionInfo()
## R version 3.5.1 (2018-07-02)
## Platform: x86_64-apple-darwin15.6.0 (64-bit)
## Running under: macOS 10.15.5
##
## Matrix products: default
## BLAS: /Library/Frameworks/R.framework/Versions/3.5/Resources/lib/libRblas.0.dylib
## LAPACK: /Library/Frameworks/R.framework/Versions/3.5/Resources/lib/libRlapack.dylib
##
## locale:
## [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
##
## attached base packages:
## [1] parallel stats4 stats graphics grDevices utils datasets
## [8] methods base
##
## other attached packages:
## [1] OmicCircos_1.20.0 GenomicRanges_1.34.0 GenomeInfoDb_1.18.2
## [4] IRanges_2.16.0 S4Vectors_0.20.1 BiocGenerics_0.28.0
## [7] scales_1.1.0 plyr_1.8.6 tidyr_1.0.2
## [10] epiR_1.0-14 survival_3.1-12 dplyr_0.8.5
## [13] devtools_2.3.0 usethis_1.6.0 dunn.test_1.3.5
## [16] reshape_0.8.8 ggpubr_0.2.5 magrittr_1.5
## [19] ggrepel_0.8.2 ggplot2_3.3.0 qtl2_0.21-5
##
## loaded via a namespace (and not attached):
## [1] pkgload_1.0.2 bit64_0.9-7 splines_3.5.1
## [4] assertthat_0.2.1 blob_1.2.1 GenomeInfoDbData_1.2.0
## [7] yaml_2.2.1 remotes_2.1.1 sessioninfo_1.1.1
## [10] pillar_1.4.3 RSQLite_2.2.0 backports_1.1.6
## [13] lattice_0.20-41 glue_1.4.0 digest_0.6.25
## [16] ggsignif_0.6.0 XVector_0.22.0 colorspace_1.4-1
## [19] htmltools_0.4.0 Matrix_1.2-18 pkgconfig_2.0.3
## [22] zlibbioc_1.28.0 purrr_0.3.4 processx_3.4.2
## [25] tibble_3.0.1 ellipsis_0.3.0 withr_2.2.0
## [28] cli_2.0.2 crayon_1.3.4 memoise_1.1.0
## [31] evaluate_0.14 ps_1.3.2 fs_1.4.1
## [34] fansi_0.4.1 pkgbuild_1.0.6 tools_3.5.1
## [37] data.table_1.12.8 prettyunits_1.1.1 lifecycle_0.2.0
## [40] stringr_1.4.0 munsell_0.5.0 callr_3.4.3
## [43] compiler_3.5.1 rlang_0.4.5 grid_3.5.1
## [46] RCurl_1.98-1.2 bitops_1.0-6 rmarkdown_2.1
## [49] testthat_2.3.2 gtable_0.3.0 DBI_1.1.0
## [52] R6_2.4.1 knitr_1.28 bit_1.1-15.2
## [55] rprojroot_1.3-2 desc_1.2.0 stringi_1.4.6
## [58] BiasedUrn_1.07 Rcpp_1.0.4.6 vctrs_0.2.4
## [61] tidyselect_1.0.0 xfun_0.13
The data used in these analyses is available to download at Dryad.
Load in the master results file (mRes
).
# Master file
mRes <- readRDS(paste0("../", data.dir, "cluster_Inputs/masterRes_3_12_20.rds"))
Generate all plots (A-D) in Figure 3. This includes the heritability x phenotypic variance explained scatter plot as well as the cis/trans status x phenotypic variance explained boxplot for mRNA and miRNA.
mRNA_heritPhenVar <- ggplot(mRes$addCov_FP, aes(hsq, qtlEffect, colour = lod))
a1 <- mRNA_heritPhenVar +
geom_point() +
stat_smooth(method = "lm", formula = y ~ x + I(x^2), size = 1,
color = "firebrick3") +
theme_bw() +
labs(y ="Phenotypic variance explained",
x = "Narrow-sense (additive) heritability",
colour = "LOD")+
theme(text = element_text(size = 12),
axis.text.x = element_text(margin = margin(t = 5,
r = 0,
b = 0,
l = 0)),
axis.text.y = element_text(margin = margin(t = 0,
r = 0,
b = 0,
l = 5)))+
scale_x_continuous(limits = c(0,1)) +
scale_y_continuous(limits = c(0,1))
mRNA_phenVarBox <- ggplot(mRes$addCov_FP, aes(x = CT_strict, y = qtlEffect, fill = CT_strict)) +
geom_boxplot(alpha = 0.8)
a2 <- mRNA_phenVarBox +
theme_bw() +
theme(text = element_text(size = 12),
axis.text.x = element_text(margin = margin(t = 5,
r = 0,
b = 0,
l = 0)),
axis.text.y = element_text(margin = margin(t = 0,
r = 0,
b = 0,
l = 5))) +
ylim(0,1)+
scale_fill_manual(values = c("dodgerblue","firebrick2"))+
labs(y ="Phenotypic variance explained ", x = "eQTL Status")+
scale_x_discrete(labels=c("cis" = paste0("Cis\n(n=",table(mRes$addCov_FP$CT_strict)[[1]],")"),
"trans" = paste0("Trans\n(n=",table(mRes$addCov_FP$CT_strict)[[2]],")")))+
stat_compare_means(comparisons = list(comp1 = c(1,2)),
method = "wilcox",
label = "p.format",
paired = FALSE)+
guides(fill = FALSE)
mir_heritPhenVar <- ggplot(mRes$mirAddCov, aes(hsq, qtlEffect, colour = lod))
b1 <- mir_heritPhenVar +
geom_point() +
stat_smooth(method = "lm",
formula = y ~ x,
size = 1,
color = "firebrick3") +
theme_bw() +
labs(y ="Phenotypic variance explained", x = "Narrow-sense (additive) heritability", colour = "LOD")+
theme(text = element_text(size = 12),
axis.text.x = element_text(margin = margin(t = 5,
r = 0,
b = 0,
l = 0)),
axis.text.y = element_text(margin = margin(t = 0,
r = 0,
b = 0,
l = 5)))+
scale_x_continuous(limits = c(0,1)) +
scale_y_continuous(limits = c(0,1))
mir_phenVarBox <- ggplot(mRes$mirAddCov, aes(x = CT, y = qtlEffect, fill = CT)) +
geom_boxplot(alpha = 0.8)
b2 <- mir_phenVarBox +
theme_bw() +
theme(text = element_text(size = 12),
axis.text.x = element_text(margin = margin(t = 5,
r = 0,
b = 0,
l = 0)),
axis.text.y = element_text(margin = margin(t = 0,
r = 0,
b = 0,
l = 5))) +
ylim(0,1)+
scale_fill_manual(values = c("dodgerblue","firebrick2"))+
labs(y ="Phenotypic variance explained ", x = "eQTL Status")+
scale_x_discrete(labels=c("cis" = paste0("Cis\n(n=",table(mRes$mirAddCov$CT)[[1]],")"),
"trans" = paste0("Trans\n(n=",table(mRes$mirAddCov$CT)[[2]],")")))+
stat_compare_means(comparisons = list(comp1 = c(1,2)),
method = "wilcox",
label = "p.format",
paired = FALSE)+
guides(fill = FALSE)
plot <- ggarrange(a1,a2,b1,b2, ncol = 2, nrow = 2, widths = c(2,1))
plot
#ggsave(paste0("../", fig.dir, "figure3.png"), plot,
# h = 6,
# w = 12,
# dpi = 600)