220624 figures

article/data1/files/QTL_plot.R

+library(ggplot2)
+library(ggrepel)
+library(tidyr)
+library(dplyr)
+library(grid)
+
+
+#Modified from https://shiring.github.io/ggplot2/2017/02/12/qtl_plots
+
+qtl_plot <- function(input,              # data frame input from scanone
+                     mult.pheno = FALSE, # multiple phenotypes?
+                     model = "normal",   # model used in scanone
+                     chrs = NA,          # chromosomes to display
+                     lod = NA,           # LOD threshold
+                     rug = FALSE,        # plot marker positions as rug?
+                     strip.pos=NA,       # position of facet labels
+                     ncol = NA,          # number of columns for facetting
+                     labels = NA,         # optional dataframe to plot QTL labels
+                     ylab = "LOD",       #optional name for y axis
+                     title = "",         #optional name for the graph
+                     x.label = element_text(),   #optional: use element_blank to remove labels
+                     size = 1,
+                     int = NA            #optional confidence interval df: 3 columns chr, pos and lod with 3 lines: first limit, peak and second limit
+) {
+  
+  # if we have multiple phenotypes and/or a 2part model, gather input
+  if (mult.pheno & model == "2part") {
+    input <- gather(input, group, lod, grep("pheno", colnames(input)))
+  } else if (mult.pheno) {
+    input <- gather(input, group, lod, grep("pheno", colnames(input)))
+  } else if (model == "2part") {
+    input <- gather(input, method, lod, lod.p.mu:lod.mu)
+  }
+  
+  # if not all chromosomes should be displayed, subset input
+  if (!is.na(chrs)[1]) {
+    input <- input[as.character(input$chr) %in% chrs, ]
+  }
+  
+  # if there is more than one LOD column, gather input
+  if (!any(colnames(input) == "lod")) {
+    input$lod <- input[, grep("lod", colnames(input))]
+  }
+  
+  # if no number of columns for facetting is defined, plot all in one row
+  if (is.na(ncol)) {
+    ncol <- length(unique(input$chr))
+  }
+  
+  # facet labels positions
+  if (is.na(strip.pos)) {
+    strip.pos <- "top"
+  }
+  
+  # if labels are set and there is no name column, set from rownames
+  if (!is.na(labels)[1]) {
+    if (is.null(labels$name)) {
+      labels$name <- rownames(labels)
+    }
+  }
+  
+  # plot input data frame position and LOD score
+  plot <- ggplot(input, aes(x = pos, y = lod)) + {
+    
+    # if LOD threshold is given, plot as horizontal line
+    if (!is.na(lod)[1] & length(lod) == 1) geom_hline(yintercept = lod, linetype = "solid")
+  } + {
+    
+    if (!is.na(lod)[1] & length(lod) > 1) geom_hline(data = lod, aes(yintercept = lod, linetype = group))
+  } + {
+    
+    # plot rug on bottom, if TRUE
+    if (rug) geom_rug(size = 0.1, sides = "b")
+  } + {
+    
+    # if input has column method but not group, plot line and color by method
+    if (!is.null(input$method) & is.null(input$group)) geom_line(aes(color = method), size = size, alpha = 0.6)
+  } + {
+    
+    # if input has column group but not method, plot line and color by group
+    if (!is.null(input$group) & is.null(input$method)) geom_line(aes(color = group), size = size, alpha = 0.6)
+  } + {
+    
+    # if input has columns method and group, plot line and color by method & linetype by group
+    if (!is.null(input$group) & !is.null(input$method)) geom_line(aes(color = method, linetype = group), size = size, alpha = 0.6)
+  } + {
+    
+    # set linetype, if input has columns method and group
+    if (!is.null(input$group) & !is.null(input$method)) scale_linetype_manual(values = c("solid", "twodash", "dotted"))
+  } + {
+    
+    # if input has neither columns method nor group, plot black line
+    if (is.null(input$group) & is.null(input$method)) geom_line(size = size, alpha = 0.6)
+  } + {
+    
+    # if QTL positions are given in labels df, plot as point...
+    if (!is.na(labels)[1]) geom_point(data = labels, aes(x = pos, y = lod))
+  } + {
+    
+    # ... and plot name as text with ggrepel to avoid overlapping
+    if (!is.na(labels)[1]) geom_text_repel(data = labels, aes(x = pos, y = lod, label = name), nudge_y = 0.5)
+  } + 
+    # facet by chromosome
+    facet_wrap(~ chr, ncol = ncol, scales = "free_x",strip.position=strip.pos) +
+    # minimal plotting theme
+    theme_minimal() +
+    # increase strip title size
+    theme(strip.text = element_text(face = "bold", size = 12),
+          plot.title = element_text(hjust = 0.5, size = 12, face = "bold"),
+          axis.text.x = x.label,
+          panel.spacing = unit(0, "lines"), #space between chromosomes graphs
+          panel.border = element_rect(colour = "#d9d9d9", fill=NA, size=0.4, linetype = "solid"),
+          panel.grid = element_blank()) + 
+    # use RcolorBrewer palette
+    scale_color_brewer(palette = "Set1") +
+    # Change plot labels
+    labs(x = "Chromosome",
+         y = ylab,
+         color = "",
+         linetype = "") +
+    ggtitle(title) 
+  
+  
+    if(is.data.frame(int)==TRUE){
+      if(!is.na(chrs)[1] & length(chrs)==1){
+        #add interval for peak
+        plot <- plot + geom_segment(x=int[1,2],y=min(input[,3]),xend=int[1,2],yend=int[1,3],linetype="dashed",color="darkblue") +
+        geom_segment(x=int[3,2],y=min(input[,3]),xend=int[3,2],yend=int[3,3],linetype="dashed",color="darkblue") +
+        geom_segment(x=int[2,2],y=min(input[,3]),xend=int[2,2],yend=int[2,3],color="darkblue")
+      } else {
+        stop("only one chromosome must be displayed to plot confidence interval")
+      }
+    } 
+  return(plot)
+}