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.
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
) 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")
#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()
#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()
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)
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)