library(ggplot2)
library(cowplot)

Setup data

#Gene of interest list
GOIs <-readRDS("GOIs")
#import data
cpm.SL <- readRDS("cpm.SL")
tmt.SL <- readRDS("tmt.SL")
#Label plots
tmt.SL$plot <- "A"
cpm.SL$plot <- "B"
#trim
tmt.SL <- tmt.SL[,c("plot", "Gene", "l2.B73xMo17.MP", "l10.ttest.MP")]
cpm.SL <- cpm.SL[,c("plot", "Gene", "l2.B73xMo17.MP", "l10.ttest.MP")]
#combine
tmt.cpm.SL <- rbind(tmt.SL,cpm.SL)

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)
}

Set colorblind-friendly palette for plotting

colors.p <- c("Other"="#D8D8D8", "Plastid-localized"="black", "Cytosolic ribosome"="#B50000", "PE plastid ribosome"="#6A00EA", "NE plastid ribosome"="#E5A3FF", "PhAPGs"="#008275", "PhANGs"="#81E401")

Plot!

plot<-unique(tmt.cpm.SL[,"plot"])
plotList<-vector(mode = "list", length = length(plot))

for(i in 1:length(plot)){
  df<-tmt.cpm.SL[which(tmt.cpm.SL[,"plot"]==plot[i]),]
  #Set plotting limits for x-axis
  df$l2.B73xMo17.MP[df$l2.B73xMo17.MP <= -1] <- -1
  df$l2.B73xMo17.MP[df$l2.B73xMo17.MP >= 1] <- 1
  #Set GOIs
  df <- df
  df$Category <- "Other"
  df[(df$Gene %in% GOIs$plastid),"Category"] <- "Plastid-localized"
  df[(df$Gene %in% GOIs$PhANGs),"Category"] <- "PhANGs"
  df[(df$Gene %in% GOIs$PhAPGs),"Category"] <- "PhAPGs"
  df[(df$Gene %in% GOIs$NE.PtRibo),"Category"] <- "NE plastid ribosome"
  df[(df$Gene %in% GOIs$PE.ribo),"Category"] <- "PE plastid ribosome"

  df$Category <- factor(df$Category, levels = c("Other", "Plastid-localized", "PhANGs", "PhAPGs", "NE plastid ribosome", "PE plastid ribosome"))
  df <- df[(order(df$Category)),]
  
  #scatter plot with marginal density curves
  df.scat <- ggplot(df, mapping=aes(x=l2.B73xMo17.MP, y=l10.ttest.MP, color=Category, shape=Category))+
    geom_point(aes(color=Category), size=1.5)+
    scale_color_manual(values=colors.p[c("Other", "Plastid-localized", "PhANGs", "PhAPGs", "NE plastid ribosome", "PE plastid ribosome")], name="Category")+
    scale_shape_manual(values=c(1,19,19,19,19,19), guide="none")+
    theme_classic()+
    xlim(c(-1,1))+
    geom_vline(xintercept=0, linetype="dashed", color = "gray40", size=0.5)+
    geom_hline(yintercept= 1.30103, linetype="dashed", color = "gray40", size=0.5)+
    ylab("-log10(p-value)")+
    ggtitle(plot[i])+
    theme(aspect.ratio=1)+
    coord_fixed(expand = F)+
    theme(plot.margin=grid::unit(c(-3,-3,0,-3), "mm"))+
    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(size=10, vjust = -3, hjust = -0.2), legend.spacing.x = unit(0.1, 'cm'), legend.spacing.y = unit(-3, 'cm'))+
    guides(col=guide_legend(ncol = 2))+
    xlab("log2(BxM/MP)")
  
  legend <- get_legend(df.scat)
  df.scat <- df.scat+theme(legend.position = "none")
  ##Add density curves
  ###set up marginal density plot X axis
  densityplotX<-axis_canvas(df.scat, axis="x")+
    geom_density(aes(df[,"l2.B73xMo17.MP"], color=df[,"Category"]), size=0.5)+
    scale_color_manual(values=colors.p[c("Other", "Plastid-localized", "PhANGs", "PhAPGs", "NE plastid ribosome", "PE plastid ribosome")], name="Category")+
    theme_void()+
    theme(legend.position = "none")+
    theme(plot.margin=grid::unit(c(-3,-3,0,-3), "mm"))
  ###set up marginal density plot Y axis
  densityplotY<-axis_canvas(df.scat, axis="y", coord_flip = T)+
    geom_density(aes(df[,"l10.ttest.MP"], color=df[,"Category"]), size=0.5)+
    scale_color_manual(values=colors.p[c("Other", "Plastid-localized", "PhANGs", "PhAPGs", "NE plastid ribosome", "PE plastid ribosome")], name="Category")+
    scale_linetype_manual(values=c("dotted", "solid"))+
    theme_void()+
    theme(legend.position = "none")+
    coord_flip(expand = F)+
    theme(plot.margin=grid::unit(c(-3,-3,0,-3), "mm"))
  #put two plots together
  q <- insert_xaxis_grob(df.scat, densityplotX, grid::unit(0.275, "null"), position = "top")
  q <- insert_yaxis_grob(q, densityplotY, grid::unit((0.275), "null"), position = "right")
  if(i==1){
    plotList<-list()
    plotList[[1]]<-q
  }
  else{
    plotList <- c(plotList, list(q))
  }
}
plotList <- c(plotList, list(legend))
des <- "AB
        CC"
r<-patchwork::wrap_plots(A=plotList[[1]], B=plotList[[2]], C=plotList[[3]], nrow = 2, heights = c(3,1.2), design = des)
ggsave(plot=r, filename ="ExpressionHeterosis.jpg", width = 4.5, height = 3)
r