Code to generate Figure 2 of “Genetic architecture modulates diet induced hepatic mRNA and miRNA expression profiles

Figure 2: High-resolution hepatic eQTL and mirQTL in Diversity Outbred Mice demonstrate complex regulation of expression traits. Results from expression quantitative trait loci (eQTL) analyses using a genome scan model treating the high-fat cholic acid diet (HFCA) and high-protein diet (HP) groups as an additive covariate. (A) An eQTL plot visualizes all significant associations between gene expression and structural SNP variants from a genome scan. The absolute genomic positions of the SNP and gene transcript probe set in the mouse genome are shown on the X and Y-axes, respectively. The blue line represents cis-eQTL, positions where gene expression variance is associated with a proximal (±4 Mb) SNP variant. (B) Resolution of cis-eQTL is determined by the distance from the SNP to the transcript probe set (x axis) – resolution across all 5,460 eQTL-probe pairs on the same chromosome shows 95.4% of cis-eQTL are within ±4 Mb of their probe set and tend to have highly significant log of the odds scores (LOD – y axis). Cis-eQTL resolution in the DO is estimated to be 0.32 Mb. (C) Narrow-sense heritability is calculated for the probe set expression of each significant eQTL in the additive QTL analysis. (D) mirQTL plot of the results from the genome scan of miRNA transcription data with diet as an additive covariate. (E) Resolution of 18 cis-mirQTL, showing 100% of cis-mirQTL are within ±4 Mb. Cis-mirQTL resolution in the DO is estimated to be 0.25 Mb. (F) Distribution of narrow-sense heritability for the miRNA expressions with significant mirQTL in the additive genome scan.

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 
#setwd("/Users/kristenjames/Box/DO2_figshare/markdown/")

data.dir <-  "data/"
results.dir <-  "results/"
fig.dir <- "figures/"
script.dir <- "scripts/" 

# Source the CT_plotter function
source(paste0("../", script.dir, "CT_plotter.R"))
source(paste0("../", script.dir, "CT_plotter_miRNA.R"))
# Source the Precision_plotter function
source(paste0("../", script.dir, "Precision_plotter.R"))
source(paste0("../", script.dir, "Precision_plotter_miRNA.R"))

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) and all other necessary files.

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

Plot- Trans-band

A) Whole-genome plot with trans-bands mRNA

#png(paste0("../",fig.dir,"figure2A.png"), width = 16, height = 16, units= "in" , res = 500)
CTplot(mRes$addCov_FP, ptSize = 1, fntSize = 20, xlab_shift = 10, ylab_shift = 10)

#dev.off()

Plot- Precision

B) Precision x LOD, mRNA

#png(paste0("../",fig.dir,"figure2B.png"), width = 16, height = 12, units= "in" , res = 500)
precPlot(mRes$addCov_FP, run = "mRNA Additive", fntSize = 12, xlab_shift = 5, ylab_shift = 5)
## Warning: Removed 1 rows containing missing values (geom_point).

#dev.off()

Plot- Heritability

C) Narros-sense heritability x Number of eQTL histogram

mRNA_heritHisto <- ggplot(mRes$addCov_FP, aes(x = hsq)) +
  geom_histogram(binwidth = 0.01, color="black", fill="grey")

plotHist <- mRNA_heritHisto + 
  theme_bw() + 
  theme(text = element_text(size = 16),
        axis.title.x = element_text(
          margin = margin(t = 5,
                          r = 0, 
                          b = 0, 
                          l = 0)),
        axis.title.y = element_text(
          margin = margin(t = 0, 
                          r = 5,
                          b = 0, 
                          l = 0)))+
  labs(y ="Number of eQTL", 
       x = "Narrow-sense heritability")

#png(paste0("../",fig.dir,"figure2C.png"), width = 16, height = 12, units= "in" , res = 500)
plotHist

#dev.off()

Plot- Trans-band miRNA

D) Whole-genome plot with trans-bands miRNA

#png(paste0("../",fig.dir,"figure2D.png"), width = 16, height = 16, units= "in" , res = 500)
CTplot_mir(mRes$mirAddCov, ptSize = 1, fntSize = 20, xlab_shift = 10, ylab_shift = 10)

#dev.off()

Plot- Precision miRNA

E) Precision x LOD, miRNA

#png(paste0("../",fig.dir,"figure2E.png"), width = 16, height = 12, units= "in" , res = 500)
precPlot_miRNA(mRes$mirAddCov, run = "miRNA Additive", fntSize = 12, xlab_shift = 5, ylab_shift = 5)

#dev.off()

Plot- Heritability, miRNA

F) Narrow-sense heritability x Number of mirQTL histogram

miRNA_heritHisto <- ggplot(mRes$mirAddCov, aes(x = hsq)) +
  geom_histogram(binwidth = 0.01, color="black", fill="grey")

plotHist <- miRNA_heritHisto + 
  theme_bw() + 
  theme(text = element_text(size = 16),
        axis.title.x = element_text(
          margin = margin(t = 5,
                          r = 0, 
                          b = 0, 
                          l = 0)),
        axis.title.y = element_text(
          margin = margin(t = 0, 
                          r = 5,
                          b = 0, 
                          l = 0)))+
  labs(y ="Number of mirQTL", 
       x = "Narrow-sense heritability")

#png(paste0("../",fig.dir,"figure2F.png"), width = 16, height = 12, units= "in" , res = 500)
plotHist

#dev.off()