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

Figure 7: Analysis of genome scans from HFCA-fed mice and HP-fed mice reveal environmentally driven differences in genetic architecture. (A) Global eQTL architecture for the HFCA model shows numerous, dense trans-bands, in contrast to the HP model with notably less co-localization of trans-eQTL (B). Both diet models (C and D) show similar precisions from eQTL with probesets on the same chromosome; resolution of these 3,987 eQTL from the HFCA-fed mice and 2,273 eQTL from HP-fed mice are estimated to be 0.29 Mb and 2.06 Mb with 93.81% and 92.32% of cis-eQTL occurring ±4 Mb from their probesets, respectively. This indicates that structural variants associated with local genetic regulation tend to occur close to the gene itself, regardless of environmental perturbations. (E and F) A large number of cis-eQTL overlap, despite environmental differences, while the lack of overlap in the trans-eQTL might be indicative of environmental effects.

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)
library(tidyr)
library(plyr)
library(scales)
library(OmicCircos)
library(VennDiagram)

getwd() # paths set such that wd is in markdown folder 
#setwd("/Users/kristenjames/Box/DO2_figshare/markdown/")

data.dir <-  "data/"
results.dir <-  "results/"
fig.dir <- "figures/"
script.dir <- "scripts/" 

# Source the CT_plotter function
source(paste0("../", script.dir, "CT_plotter.R"))
source(paste0("../", script.dir, "CT_plotter_miRNA.R"))
# Source the Precision_plotter function
source(paste0("../", script.dir, "Precision_plotter.R"))
source(paste0("../", script.dir, "Precision_plotter_miRNA.R"))

Version of R and packages used:

sessionInfo()

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. To generate Figure 9, this includes the congruent eQTL between the HFCA and HP-fed mice (matchDF).

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

# Congruent eQTL between HFCA and HP
matchDF <- readRDS(paste0("../", data.dir, "scan_Overlaps/HP_HFCAoverlap_KJ.rds"))

Plot Trans-band HFCA

A) Whole-genome plot with trans-bands mRNA

#png(paste0("../",fig.dir,"figure9A.png"), width = 16, height = 16,units= "in" , res = 500)
CTplot(mRes$HFCA_FP, ptSize = 1, fntSize = 20, xlab_shift = 10, ylab_shift = 10)

#dev.off()

Plot Trans-band HP

B) Whole-genome plot with trans-bands mRNA

#png(paste0("../",fig.dir,"figure9B.png"), width = 16, height = 16, units= "in" , res = 500)
CTplot(mRes$HP_FP, ptSize = 1, fntSize = 20, xlab_shift = 10, ylab_shift = 10)

#dev.off()

Plot Precision HFCA

C)

#png(paste0("../",fig.dir,"figure9C.png"), width = 16, height = 12, units= "in" , res = 500)
precPlot(mRes$HFCA_FP, run = "mRNA HFCA", fntSize = 12, xlab_shift = 5, ylab_shift = 5)
## Warning: Removed 1 rows containing missing values (geom_point).

#dev.off()

Plot Precision HP

D)

#png(paste0("../",fig.dir,"figure9D.png"), width = 16, height = 12, units= "in" , res = 500)
precPlot(mRes$HP_FP, run = "mRNA HP", fntSize = 12, xlab_shift = 5, ylab_shift = 5)

#dev.off()

Plot Venn Diagrams

E) Concordant HFCA and HP cis-eQTL and (F) trans-eQTL

table(matchDF$CTStrict)
# cis:cis = 1701
# trans:trans = 35
# cis:trans = 71
# trans:cis = 29
table(mRes$HFCA_FP$CT_strict) 
# cis = 3738
# trans = 1832
table(mRes$HP_FP$CT_strict) 
# cis = 2098
# trans = 1317

# HFCA x HP cis
grid.newpage()
HFCA_HP_cisPlot <-  draw.pairwise.venn(area1 = 3738, #HFCA
                   area2 = 2098,                  #HP
                   cross.area = 1701,             #overlap
                   fill = c("red","blue"),
                   lty = "blank",
                   #category =c("HFCA","HP"),
                   cex=2,
                   fontfamily = "serif")
# Save
#png(file = paste0("../",fig.dir,"figure9E.png"))
grid.draw(HFCA_HP_cisPlot)

#dev.off()

# HFCA HP trans
grid.newpage()
HFCA_HP_transPlot <-  draw.pairwise.venn(area1 = 1832, #HFCA
                      area2 = 1317,                  #HP
                      cross.area = 35,             #overlap
                      fill = c("red","blue"),
                      lty = "blank",
                      #category =c("HFCA","HP"),
                      cex=2,
                      fontfamily = "serif")
# Save
#png(file = paste0("../",fig.dir,"figure9F.png"))
grid.draw(HFCA_HP_transPlot)

#dev.off()