library(ggplot2)
library(cowplot)

Setup Data For Panel A

#import
tmt.acs <- readRDS("tmt.acs")# mutant data
tmt.6H <- readRDS("tmt.6H")# 6-hybrid data
#Gene of interest list
GOIs <-readRDS("GOIs")
#define genes as above mid-parent (AMP), below mid-parent (BMP), or mid-parent (MP)
AMP <- tmt.6H[tmt.6H$l2.A682xB73.MP > 0 & tmt.6H$ttest.A682xB73.MP < 0.05,] #AMP in A682xB73
BMP <- tmt.6H[tmt.6H$l2.A682xB73.MP < 0 & tmt.6H$ttest.A682xB73.MP < 0.05,] #BMP in A682xB73
tmt.acs$A682B <- "MP"
tmt.acs$A682B[tmt.acs$Gene %in% AMP$Gene] <- "AMP"
tmt.acs$A682B[tmt.acs$Gene %in% BMP$Gene] <- "BMP"
tmt.acs$A682B <- factor(tmt.acs$A682B, levels = c("MP","AMP","BMP"))
tmt.acs <- tmt.acs[order(tmt.acs$A682B),]
#calculate -log10(pvalue) for y axis of volcano plot
tmt.acs$l10.pvalue <- -log10(tmt.acs$ttest)
#merge mutant and hybrid data
tmt.6H.acs <- merge.data.frame(tmt.acs[,c("Accession", "l2.ACS.B73")], tmt.6H[,c("Accession", "l2.A682xB73.MP")], by.x = "Accession", by.y = "Accession", all=F)
#Convert accession to gene
tmt.6H.acs$Gene <- substr(tmt.6H.acs$Accession, start = 1, stop = 14)
#set plotting axis limits
tmt.acs$l2.ACS.B73[tmt.acs$l2.ACS.B73 < -2] <- -2
tmt.acs$l2.ACS.B73[tmt.acs$l2.ACS.B73 > 2] <- 2
tmt.acs$l10.pvalue[tmt.acs$l10.pvalue > 7.5] <- 7.5

Setup Data For Panel B

#merge acs with A682xB73 data
tmt.acs.a682 <- merge.data.frame(tmt.acs[,c("Accession", "l2.ACS.B73")], tmt.6H[,c("Accession", "l2.A682xB73.MP")], by.x = "Accession", by.y = "Accession", all=F)
#convert accession to gene
tmt.acs.a682$Gene <- substr(tmt.acs.a682$Accession, start = 1, stop = 14)
#Set plotting axis limits
tmt.acs.a682$l2.ACS.B73[tmt.acs.a682$l2.ACS.B73 <= -.75] <- -.75
tmt.acs.a682$l2.ACS.B73[tmt.acs.a682$l2.ACS.B73 >= .75] <- .75
tmt.acs.a682$l2.A682xB73.MP[tmt.acs.a682$l2.A682xB73.MP <= -.75] <- -.75
tmt.acs.a682$l2.A682xB73.MP[tmt.acs.a682$l2.A682xB73.MP >= .75] <- .75
#Set GOIs
tmt.acs.a682 <- tmt.acs.a682
tmt.acs.a682$Category <- "Other"
tmt.acs.a682[(tmt.acs.a682$Gene %in% GOIs$CytoRibo),"Category"] <- "Cytosolic ribosome"
tmt.acs.a682[(tmt.acs.a682$Gene %in% GOIs$PhANGs),"Category"] <- "PhANGs"
tmt.acs.a682[(tmt.acs.a682$Gene %in% GOIs$PhAPGs),"Category"] <- "PhAPGs"
tmt.acs.a682[(tmt.acs.a682$Gene %in% GOIs$NE.PtRibo),"Category"] <- "NE plastid ribosome"
tmt.acs.a682[(tmt.acs.a682$Gene %in% GOIs$PE.ribo),"Category"] <- "PE plastid ribosome"
tmt.acs.a682[(tmt.acs.a682$Gene %in% GOIs$EtBio),"Category"] <- "Ethylene biosynthesis"
#factor
tmt.acs.a682$Category <- factor(tmt.acs.a682$Category, levels = c("Other", "Cytosolic ribosome", "PhANGs", "PhAPGs", "NE plastid ribosome", "PE plastid ribosome", "Ethylene biosynthesis"))
tmt.acs.a682 <- tmt.acs.a682[(order(tmt.acs.a682$Category)),]

Define function to retrieve legend from plot

get_legend<-function(myggplot){
  tmp <- ggplot_gtable(ggplot_build(myggplot))
  leg <- which(sapply(tmp$grobs, function(x) x$name) == "guide-box")
  legend <- tmp$grobs[[leg]]
  return(legend)
}

Plot Panel A

#Volcano plot
tmt.acs.vol <- ggplot(tmt.acs, mapping = aes(x=l2.ACS.B73, y=l10.pvalue, color=A682B))+
  geom_point(size=1)+
  scale_color_manual(values = c("MP"="grey", "AMP"="red", "BMP"="blue"), name="")+
  theme_bw()+
  xlab("log2 acs/B73")+
  ylab("-log10 p-value")+
  ggtitle("A")+
  geom_vline(xintercept = 0)+
  geom_hline(yintercept = -log10(0.05))+
  theme(aspect.ratio=1, legend.position = "bottom", axis.title = element_text(size=10), axis.text=element_text(size=10), legend.text = element_text(size=10), plot.title = element_text(hjust = -.05, vjust = -2))+
  coord_cartesian(expand = F)
tmt.acs.vol.l <- get_legend(tmt.acs.vol)
tmt.acs.vol <- tmt.acs.vol+theme(legend.position = "none")

Plot Panel B

#Set colorblind-friendly palette for plotting
colors.p <- c("Other"="black", "Cytosolic ribosome"="#B50000", "PE plastid ribosome"="#6A00EA", "NE plastid ribosome"="#E5A3FF", "PhAPGs"="#008275", "PhANGs"="#81E401", "Ethylene biosynthesis"="#EF8F6D")
#Scatter plot with marginal density curves
    tmt.acs.a682.scat <- ggplot(tmt.acs.a682, mapping=aes(x=l2.A682xB73.MP, y=l2.ACS.B73, color=Category, shape=Category))+
    geom_point(aes(color=Category), size=1)+
  scale_color_manual(values=colors.p, name="Category")+
      scale_shape_manual(values=c(1,19,19,19,19,19,19), guide="none")+
    theme_classic()+
      xlim(c(-.75,.75))+
      ylim(c(-.75,.75))+
      geom_vline(xintercept=0, linetype="dashed", color = "gray60", size=0.5)+
      geom_hline(yintercept=0, linetype="dashed", color = "gray60", size=0.5)+
      ylab("log2 acs/B73")+
      xlab("log2 A682xB73/MP")+
      ggtitle("B")+
      theme(aspect.ratio=1)+
      coord_fixed(expand = F)+
      theme(axis.text=element_text(size=10), legend.text=element_text(size=10),axis.title = element_text(size = 10),
        legend.position="bottom", legend.title = element_blank(), plot.title = element_text( vjust = -2))+
      guides(col=guide_legend(ncol = 2))
    tmt.acs.a682.scat.l <- get_legend(tmt.acs.a682.scat)
    tmt.acs.a682.scat <- tmt.acs.a682.scat+theme(legend.position = "none")

  ##Add density curves
    ###set up marginal density plot X axis
    densityplotX<-axis_canvas(tmt.acs.a682.scat, axis="x")+
      geom_density(aes(tmt.acs.a682[,"l2.A682xB73.MP"], color=tmt.acs.a682[,"Category"]), size=0.5)+
    scale_color_manual(values=colors.p, name="Category")+
      geom_vline(xintercept=0, linetype="dashed", color = "gray60", size=0.5)+
      theme_void()+
      theme(legend.title = element_blank())
    ###set up marginal density plot Y axis
    densityplotY<-axis_canvas(tmt.acs.a682.scat, axis="y", coord_flip = T)+
      geom_density(aes(tmt.acs.a682[,"l2.ACS.B73"], color=tmt.acs.a682[,"Category"]), size=0.5)+
      scale_color_manual(values=colors.p, name="Category")+
      scale_linetype_manual(values=c("dotted", "solid"))+
      geom_vline(xintercept=0, linetype="dashed", color = "gray60", size=0.5)+
      theme_void()+
      theme(legend.title = element_blank())+
      coord_flip()
    ###put two plots together
      tmt.acs.a682.scat <- insert_xaxis_grob(tmt.acs.a682.scat, densityplotX, grid::unit(0.275, "null"), position = "top")
      tmt.acs.a682.scat <- insert_yaxis_grob(tmt.acs.a682.scat, densityplotY, grid::unit(0.275, "null"), position = "right")
#Define layout for panels
  des <- "AB
          CD"
patchwork::wrap_plots(A=tmt.acs.vol,B=tmt.acs.a682.scat,C=tmt.acs.vol.l,  D=tmt.acs.a682.scat.l,nrow = 2, heights = c(3,1), design = des)