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

Figure 4: Alleles from Wild Derived strains contribute to eQTL and mir-eQTL. (A) Center-scaled Best Linear Unbiased Predictor (BLUP) coefficients from 8 DO founder strains are taken at the peak SNP of all significant eQTL from the additive model. Their relative effect sizes are visualized and mapped to the eQTL’s physical location on a circularized mouse genome. The 6th and 7th tracks reveal a clear pattern of peak contribution (dark red and blue bands) from wild-type founder strains CAST/EiJ and PWK/PhJ that is consistent across the genome. (B) Boxplot representation of center-scaled BLUP coefficients from all significant eQTL in the additive model. Kruskal-Wallis rank sum test indicates significant differences (p < 2.2 x 10-16) among founder effects. Dunn’s post-hoc test confirms significantly (BH-adjusted p < 1.0 x 10-4) larger effects from the CAST/EiJ and PWK/PhJ strains relative to all other strains. (C) BLUP coefficients from 8 DO founder strains are taken at the peak SNP of all significant mirQTL from the additive model. Their relative effect sizes are visualized and mapped to the mirQTL’s physical location on a circularized mouse genome. Patterns of peak contribution by any particular strain is visually unclear. (D) Boxplot of center-scaled BLUP coefficients from all significant mirQTL in the additive model. Kruskal-Wallis rank sum test does not indicate significant differences (p < 0.064) among founder effects. Dunn’s post-hoc test shows a single significant (BH-adjusted p ≤ 0.05) difference with PWK/PhJ showing higher contributions over NOD/LtJ.

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) and all other necessary files. To generate Figure 4, this includes the circos plot dataframes (data, data_miRNA), the additive model BLUP results (normDat), and the circos plot’s angle values (offset)

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

# Plotting DF
data <- readRDS(paste0("../", data.dir, "circos_Plots/eQTL_CircosRdy.rds"))
data_miRNA <- readRDS(paste0("../", data.dir, "circos_Plots/mirQTL_CircosRdy.rds"))

# Additive mRNA BLUP results
#normDat <-  readRDS(paste0("../", data.dir, "circos_Plots/AddCovFP_BLUP_SNP.rds"))
#test
normDat <-  readRDS(paste0("../", data.dir, "BLUP/addBLUP_SNP_diet.rds"))
normDat <- normDat[,1:8]

# Circos Plot segAnglePo
load("../data/circos_Plots/offset.RData")

Plot Circos mRNA

  1. Circos Plot for the mRNA addBLUP SNPs
#Start with the outer ring
mouseRing <- data.frame(seg.name =c(paste0("chr",1:19), "chrX"),
                        seg.Start = offset$offset[-nrow(offset)],
                        seg.End =offset$offset[-1],
                        the.v = factor(NA),
                        NO = factor(NA),
                        stringsAsFactors = F)

#Probably doesn't play nice with absolute genomic position - retool the positions
mouseRing$seg.End <- mouseRing$seg.End - mouseRing$seg.Start
mouseRing$seg.Start <- 0
db <- segAnglePo(mouseRing, seg = mouseRing$seg.name)
colors <- rainbow(nrow(mouseRing), alpha = 0.5)

#pdf(paste0("../",fig.dir,"figure4A.pdf"))
par(mar = c(2,2,2,2))
plot(c(1,800), 
     c(1,800), 
     type="n", 
     axes=FALSE,
     xlab="", 
     ylab="",
     main="")

#Plot theouter ring, Radius 400, type = chromosome segments, W = width of the circle, scale = draw a scale
circos(R = 400, cir=db, 
       type="chr", 
       col="white", 
       print.chr.lab = T,
       W=4, 
       scale=TRUE)
#col6 through end - only works when you add the BLUPS
circos(R = 160, cir=db,
       mapping = data,
       col.v = 6,
       type="heatmap",
       col.bar = T,
       print.chr.lab = T,
       W=260, 
       scale=TRUE)

#dev.off()

Plot Circos miRNA

  1. Circos Plot for the miRNA addBLUP SNPs
#pdf(paste0("../",fig.dir,"figure4C.pdf"))
par(mar = c(2,2,2,2))
plot(c(1,800), 
     c(1,800), 
     type="n", 
     axes=FALSE,
     xlab="", 
     ylab="",
     main="")

#Plot theouter ring, Radius 400, type = chromosome segments, W = width of the circle, scale = draw a scale
circos(R = 400, cir=db, 
       type="chr", 
       col="white", 
       print.chr.lab = T,
       W=4, 
       scale=TRUE)
#col6 through end
circos(R = 160, cir=db,
       mapping = data_miRNA,
       col.v = 4,
       type="heatmap",
       col.bar = T,
       print.chr.lab = T,
       W=260, 
       scale=TRUE)

#dev.off()

Plot BLUPs mRNA

  1. Scaled-BLUP founder effects for mRNA
inDat <- normDat
founders <- c("A/J", "C57BL/6J", "129S1/SvImJ", "NOD/LtJ", "NZO/HlLtJ", "CAST/EiJ", "PWK/PhJ","WSB/EiJ")
names(inDat) <- founders
inDat <- t(scale(t(inDat)))

mtest <-  melt(inDat)
names(mtest)[2] <- "variable"
mtest$variable <- factor(mtest$variable, levels = founders)

p4 <- ggplot(mtest, aes(x=variable, y=abs(value), 
                        fill = variable, alpha = 0.95))+
  geom_boxplot()+
  scale_fill_manual(values=c("#FFDC00",
                             "#888888",
                             "#F08080",
                             "#0064C9",
                             "#7FDBFF",
                             "#2ECC40",
                             "#FF4136",
                             "#B10DC9"))+ 
  theme_bw()+
  theme(text = element_text(size = 16),
        axis.text.x = element_text(margin = margin(t = 5,
                                                   r = 0,
                                                   b = 0,
                                                   l = 0),
                                   angle = 45, 
                                   hjust = 1),
        axis.text.y = element_text(margin = margin(t = 0,
                                                   r = 0,
                                                   b = 0,
                                                   l = 5)), 
        legend.position = "none") + 
  labs(x = "Strain", 
       y = "Genome-wide\nScaled Effect Size")

p4

#ggsave(paste0("../",fig.dir,"figure4B.png"),p4,
#       h = 5,
#       w = 14,
#       dpi = 600)

Plot BLUP miRNA

  1. Scaled-BLUP founder effects for miRNA
names(data_miRNA)[4:11] <- founders
mtest <-  melt(data_miRNA[,c(4:11)])
## Using  as id variables
p4_mir <- ggplot(mtest, aes(x=variable, y=abs(value), 
                        fill = variable, alpha = 0.95))+
  geom_boxplot()+
  scale_fill_manual(values=c("#FFDC00",
                             "#888888",
                             "#F08080",
                             "#0064C9",
                             "#7FDBFF",
                             "#2ECC40",
                             "#FF4136",
                             "#B10DC9"))+ 
  theme_bw()+
  theme(text = element_text(size = 16),
        axis.text.x = element_text(margin = margin(t = 5,
                                                   r = 0,
                                                   b = 0,
                                                   l = 0),
                                   angle = 45, 
                                   hjust = 1),
        axis.text.y = element_text(margin = margin(t = 0,
                                                   r = 0,
                                                   b = 0,
                                                   l = 5)), 
        legend.position = "none") + 
  labs(x = "Strain", 
       y = "Genome-wide\nScaled Effect Size")

p4_mir

#ggsave(paste0("../",fig.dir,"figure4D.png"),p4_mir,
#      h = 5,
#       w = 14,
#       dpi = 600)