Code to generate Figure 9 of “Genetic architecture modulates diet induced hepatic mRNA and miRNA expression profiles”
Figure 9: Significant genotype by diet interactions on the mRNA expression. (A) The top 20 genes with the greatest percent increase or decrease of the Best Linear Unbiased Predictors (BLUP) founder allele coefficient for eQTLs having significant allele-diet interaction in HFCA diet with respect to HP diet. (B) BLUPs coefficient plot of eight founder mice strains to the Zfp982 QTL in HFCA diet and (C) HP diet. Color represents the eight founder mice strain as indicated.
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)
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] scales_1.1.0 plyr_1.8.6 tidyr_1.0.2 epiR_1.0-14
## [5] survival_3.1-12 dplyr_0.8.5 devtools_2.3.0 usethis_1.6.0
## [9] dunn.test_1.3.5 reshape_0.8.8 ggpubr_0.2.5 magrittr_1.5
## [13] 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 stringr_1.4.0 munsell_0.5.0
## [25] ggsignif_0.6.0 gtable_0.3.0 memoise_1.1.0 evaluate_0.14
## [29] knitr_1.28 callr_3.4.3 ps_1.3.2 parallel_3.5.1
## [33] fansi_0.4.1 Rcpp_1.0.4.6 backports_1.1.6 desc_1.2.0
## [37] pkgload_1.0.2 fs_1.4.1 bit_1.1-15.2 digest_0.6.25
## [41] stringi_1.4.6 processx_3.4.2 rprojroot_1.3-2 grid_3.5.1
## [45] cli_2.0.2 tools_3.5.1 tibble_3.0.1 RSQLite_2.2.0
## [49] crayon_1.3.4 pkgconfig_2.0.3 Matrix_1.2-18 ellipsis_0.3.0
## [53] data.table_1.12.8 prettyunits_1.1.1 assertthat_0.2.1 rmarkdown_2.1
## [57] 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 9A, this includes the coefficients of the founder effects for congruent eQTL between the HFCA-fed and HP-fed mice (matchDF
), the results of the allele-diet ANOVA (P_SNP_Diet
), and the annotation file (annot
). To generate Figure 9B and C, this includes Best Linear Unbiased Predictions (BLUP
) of the top 20 genes whose founder effects are most effected by diet (pList_HFCA
, pList_HP
), and the physical genetic map (pmap
).
# Master file
mRes <- readRDS(paste0("../", data.dir, "cluster_Inputs/masterRes_3_12_20.rds"))
# Founder effects of congruent eQTL between the HFCA and HP models
matchDF <- readRDS(paste0("../", data.dir, "allele_diet_Interaction/scanCoef_cors_HP_HFCA.rds"))
# All significant terms
P_SNP_Diet <- readRDS(paste0("../", data.dir, "allele_diet_Interaction/dietGenoIntP_KJ_V3.rds"))
# Annotation file
annot <- readRDS(paste0("../", data.dir, "cluster_Inputs/DO2annot_final.rds"))
# BLUP top 20 results
pList_HFCA <- readRDS(paste0("../", data.dir, "BLUP/scan1blupTop20_HFCA.rds"))
pList_HP <- readRDS(paste0("../", data.dir, "BLUP/scan1blupTop20_HP.rds"))
# Physical map
pmap <- readRDS(paste0("../", data.dir, "cluster_Inputs/DO2_pmap.rds"))
Calculate the percent difference of the founder effects at congruent eQTL of HFCA-fed mice using HP-fed mice as the reference.
# Take what you want; transcript column and then HP A-H and HFCA A-H (A-H represent the founders)
matchDF_A <- matchDF %>% select(A_QTL, AA, AB, AC, AD, AE, AF, AG, AH) # HFCA
matchDF_B <- matchDF %>% select(B_QTL, BA, BB, BC, BD, BE, BF, BG, BH) # HP
# Percent Difference
PercentDifference <- ((matchDF_A[,2:9] - matchDF_B[,2:9])/matchDF_B[,2:9]) * 100
# Assign row names; need to give a unique identifier since there are repeats
Num <- seq(1,nrow(PercentDifference),by=1)
rownames(PercentDifference)<- paste0(matchDF_A$A_QTL, "_", Num)
# Get PercentDifference dataframe ready for plotting; need to merge with annotations
PercentDifference$Transcript_ID <- row.names(PercentDifference)
# Scrub the "_count"
PercentDifference$Transcript_ID <- gsub("_.*","",PercentDifference$Transcript_ID)
# Grab transcript cluster ID and symbol column from annotation object
annot2 <- annot[,c(1,8,9,10)]
# Merge
PercentDifference <- merge(PercentDifference, annot2, by.x = "Transcript_ID", by.y = "transcript_cluster_id", all.x = T, all.y = F)
# Find the overlap with the eQTL with a significant allele-diet interaction term from the ANOVA
# Filter for BH-adjusted p<0.05 results
IntDF <- P_SNP_Diet[(which(P_SNP_Diet$P_adjust <0.05)),] # adjusted
dim(IntDF) #390
# Merge sig allele-diet results with annotation
P_SNP_Diet_annot <- left_join(IntDF, annot2, by = c("pheno" ="transcript_cluster_id"))
#Filter for significant results!
PercentDifference_Sig <- PercentDifference[PercentDifference$SYMBOL %in% P_SNP_Diet_annot$SYMBOL,] #67
A) Identify the top 20 genes whose founder effects are the most different when compared by diet.
# Find the top 20 most different genes
temp <- apply(abs(PercentDifference_Sig[2:9]),1, sum) # sum the absolute value of the % differences
temp <- as.data.frame(temp)
temp$TCIDs <- PercentDifference_Sig[,1] # make a transcript ID column
temp <- temp[order(temp$temp, decreasing = T),] # order from greatest to least
temp20 <- temp[1:20,] # select the top 20
# Filter
PercentDifference2 <- PercentDifference_Sig[PercentDifference_Sig$Transcript_ID %in% temp20$TCIDs,]
# Plotting DataFrame
gData <- gather(PercentDifference2, key = "Founder_allele", value = "Percent_difference", -Transcript_ID, -SYMBOL, -ENTREZID, -GENENAME)
# Check the structure of the variables
str(gData)
gData$Percent_difference <- as.numeric(gData$Percent_difference)
gData$Founder_allele <- as.factor(gData$Founder_allele)
gData$Founder_allele <- revalue(gData$Founder_allele, c("AA"="A", "AB"="B", "AC"="C", "AD"="D", "AE"="E", "AF"="F", "AG"="G", "AH"="H"))
# Make graph
p = ggplot(gData, aes(x=Founder_allele, y=SYMBOL)) +
geom_tile(aes(fill=Percent_difference), colour = "black") +
xlab("") + ylab("") +
theme(text=element_text(size=14)) +
scale_fill_gradient2(low="#2C7BB6",
mid="white",
high="#D7191C") +
labs(fill = "% Difference")
p
# Save!
#ggsave(paste0("../",fig.dir,"figure11A.png"), p,
# width = 6.5,
# height = 5,
# dpi = 300)
B) BLUP plot for HFCA-fed mice and (C) HP-fed mice
# Generate blup plots for manuscript. We will plot Zfp982
HFCA_Zfp_plot <- plot_coefCC(pList_HFCA$Zfp982, map = pmap, bgcol="white", legend = "topleft", main = "Zfp982 in HFCA-fed Mice")
HP_Zfp_plot <- plot_coefCC(pList_HP$Zfp982, map = pmap, bgcol="white", legend = "topleft", main = "Zfp982 in HP-fed Mice")
# Save
#pdf(paste0("../",fig.dir,"figure11B.pdf"), width = 6 , height = 5)
#plot_coefCC(pList_HFCA$Zfp982, map = pmap, bgcol="white", legend = "topleft", main = "Zfp982 in HFCA-fed Mice")
#dev.off()
#pdf(paste0("../",fig.dir,"figure11C.pdf"), width = 6 , height = 5)
#plot_coefCC(pList_HP$Zfp982, map = pmap, bgcol="white", legend = "topleft", main = "Zfp982 in HP-fed Mice")
#dev.off()