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.

Set up

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 Data

Load in the master results file (mRes).

# Master file
mRes <- readRDS(paste0("../", data.dir, "cluster_Inputs/masterRes_3_12_20.rds"))

Plot

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)