Code to generate all figures in the Supplemental Material for, “Genetic architecture modulates diet induced hepatic mRNA and miRNA expression profiles”
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)
library(VennDiagram)
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()
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 the Supplemental Figures, this includes the correlation results (res
), the eQTL Best Linearized Unbiased Predictor (BLUP
) results by the HFCA diet (HFCABLUP_SNP
) and the HP diet (HPBLUP_SNP
), as well as the mirQTL BLUP results for the HFCA diet (mirCoef
) and HP diet (mirCoefHP
). Add_Int_overlap
describes the concordant eQTL between the additive and interactive models and mirAdd_Int_overlap
describes the concordant mirQTL. CT_plotter
and CT_plotter_miRNA
are functions that generate global expression plots.
# Master file
mRes <- readRDS(paste0("../", data.dir, "cluster_Inputs/masterRes_3_12_20.rds"))
# Correlation results (SFig2)
res <- readRDS(paste0("../", data.dir, "correlation/DO2_mesMirLong_kj.rds"))
# HFCA and HP eQTL BLUP results (SFig6)
HFCABLUP_SNP <- readRDS(paste0("../", data.dir, "BLUP/hfcaBLUP_SNP.rds"))
HFCABLUP_SNP <- HFCABLUP_SNP[,1:8] # get rid of intercept colummn (column 9)
HPBLUP_SNP <- readRDS(paste0("../", data.dir, "BLUP/hpBLUP_SNP.rds"))
HPBLUP_SNP <- HPBLUP_SNP[,1:8] # get rid of intercept colummn (column 9)
# HFCA and HP mirQTL BLUP results (SFig6)
mirCoef <- readRDS(paste0("../", data.dir, "BLUP/mirHFCABLUP_SNP.rds"))
mirCoefHP <- readRDS(paste0("../", data.dir, "BLUP/mirHPBLUP_SNP.rds"))
# Congruent eQTL in additive and interactive models (SFig8)
Add_Int_overlap <- readRDS(paste0("../", data.dir, "scan_Overlaps/Add_Int_Overlap.rds"))
# Congruent mirQTL in additive and interactive models (SFig8)
mirAdd_Int_overlap <- readRDS(paste0("../", data.dir, "scan_Overlaps/mirAdd_Int_overlap_KJ.rds"))
# Source the CT_plotter function (SFig8)
source(paste0("../", script.dir, "CT_plotter.R"))
source(paste0("../", script.dir, "CT_plotter_miRNA.R"))
Supplemental Figure 1: Distribution of LOD scores significant at p < 0.05 in the genome scan with diet as an additive covariate, as determined by the null distribution of max LODs from 1,000 permuted genome scans.
p <- ggplot(mRes$addCov_FP, aes( x = lod)) +
geom_histogram(colour = "black",
fill = "dodgerblue3",
binwidth = 1.5) +
labs(x = "LOD",
y = "n")
p
#ggsave(paste0("../", fig.dir, "SupFig1.png"), p, dpi = 600,
# h = 6,
# w = 9)
Supplemental Figure 2: Magnitude of correlation is greater in non-mapping miRNA-mRNA pairs than in those with mRNA with eQTL. Expression data was used to generate correlations between pairs of miRNA (n = 246) and mRNA (n = 24,004). The results were filtered for BH-corrected significance (α ≤ 0.05) and Spearman’s rho was plotted against the eQTL status of the mRNA from the additive genome scan. Wilcoxon rank sum test found a significant difference between the absolute correlations of cis and trans-mapping eQTL (p < 2.22 x 10-16) with the median of trans-eQTL being 0.234, whereas cis was 0.233. A less significant trend was observed that non-mapping mRNA had lower miRNA correlations than mapping mRNA (p < 1.1 x 10-5). The median value for no eQTL was 0.232.
# Check structure of status_Add
str(res$status_Add)
# Factor this as desired for plotting and stats in mind
res$status_Add <- factor(res$status_Add, levels = c("cis","trans", "No eQTL"))
#Visualize
sfig2_Add <- ggplot(res, 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)+
theme_bw()+
labs(y = "Spearman's Rho (Absolute Value)",x = "eQTL Status")+
scale_x_discrete(limits = c("cis", "trans", "No eQTL"),
labels=c("cis" = paste0("Cis\n(n=",table(res$status_Add)[[1]],")"),
"trans" = paste0("Trans\n(n=",table(res$status_Add)[[3]],")"),
"No eQTL" = paste0("No eQTL\n(n=",table(res$status_Add)[[2]],")")))+
theme(axis.text = element_text(size = 16),
axis.title.x = element_text(margin = margin(t = 20),
size = 16),
axis.title.y = element_text(margin = margin(r = 20),
size = 16),
legend.position="none")
sfig2_Add
#ggsave(paste0("../", fig.dir, "SupFig2.png"), sfig2_Add, dpi = 600,
# h = 9,
# w = 6)
Supplemental Figure 3: eQTL and mirQTL have similar effect size in Diversity Outbred mice fed different diets (A, B) Wilcoxon rank sum tests reveal a significant difference in the phenotypic variance explained between cis and trans-eQTL of mice fed a HFCA diet or a HP diet. The median phenotypic variance explained by cis-eQTL in the HFCA is 0.349 whereas it is 0.244 in trans. Similarly, in the HP, the variance explained by cis-eQTL is 0.356 and 0.282 in trans. (C, D) Wilcoxon rank sum tests reveal a significant difference in the phenotypic variance explained between cis and trans-mirQTL and a similar but insignificant (p ≥ 0.05) difference in the phenotypic variance explained between the cis and trans-mirQTL of mice fed a HP diet. The median phenotypic variance explained by cis-mirQTL in the HFCA is 0.465 whereas it is 0.251 in trans. Similarly, in the HP, the variance explained by cis-mirQTL is 0.525 and 0.302 in trans.
# HFCA; eQTL
# ggplot + geom_boxplot
mRNA_phenVarBox_HFCA <- ggplot(mRes$HFCA_FP,
aes(x = CT_strict,
y = qtlEffect,
fill = CT_strict)) +
geom_boxplot(alpha = 0.8)
# Add theme for manuscript
a1_HFCA <- mRNA_phenVarBox_HFCA +
theme_bw() +
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))) +
ylim(0,1)+
scale_fill_manual(values = c("dodgerblue","firebrick2"))+
labs(y ="Phenotypic variance explained ", x = "eQTL Status")+
scale_x_discrete(labels=c("cis" =
paste0("Cis\n(n=",
table(mRes$HFCA_FP$CT_strict)[[1]],")"),
"trans" =
paste0("Trans\n(n=",table(mRes$HFCA_FP$CT_strict)[[2]],")")))+
stat_compare_means(comparisons = list(comp1 = c(1,2)),
method = "wilcox",
label = "p.format",
paired = FALSE)+
guides(fill = FALSE)
# HFCA; mirQTL
# ggplot + geom_boxplot
mir_phenVarBox_HFCA <- ggplot(mRes$mirHFCA, aes(x = CT, y = qtlEffect, fill = CT)) +
geom_boxplot(alpha = 0.8)
# Add theme for manuscript
b1_HFCA <- mir_phenVarBox_HFCA +
theme_bw() +
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))) +
ylim(0,1)+
scale_fill_manual(values = c("dodgerblue","firebrick2"))+
labs(y ="Phenotypic variance explained ", x = "mirQTL Status")+
scale_x_discrete(labels=c("cis" = paste0("Cis\n(n=",table(mRes$mirHFCA$CT)[[1]],")"),
"trans" = paste0("Trans\n(n=",table(mRes$mirHFCA$CT)[[2]],")")))+
stat_compare_means(comparisons = list(comp1 = c(1,2)),
method = "wilcox",
label = "p.format",
paired = FALSE)+
guides(fill = FALSE)
# HP; eQTL
# ggplot + geom_boxplot
mRNA_phenVarBox_HP <- ggplot(mRes$HP_FP, aes(x = CT_strict, y = qtlEffect, fill = CT_strict)) +
geom_boxplot(alpha = 0.8)
# Manuscript theme
a2_HP <- mRNA_phenVarBox_HP +
theme_bw() +
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))) +
ylim(0,1)+
scale_fill_manual(values = c("dodgerblue","firebrick2"))+
labs(y ="Phenotypic variance explained ", x = "eQTL Status")+
scale_x_discrete(labels=c("cis" = paste0("Cis\n(n=",table(mRes$HP_FP$CT_strict)[[1]],")"),
"trans" = paste0("Trans\n(n=",table(mRes$HP_FP$CT_strict)[[2]],")")))+
stat_compare_means(comparisons = list(comp1 = c(1,2)),
method = "wilcox",
label = "p.format",
paired = FALSE)+
guides(fill = FALSE)
# HP; mirQTL
# ggplot + geom_boxplot
mir_phenVarBox_HP <- ggplot(mRes$mirHP, aes(x = CT, y = qtlEffect, fill = CT)) +
geom_boxplot(alpha = 0.8)
# Manuscript theme
b2_HP <- mir_phenVarBox_HP +
theme_bw() +
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))) +
ylim(0,1)+
scale_fill_manual(values = c("dodgerblue","firebrick2"))+
labs(y ="Phenotypic variance explained ", x = "mirQTL Status")+
scale_x_discrete(labels=c("cis" = paste0("Cis\n(n=",table(mRes$mirHP$CT)[[1]],")"),
"trans" = paste0("Trans\n(n=",table(mRes$mirHP$CT)[[2]],")")))+
stat_compare_means(comparisons = list(comp1 = c(1,2)),
method = "wilcox",
label = "p.format",
paired = FALSE)+
guides(fill = FALSE)
# Save!
out_PhenoVariance_diet <- ggarrange(a1_HFCA,a2_HP,b1_HFCA,b2_HP, ncol = 2, nrow = 2, widths = c(2,2))
out_PhenoVariance_diet
#ggsave(paste0("../", fig.dir, "SupFig5.png"),out_PhenoVariance_diet,
# h = 8,
# w = 14,
# dpi = 600)
Supplemental Figure 4: Alleles from Wild Derived strains contribute to eQTL and mir-eQTL. (A) Boxplot representation of center-scaled BLUP coefficients from all significant eQTL in the HFCA-fed mice. 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. (B) Boxplot representation of center-scaled BLUP coefficients from all significant mirQTL in the HFCA-fed mice. Kruskal-Wallis rank sum test indicates no significant differences (p = 0.34) among founder effects. C) Boxplot of center-scaled BLUP coefficients from all significant eQTL in the HP-fed mice. Kruskal-Wallis rank sum test indicates significant differences (p < 2.2 x 10-16) among founder effects. Dunn’s post-hoc test confirms significantly (p < 1.0 x 10-4, BH-adjusted) larger effects from the CAST/EiJ and PWK/PhJ strains relative to all other strains. (D) Boxplot of center-scaled BLUP coefficients from all significant mirQTL in the HP-fed mice. Kruskal-Wallis rank sum test does not indicate significant differences (p < 0.064) among founder effects.
founders <- c("A/J", "C57BL/6J", "129S1/SvImJ", "NOD/LtJ", "NZO/HlLtJ", "CAST/EiJ", "PWK/PhJ","WSB/EiJ")
inDat <- HFCABLUP_SNP # load this in by going to figure 3, hfcaBLUP_SNP.rds
names(inDat) <- founders # assign column names to founder mice
# Scale/center for plots
inDat <- t(scale(t(inDat))) # scale centers and/or scales the columns of a numeric matrix. Here we scale by marker.
# Prepare for plotting
mtest <- melt(inDat) # make data long
names(mtest)[2] <- "variable"
mtest$variable <- factor(mtest$variable, levels = founders)
p1 <- 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")
#mirQTL Plots-----------------------------------------------------
mirCoef <- t(scale(t(mirCoef)))
colnames(mirCoef) <- founders
mtest <- melt(mirCoef)
names(mtest)[2] <- "variable"
mtest$variable <- factor(mtest$variable, levels = founders)
p2 <- 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")
ggarrange(p1,p2, ncol = 2)
#ggsave(paste0("../", fig.dir, "SupFig6AB.png"),ggarrange(p1,p2, ncol = 2),
# w = 14,
# h = 5,
# dpi = 600)
#---------------------------
inDat <- HPBLUP_SNP
names(inDat) <- founders
inDat <- t(scale(t(inDat)))
mtest <- melt(inDat)
names(mtest)[2] <- "variable"
mtest$variable <- factor(mtest$variable, levels = founders)
p3 <- 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")
#mirQTL Plots-----------------------------------------------------
mirCoefHP <- t(scale(t(mirCoefHP)))
colnames(mirCoefHP) <- founders
mtest <- melt(mirCoefHP)
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")
ggarrange(p3,p4, ncol = 2)
## Warning: Removed 48 rows containing non-finite values (stat_boxplot).
#ggsave(paste0("../", fig.dir, "SupFig6CD.png"),ggarrange(p3,p4, ncol = 2),
# w = 14,
# h = 5,
# dpi = 600)
Supplemental Figure 5: Magnitude of correlation of miRNA-mRNA pairs with eQTL status differs by diet. Expression data was used to generate correlations between pairs of miRNA (n = 246) and mRNA (n = 24,004). The results were filtered for FDR-corrected significance (α ≤ 0.05) and Spearman’s rho was plotted against the eQTL status of the mRNA from the HFCA model (A) and the HP model (B). Wilcoxon rank sum test reveals a significant difference between the absolute correlations of cis and trans-mapping eQTL in the HFCA and HP model (HFCA, p < 2.2 x 10-16; HP, p < 0.0005).
# HFCA
Sfig7_HFCA <- ggplot(res, aes(x=status_HFCA, y=abs(rho), fill = status_HFCA))+
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)+
theme_bw()+
labs(y = "Spearman's Rho (Absolute Value)",x = "eQTL Status")+
scale_x_discrete(limits = c("cis", "trans", "No eQTL"),
labels=c("cis" = paste0("Cis\n(n=",table(res$status_HFCA)[[1]],")"),
"trans" = paste0("Trans\n(n=",table(res$status_HFCA)[[3]],")"),
"No eQTL" = paste0("No eQTL\n(n=",table(res$status_HFCA)[[2]],")")))+
theme(axis.text = element_text(size = 12),
axis.title.x = element_text(margin = margin(t = 20),
size = 12),
axis.title.y = element_text(margin = margin(r = 20),
size = 12),
legend.position="none")
# HP
Sfig7_HP <- ggplot(res, aes(x=status_HP, y=abs(rho), fill = status_HP))+
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)+
theme_bw()+
labs(y = "Spearman's Rho (Absolute Value)",x = "eQTL Status")+
scale_x_discrete(limits = c("cis", "trans", "No eQTL"),
labels=c("cis" = paste0("Cis\n(n=",table(res$status_HP)[[1]],")"),
"trans" = paste0("Trans\n(n=",table(res$status_HP)[[3]],")"),
"No eQTL" = paste0("No eQTL\n(n=",table(res$status_HP)[[2]],")")))+
theme(axis.text = element_text(size = 12),
axis.title.x = element_text(margin = margin(t = 20),
size = 12),
axis.title.y = element_text(margin = margin(r = 20),
size = 12),
legend.position="none")
# Save the diet figures together
out_correlationCisTransNone <- ggarrange(Sfig7_HFCA, Sfig7_HP, ncol = 2, nrow = 1, widths = c(2,2))
out_correlationCisTransNone
#ggsave(paste0("../", fig.dir, "SupFig7.png"),out_correlationCisTransNone,
# h = 8,
# w = 12,
# dpi = 600)
Supplemental Figure 6: Overlap between results of additive and interactive models revealed large architectural differences in trans-eQTLs compared to cis-eQTLs. Differences in the genetic architecture of eQTL and mirQTL in the additive and interactive QTL analyses. Global architecture of significant eQTL (A) and mirQTL (B) from the interactive scan. We observed larger overlaps between additive and interactive scans in cis-regulated QTL (C and D) and sparser overlaps between trans-regulated QTL (E and F).
# Square plots
#png(paste0("../",fig.dir,"SupFig8A.png"), width = 16, height = 12, units= "in" , res = 500)
CTplot(mRes$intCov_FP, ptSize = 1, fntSize = 20, xlab_shift = 10, ylab_shift = 10)
#dev.off()
#png(paste0("../",fig.dir,"SupFig8B.png"), width = 16, height = 12, units= "in" , res = 500)
CTplot_mir(mRes$mirIntCov, ptSize = 3, fntSize = 20, xlab_shift = 10, ylab_shift = 10)
#dev.off()
#----------------------------------------------------
# Venn Diagrams
# eQTL
# See distribution of classifications
table(Add_Int_overlap$CTStrict)
# Make simplified DF
int <- mRes$intCov_FP
add <- mRes$addCov_FP
# See dimensions
int_cis <- int[int$CT_strict == "cis",] #4415
add_cis <- add[add$CT_strict == "cis",] #5207
# Make new venn diagrams to reflect these numbers:
grid.newpage()
Add_Int_cis_plot <- draw.pairwise.venn(area1 = 5207, # additiv cis
area2 = 4415, # interactice cis
cross.area = 4321,
fill = c("red","blue"),
lty = "blank",
#category =c("Additive","Interactive"),
cex=2,
fontfamily = "serif")
# Save
#png(file = paste0("../", fig.dir, "SupFig8C.png"))
grid.draw(Add_Int_cis_plot)
#dev.off()
# Repeat for trans
add_trans <- add[add$CT_strict == "trans",] #1485
int_trans <- int[int$CT_strict == "trans",] #1740
# Plot
grid.newpage()
Add_Int_trans_plot <- draw.pairwise.venn(area1 = 1485, # additiv trans
area2 = 1740, # interactice trans
cross.area = 420,
fill = c("red","blue"),
lty = "blank",
#category =c("Additive","Interactive"),
cex=2,
fontfamily = "serif")
# Save
#png(file = paste0("../", fig.dir, "SupFig8E.png"))
grid.draw(Add_Int_trans_plot)
#dev.off()
#----------------------
# mirQTL
table(mirAdd_Int_overlap$CT)
# Simplified DF
mir_add <- mRes$mirAddCov
mir_int <- mRes$mirIntCov
# See dimensions
add_cis <- mir_add[mir_add$CT == "cis",] #18
int_cis <- mir_int[mir_int$CT == "cis",] #9
# Make new venn diagrams to reflect these numbers:
grid.newpage()
mirAdd_Int_cis_plot <- draw.pairwise.venn(area1 = 18, # additiv cis
area2 = 9, # interactice cis
cross.area = 8,
fill = c("red","blue"),
lty = "blank",
#category =c("Additive","Interactive"),
cex=2,
fontfamily = "serif")
# Save
#png(file = paste0("../", fig.dir, "SupFig8D.png"))
grid.draw(mirAdd_Int_cis_plot)
#dev.off()
# Repeat for trans
mir_add_trans <- mir_add[mir_add$CT == "trans",] # 19
mir_int_trans <- mir_int[mir_int$CT == "trans",] # 14
# Plot
grid.newpage()
mirAdd_Int_trans_plot <- draw.pairwise.venn(area1 = 19, # additiv trans
area2 = 14, # interactice trans
cross.area = 9,
fill = c("red","blue"),
lty = "blank",
#category =c("Additive","Interactive"),
cex=2,
fontfamily = "serif")
# Save
#png(file = paste0("../", fig.dir, "SupFig8F.png"))
grid.draw(mirAdd_Int_trans_plot)
#dev.off()