Code to generate Figure 6 of “Genetic architecture modulates diet induced hepatic mRNA and miRNA expression profiles”
Figure 6: Magnitude of correlation is similar between miRNA and mRNA with cis, trans, and no-eQTL and reduced when mRNA are not differentially expressed. Expression data was used to generate correlations between pairs of miRNA (n = 246) and mRNA (n = 24,004) and the results were filtered for FDR-corrected significance (α ≤ 0.05) resulting in 1,424,040 significant correlations. miRNA-mRNA pairs were classified by their eQTL status in the additive model. 1,055,010 eQTL with transcripts that passed the quality threshold for the DE analysis were included. (A) Spearman’s rho was plotted against the DE status of the mRNA from the additive genome scan. Median ± standard error (SE) for DEG and non-DEG include 0.245 ± 0.0001 and 0.207 ± 0.0002, respectively. (B) Spearman’s rho was plotted against the eQTL status of mRNA that were differentially expressed by diet and that were (C) not differentially expressed by diet. For differentially expressed mRNA with cis, trans, and no-eQTL the median ± SE are 0.238 ± 0.0002, 0.245 ± 0.0005 , 0.249 ± 0.001, respectively. For non-differentially expressed mRNA with cis, trans, and no-eQTL the median ± SE are 0.207 ± 0.003, 0.207 ± 0.0007, 0.207 ± 0.002.
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 differential expression results for mRNA (mRNA_DE
), the dataframe containing mRNA that passed the RMA threshold and were included in the differential expression analysis (expr_filt.t
), and the miRNA-mRNA correlation results (res
).
# Master file
mRes <- readRDS(paste0("../", data.dir, "cluster_Inputs/masterRes_3_12_20.rds"))
# Load mRNA-miRNA correlation results
res <- readRDS(paste0("../", data.dir,"correlation/DO2_mesMirLong_kj.rds")) # note, res has multimappers removed; object with the Rho values we need for this analysis
# Factor this as desired for plotting and stats in mind
str(res$status_Add)
res$status_Add <- factor(res$status_Add, levels = c("cis","trans", "No eQTL"))
# Load differentially expressed genes results. We will use them downstream.
mRNA_DE <- readRDS(paste0("../", data.dir,"differential_Expression/mRNA_diffs_filtered.rds"))
# Load expression DF
expr_filt <- readRDS(paste0("../", data.dir,"differential_Expression/mRNA_passThresh.rds"))
# Transpose
expr_filt.t <- as.data.frame(as.matrix(t(expr_filt))) # we want transcripts to be rows
# We already ran the DEG analysis. Use those results to 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")
expr_filt.t$hasDEG <- factor(expr_filt.t$hasDEG, levels = c("YES", "NO"))
# Great, now we need to add a "hasDEG" column to the res object from above.
# Recall, expr_filt.t is our list of DEG that pass the RMA filtering threshold
tmp <- expr_filt.t %>% select(hasDEG) # this contains all transcripts that pass RMA thresh
# Make a rowname with transcript ID for easier merging
tmp$transcript <- row.names(tmp)
# left join
res2 <- left_join(res, tmp, by = c("mRNA" = "transcript"))
dim(res2) #1424040 x 7
res3 <- res2[complete.cases(res2),]
dim(res3) #1055010 x 7
# Get ready to plot
# DE Yes
cor_DEyes <- res3[res3$hasDEG == "YES",]
# DE No
cor_DEno <- res3[res3$hasDEG == "NO",]
Generate plot containing 3 panes side by side by side.
A) DE-yes/DE-no by Spearman’s rho
B) cis/trans/none by Spearman’s rho (DE yes)
C) cis/trans/none x by Spearman’s rho (non-DE)
# Plot
# Has DEG and No DEG plot
fig_DEGrho <- ggplot(res3, aes(x=hasDEG, y=abs(rho), fill = hasDEG))+ # Use res3 (complete cases) for DEG analysis.
geom_boxplot(fill = c("gray60","gray96")) +
stat_compare_means(method = "t.test",
label = "p.format",
label.x.npc = "center",
paired = FALSE,
size = 3) +
theme_bw() +
labs(y = "Spearman's Rho (Absolute Value)", x = "DEG Status") +
scale_x_discrete(limits = c("YES", "NO"),
labels=c("YES" = paste0("Yes\n(n=",table(res3$hasDEG)[[1]],")"),
"NO" = paste0("No\n(n=",table(res3$hasDEG)[[2]],")"))) +
theme(axis.text = element_text(size = 11),
axis.title.x = element_text(margin = margin(t = 10),
size = 11),
axis.title.y = element_blank(),
legend.position="none") +
ylim(0.1,1.2)
# DE yes
fig_DEyes <- ggplot(cor_DEyes, aes(x=status_Add, y=abs(rho), fill = status_Add))+
geom_boxplot(fill = c("dodgerblue","firebrick2","#00BA38"))+
stat_compare_means(comparisons = list(comp1 = c(1,2),
comp2 = c(1,3),
comp3 = c(2,3)),
method = "wilcox",
label = "p.format",
paired = FALSE,
size = 3)+
theme_bw()+
labs(y = "Spearman's Rho (Absolute Value)",x = "DEG")+
scale_x_discrete(limits = c("cis", "trans", "No eQTL"),
labels=c("cis" = paste0("Cis\n(n=",table(cor_DEyes$status_Add)[[1]],")"),
"trans" = paste0("Trans\n(n=",table(cor_DEyes$status_Add)[[2]],")"),
"No eQTL" = paste0("No eQTL\n(n=",table(cor_DEyes$status_Add)[[3]],")")))+
theme(axis.text = element_text(size = 11),
axis.title.x = element_text(margin = margin(t = 10),
size = 11),
axis.title.y = element_blank(),
legend.position="none") +
ylim(0.1,1.2)
# DE no
fig_DEno <- ggplot(cor_DEno, aes(x=status_Add, y=abs(rho), fill = status_Add))+
geom_boxplot(fill = c("dodgerblue","firebrick2","#00BA38"))+
stat_compare_means(comparisons = list(comp1 = c(1,2),
comp2 = c(1,3),
comp3 = c(2,3)),
method = "wilcox",
label = "p.format",
paired = FALSE,
size = 3)+
theme_bw()+
labs(y = "Spearman's Rho (Absolute Value)",x = "non-DEG")+
scale_x_discrete(limits = c("cis", "trans", "No eQTL"),
labels=c("cis" = paste0("Cis\n(n=",table(cor_DEno$status_Add)[[1]],")"),
"trans" = paste0("Trans\n(n=",table(cor_DEno$status_Add)[[2]],")"),
"No eQTL" = paste0("No eQTL\n(n=",table(cor_DEno$status_Add)[[3]],")")))+
theme(axis.text = element_text(size = 11),
axis.title.x = element_text(margin = margin(t = 10),
size = 11),
axis.title.y = element_blank(),
legend.position="none") +
ylim(0.1,1.2)
# Arrange
DEandRhoPlots2 <- ggarrange(fig_DEGrho, fig_DEyes, fig_DEno, ncol = 3, nrow = 1, widths = c(2,2,2))
# View
DEandRhoPlots2
# Save
#ggsave(paste0("../",fig.dir,"figure8.png"),DEandRhoPlots2,
# h = 6,
# w = 9,
# dpi = 600)