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

Figure 5: Effect size of eQTL is lower in genes and miRNA that are differentially expressed. The allelic effect size is significantly lower (p < 2.22x10-16) in genes that are differentially expressed in both cis and trans-eQTL. (A) For cis-eQTL, the median variance explained when the transcript is differentially expressed is 0.186 compared to 0.238 when it is not differentially expressed. (B) For trans-eQTL, the median variance explained when the transcript is differentially expressed is 0.124 compared to 0.143 when it is not differentially expressed. (C) The allelic effect size is not significantly different in cis or trans-mirQTL; however, the median of differentially expressed cis-mirQTL is 0.139 compared to 0.169 when not differentially expressed. (D) Similarly, differentially expressed mirQTL in trans had a smaller median than those not differentially expressed at 0.127 and 0.141, respectively.

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)

getwd()

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 Data

Load in the master results file (mRes) and all other necessary files. For the differential gene expression analysis, this includes the list of transcripts that pass the RMA DEG threshold (expr_filt) and the differential expression results (mRNA_DE). For the differential miRNA expression analysis, this includes the miRNA cross2 object (mirCross) and the differential expression results (filtDF).

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

# miRNA cross object
mirCross <- readRDS(paste0("../", data.dir, "cluster_Inputs/mirCross.rds"))
mirExpr <- mirCross$pheno

# Transcripts that pass the RMA threshold (4 or greater in at least 20% samples) and miRNA  were used in the differential expression analysis
# mRNA
expr_filt <- readRDS(paste0("../", data.dir,"differential_Expression/mRNA_passThresh.rds"))
expr_filt.t <- as.data.frame(as.matrix(t(expr_filt)))
# miRNA
mirExpr_filt.t <- as.data.frame(as.matrix(t(mirExpr)))

# Differential Expression Results
# mRNA
mRNA_DE <- readRDS(paste0("../", data.dir,"differential_Expression/mRNA_diffs_filtered.rds"))
# miRNA
filtDF <- readRDS(paste0("../", data.dir,"differential_Expression/miRNA_diffs_filtered_246.rds"))

Prepare mRNA data

Get mRNA data ready for plotting

# Make a new column stating if the transcript has an eQTL
expr_filt.t$eQTL_add <- ifelse(rownames(expr_filt.t) %in% mRes$addCov_FP$lodcolumn,
                         "YES","NO")
# We already ran the DEG analysis. Lets load those results and make a column in expr_filt.t stating if there is a DEG!
expr_filt.t$hasDEG <- ifelse(rownames(expr_filt.t) %in% mRNA_DE$transcript,
                         "YES","NO")
# Check the structure of hasDEG. We want factor and yes = 1, no = 2
str(expr_filt.t$hasDEG)
expr_filt.t$hasDEG <- factor(expr_filt.t$hasDEG, levels = c("YES", "NO"))
str(expr_filt.t$eQTL_add)
expr_filt.t$eQTL_add <- factor(expr_filt.t$eQTL_add, levels = c("YES", "NO"))

# Simplify the dataframe 
df <- expr_filt.t %>% select(eQTL_add, hasDEG)

# Filter to transcripts with an eQTL only
df_eQTL <- df[df$eQTL_add == "YES",] # 4794

# Now add if the eQTL is cis or trans
# Make a transcript column to make the merge easier
df_eQTL$transcript <- row.names(df_eQTL)

# Note, there are ~300 multimappers in addcovar model. To avoid problems merging, use left join and merge the DEG results to the add object. 
add <- mRes$addCov_FP #6692
add_DEG <- left_join(add, df_eQTL, by = c("lodcolumn" = "transcript"))

# Filter rows with no DEG data (we only want to include transcripts that passed our filtering standards)
add_DEG2 <- add_DEG[complete.cases(add_DEG),] #5079

# Subset by cis/trans
#cis
add_DEG2_cis <- add_DEG2[add_DEG2$CT_strict == "cis",] # 4252
#trans
add_DEG2_trans <- add_DEG2[add_DEG2$CT_strict == "trans",] #827

Prepare miRNA data

Get miRNA DE data ready for plotting

# Make a new column stating if the transcript has a mirQTL
mirExpr_filt.t$mirQTL_add <- ifelse(rownames(mirExpr_filt.t) %in% mRes$mirAddCov$lodcolumn,
                         "YES","NO")
# We already ran the DEG analysis. Lets load those results and make a column in mirExpr_filt.t stating if there is a DEG!
mirExpr_filt.t$hasDEmiRNA <- ifelse(rownames(mirExpr_filt.t) %in% filtDF$transcript,
                         "YES","NO")
# Check the structure of hasDEmiRNA. We want factor and yes = 1, no = 2
str(mirExpr_filt.t$hasDEmiRNA)
mirExpr_filt.t$hasDEmiRNA <- factor(mirExpr_filt.t$hasDEmiRNA, levels = c("YES", "NO"))
str(mirExpr_filt.t$mirQTL_add)
mirExpr_filt.t$mirQTL_add <- factor(mirExpr_filt.t$mirQTL_add, levels = c("YES", "NO"))

# Simplify the dataframe 
df_mir <- mirExpr_filt.t %>% select(mirQTL_add, hasDEmiRNA)

# Filter to transcripts with an eQTL only
df_mirQTL <- df_mir[df_mir$mirQTL_add == "YES",] # 33

# Now add if the eQTL is cis or trans
# Make a transcript column to make the merge easier
df_mirQTL$transcript <- row.names(df_mirQTL)

# Note, there are ~300 multimappers in addcovar model. To avoid problems merging, use left join and merge the DEG results to the add object. 
mir_add <- mRes$mirAddCov #37
mir_add_DEG <- left_join(mir_add, df_mirQTL, by = c("lodcolumn" = "transcript"))

# Filter for rows with no DEG data (we only want to include transcripts that passed our filtering standards)
mir_add_DEG2 <- mir_add_DEG[complete.cases(mir_add_DEG),] #37

# Subset by cis/trans
#cis
mir_add_DEG2_cis <- mir_add_DEG2[mir_add_DEG2$CT == "cis",]
#trans
mir_add_DEG2_trans <- mir_add_DEG2[mir_add_DEG2$CT == "trans",]

Plot

A-D) Plot mRNA and miRNA DE results X phenotypic variance explained

### mRNA
# ggplot + geom_boxplot
plot_cis <- ggplot(add_DEG2_cis, 
                               aes(x = hasDEG, 
                                   y = qtlEffect)) +
  geom_boxplot(fill = "dodgerblue", alpha = 0.8)

# Add theme for manuscript
plot_cis2 <- plot_cis + 
  theme_bw() + 
  theme(text = element_text(size = 12),
        axis.text.x = element_text(margin = margin(t = 3,
                                                   r = 0,
                                                   b = 0,
                                                   l = 0)),
        axis.text.y = element_text(margin = margin(t = 0,
                                                   r = 0,
                                                   b = 0,
                                                   l = 3))) + 
  ylim(0,1)+
  scale_fill_manual(values = c("dodgerblue"))+
  labs(y ="Phenotypic variance explained ", x = "DEG Status") + 
  scale_x_discrete(labels=c("DEG (yes)" =
                        paste0("DEG\n(n=",
                               table(add_DEG2_cis$hasDEG)[[1]],")"), #3131
                            "DEG (no)" =
                    paste0("DEG\n(n=",table(add_DEG2_cis$hasDEG)[[2]],")"))) + #1121
  stat_compare_means(comparisons = list(comp1 = c(1,2)),
                     method = "wilcox", 
                     label = "p.format",
                     paired = FALSE)+ 
  guides(fill = FALSE)

# View
#plot_cis2

# Now repeat plot for trans
# ggplot + geom_boxplot
plot_trans <- ggplot(add_DEG2_trans, 
                               aes(x = hasDEG, 
                                   y = qtlEffect)) +
  geom_boxplot(fill = "firebrick2", alpha = 0.8)

# Add theme for manuscript
plot_trans2 <- plot_trans + 
  theme_bw() + 
  theme(text = element_text(size = 12),
        axis.text.x = element_text(margin = margin(t = 3,
                                                   r = 0,
                                                   b = 0,
                                                   l = 0)),
        axis.text.y = element_text(margin = margin(t = 0,
                                                   r = 0,
                                                   b = 0,
                                                   l = 3))) + 
  ylim(0,1)+
  scale_fill_manual(values = c("dodgerblue","firebrick2"))+
  labs(y ="Phenotypic variance explained ", x = "DEG Status")+
  scale_x_discrete(labels=c("DEG (yes)" =
                        paste0("DEG\n(n=",
                               table(add_DEG2_cis$hasDEG)[[1]],")"), 
                            "DEG (no)" =
                    paste0("DEG\n(n=",table(add_DEG2_cis$hasDEG)[[2]],")")))+
  stat_compare_means(comparisons = list(comp1 = c(1,2)),
                     method = "wilcox", 
                     label = "p.format",
                     paired = FALSE)+ 
  guides(fill = FALSE)
# View
#plot_trans2

### miRNA
# Cis
# ggplot + geom_boxplot
plot_cis_mir <- ggplot(mir_add_DEG2_cis, 
                               aes(x = hasDEmiRNA, 
                                   y = qtlEffect)) +
  geom_boxplot(fill = "dodgerblue", alpha = 0.8)

# Add theme for manuscript
plot_cis_mir2 <- plot_cis_mir + 
  theme_bw() + 
  theme(text = element_text(size = 12),
        axis.text.x = element_text(margin = margin(t = 3,
                                                   r = 0,
                                                   b = 0,
                                                   l = 0)),
        axis.text.y = element_text(margin = margin(t = 0,
                                                   r = 0,
                                                   b = 0,
                                                   l = 3))) + 
  ylim(0,1)+
  scale_fill_manual(values = c("dodgerblue"))+
  labs(y ="Phenotypic variance explained ", x = "DEmiRNA Status") + 
  scale_x_discrete(labels=c("DEmiRNA (yes)" =
                        paste0("DEmiRNA\n(n=",
                               table(mir_add_DEG2_cis$hasDEmiRNA)[[1]],")"), #18
                            "DEG (no)" =
                        paste0("DEG\n(n=",table(mir_add_DEG2_cis$hasDEmiRNA)[[2]],")"))) + #19
  stat_compare_means(comparisons = list(comp1 = c(1,2)),
                     method = "wilcox", 
                     label = "p.format",
                     paired = FALSE)+ 
  guides(fill = FALSE)

# View
#plot_cis_mir2

####
# Trans
plot_trans_mir <- ggplot(mir_add_DEG2_trans, 
                               aes(x = hasDEmiRNA, 
                                   y = qtlEffect)) +
  geom_boxplot(fill = "firebrick2", alpha = 0.8)

# Add theme for manuscript
plot_trans_mir2 <- plot_trans_mir + 
  theme_bw() + 
  theme(text = element_text(size = 12),
        axis.text.x = element_text(margin = margin(t = 3,
                                                   r = 0,
                                                   b = 0,
                                                   l = 0)),
        axis.text.y = element_text(margin = margin(t = 0,
                                                   r = 0,
                                                   b = 0,
                                                   l = 3))) + 
  ylim(0,1)+
  scale_fill_manual(values = c("dodgerblue"))+
  labs(y ="Phenotypic variance explained ", x = "DEmiRNA Status") + 
  scale_x_discrete(labels=c("DEmiRNA (yes)" =
                        paste0("DEmiRNA\n(n=",
                               table(mir_add_DEG2_cis$hasDEmiRNA)[[1]],")"), #18
                            "DEG (no)" =
                        paste0("DEG\n(n=",table(mir_add_DEG2_cis$hasDEmiRNA)[[2]],")"))) + #19
  stat_compare_means(comparisons = list(comp1 = c(1,2)),
                     method = "wilcox", 
                     label = "p.format",
                     paired = FALSE)+ 
  guides(fill = FALSE)

# View
#plot_trans_mir2

# Save plots
# Make a large plot that includes mRNA and miRNA
out_DEG_effectSize_all <- ggarrange(plot_cis2,plot_cis_mir2, plot_trans2,plot_trans_mir2, ncol = 2, nrow = 2, widths = c(2,2))
out_DEG_effectSize_all

#ggsave(paste0("../",fig.dir,"figure7.png"),out_DEG_effectSize_all,
#       h = 8,
#       w = 12,
#       dpi = 600)