Code to generate Figure 8 of “Genetic architecture modulates diet induced hepatic mRNA and miRNA expression profiles”
Figure 8: Comparison of concordant eQTL pairs in mice fed different diets demonstrates similarities despite environment. Effect size of eQTL observed in both HFCA and HP diets are significantly correlated in both (A) cis-eQTL (r = 0.59) and (B) trans-eQTL (r = 0.50). Allele effects underlying concordant eQTL are highly correlated in both (C) cis-eQTL and (D) trans-eQTL despite different diet perturbations. Black dots represent eQTL with a significant allele-diet interaction (ANOVA, BH-adjusted p < 0.05). Color gradient represents dot density.
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)
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] stats graphics grDevices utils datasets methods base
##
## other attached packages:
## [1] epiR_1.0-14 survival_3.1-12 dplyr_0.8.5 devtools_2.3.0
## [5] usethis_1.6.0 dunn.test_1.3.5 reshape_0.8.8 ggpubr_0.2.5
## [9] magrittr_1.5 ggrepel_0.8.2 ggplot2_3.3.0 qtl2_0.21-5
##
## loaded via a namespace (and not attached):
## [1] tidyselect_1.0.0 xfun_0.13 remotes_2.1.1 purrr_0.3.4
## [5] lattice_0.20-41 splines_3.5.1 testthat_2.3.2 colorspace_1.4-1
## [9] vctrs_0.2.4 htmltools_0.4.0 yaml_2.2.1 blob_1.2.1
## [13] rlang_0.4.5 pkgbuild_1.0.6 pillar_1.4.3 glue_1.4.0
## [17] withr_2.2.0 DBI_1.1.0 bit64_0.9-7 BiasedUrn_1.07
## [21] sessioninfo_1.1.1 lifecycle_0.2.0 plyr_1.8.6 stringr_1.4.0
## [25] munsell_0.5.0 ggsignif_0.6.0 gtable_0.3.0 memoise_1.1.0
## [29] evaluate_0.14 knitr_1.28 callr_3.4.3 ps_1.3.2
## [33] parallel_3.5.1 fansi_0.4.1 Rcpp_1.0.4.6 backports_1.1.6
## [37] scales_1.1.0 desc_1.2.0 pkgload_1.0.2 fs_1.4.1
## [41] bit_1.1-15.2 digest_0.6.25 stringi_1.4.6 processx_3.4.2
## [45] rprojroot_1.3-2 grid_3.5.1 cli_2.0.2 tools_3.5.1
## [49] tibble_3.0.1 RSQLite_2.2.0 crayon_1.3.4 pkgconfig_2.0.3
## [53] Matrix_1.2-18 ellipsis_0.3.0 data.table_1.12.8 prettyunits_1.1.1
## [57] assertthat_0.2.1 rmarkdown_2.1 R6_2.4.1 compiler_3.5.1
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 8, this includes the congruent eQTL identified in the HFCA-fed and HP-fed mice models (matchDF
), the Best Linear Unbiased Predictions of the founder effects for the congruent eQT (matchDF_blup
), and the results of the allele-diet interaction study (allelDiet
).
# Master file
mRes <- readRDS(paste0("../", data.dir, "cluster_Inputs/masterRes_3_12_20.rds"))
# Congruent eQTL between HFCA and HP-fed mice models
matchDF <- readRDS(paste0("../", data.dir, "scan_Overlaps/HP_HFCAoverlap_KJ.rds"))
# BLUP results of the congruent eQTL between the HFCA and HP models
matchDF_blup <- readRDS(paste0("../", data.dir, "BLUP/blup_cors_HP_HFCA_V3.rds"))
# We want to highlight certain data points from the interactive/HFCA/HP allele-diet interaction study
allelDiet <- readRDS(paste0("../", data.dir, "allele_diet_Interaction/P_SNP_Diet_MatchMaker_Padjust2.rds"))
# From here, make a more manageable DF
sub <- matchDF %>% select(Symbol, CTStrict, A_qtlEffect, B_qtlEffect)
# And further subset by cis or trans
# Break up the CT:strict column
# Column A classification first. Just paying attention to the left part of the ":".
sub$A_qtl <- ifelse(sub$CTStrict == "cis:cis" | sub$CTStrict == "cis:trans", "cis", "trans")
# Column B classification second Just paying attention to the right part of the ":".
sub$B_qtl <- ifelse(sub$CTStrict == "cis:cis" | sub$CTStrict == "trans:cis", "cis", "trans")
# Use the new columns to subset accordingly
sub_cis <- sub[sub$A_qtl == "cis" | sub$B_qtl == "cis",] # 1801
sub_trans <- sub[sub$A_qtl == "trans" | sub$B_qtl == "trans",] # 135
# Note, the sum of of these dataframes equals 1936, 100 over the number of matches (1836) in matchDF. This is because there are 71 cis:trans and 29 trans:cis = 100 getting counted in both. This is likely due to resolution/power issues between the models, which we discuss in the manuscript.To handle this, we decided to only include transcripts with the same classification in both models.
# Test new method of only calling cis:cis cis and trans:trans trans
sub_cis <- sub[sub$A_qtl == "cis" & sub$B_qtl == "cis",] # 1701
sub_trans <- sub[sub$A_qtl == "trans" & sub$B_qtl == "trans",] # 35
# We want to highlight certain data points from the interactive/HFCA/HP allele-diet interaction study
# Cis color:
sub_cis_color <- sub_cis[sub_cis$Symbol %in% allelDiet$SYMBOL,] #color DF
# Trans color:
sub_trans_color <- sub_trans[sub_trans$Symbol %in% allelDiet$SYMBOL,] #color DF
Plot the relationship of (A) cis-eQTL and (B) trans-eQTL effect sizes (phenotypic variance explained) of congruent eQTL between the HFCA and HP-fed mice.
# ggplot cis
plot_cis <- ggplot(sub_cis, aes(A_qtlEffect, B_qtlEffect))
# add manuscript theme
plot_cis2 <- plot_cis +
geom_point(color = "dodgerblue", alpha = 0.6) +
stat_smooth(method = "lm",
formula = y ~ x,
size = 1,
color = "black",
fullrange = TRUE) +
theme_bw() +
labs(y ="HP Effect Size", x = "HFCA Effect Size")+
theme(text = element_text(size = 16),
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)) +
geom_point(data = sub_cis_color, aes(A_qtlEffect, B_qtlEffect), color = "black", alpha = 0.6) # this controls the black dots
# ggplot trans
plot_trans <- ggplot(sub_trans, aes(A_qtlEffect, B_qtlEffect))
# add manuscript theme
plot_trans2 <- plot_trans +
geom_point(color = "firebrick2", alpha = 0.6) +
stat_smooth(method = "lm",
formula = y ~ x,
size = 1,
color = "black",
fullrange = TRUE) +
theme_bw() +
labs(y ="HP Effect Size", x = "HFCA Effect Size")+
theme(text = element_text(size = 16),
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)) +
geom_point(data = sub_trans_color, aes(A_qtlEffect, B_qtlEffect), color = "black", alpha = 0.6) # this controls the black dots
# Save plots
out_effectSize <- ggarrange(plot_cis2, plot_trans2, ncol = 2, nrow = 1, widths = c(2,2))
out_effectSize
#ggsave(paste0("../",fig.dir,"figure10A.png"),out_effectSize,
# h = 5,
# w = 14,
# dpi = 600)
# Our goal is to visualize the relationship between distance of the peakSNP between congruent eQTL VS correlation of allele effect between diets. Therefore...
# ... we want to look at cis and trans separately
# Classification strategy- use CTStrict (4mb) and only consider cis=cis:cis and trans=trans:trans
# Subset
matchDF_cis_S <- matchDF_blup[matchDF_blup$CTStrict == "cis:cis",] # 1701
dim(matchDF_cis_S)
# Only want to include eQTL with clear trans regulation- trans:trans
matchDF_trans_S <- matchDF_blup[matchDF_blup$CTStrict == "trans:trans",] # 35
dim(matchDF_trans_S)
# Define color from allele-diet interaction study
# Cis color:
matchDF_cis_S_color <- matchDF_cis_S[matchDF_cis_S$Symbol %in% allelDiet$SYMBOL,] #color DF
# Trans color:
matchDF_trans_S_color <- matchDF_trans_S[matchDF_trans_S$Symbol %in% allelDiet$SYMBOL,] #color DF
Plot the relationship of the correlations of the founder alleles between congruent (C) cis-eQTL and (D) trans-eQTL in the HFCA and HP fed mice VS the distance of the peak SNPs
# Replot
# Cis
# ggplot + geom_point
plot_matchDF_cis_S <- ggplot(matchDF_cis_S, aes(x = Dist, blup_cors))
# add manuscript theme
plot_matchDF_cis_S2 <- plot_matchDF_cis_S +
geom_point(color = "dodgerblue", alpha = 0.3) +
theme_bw() +
labs(y ="Correlation", x = "Distance (Mb)")+
theme(text = element_text(size = 16),
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))) +
xlim(0,10) +
geom_point(data = matchDF_cis_S_color,
aes(x = Dist, blup_cors),
color = "black",
alpha = 0.3)
# Trans
# ggplot + geom_point
plot_matchDF_trans_S <- ggplot(matchDF_trans_S, aes(x = Dist, blup_cors))
# add manuscript theme
plot_matchDF_trans_S2 <- plot_matchDF_trans_S +
geom_point(color = "firebrick2", alpha = 0.3) +
theme_bw() +
labs(y ="Correlation", x = "Distance (Mb)")+
theme(text = element_text(size = 16),
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))) +
xlim(0,60) +
geom_point(data = matchDF_trans_S_color,
aes(x = Dist, blup_cors),
color = "black",
alpha = 0.3)
out_blup_S <- ggarrange(plot_matchDF_cis_S2, plot_matchDF_trans_S2, ncol = 2, nrow = 1, widths = c(2,2))
out_blup_S
#ggsave(paste0("../",fig.dir,"figure10B.png"),out_blup_S,
# h = 5,
# w = 14,
# dpi = 600)