Merge pull request #80 from UMCUGenetics/feature/TICplots

release/DIMSv2.6.0

pipeline/scripts/2-DIMS.R

@@ -69,6 +69,11 @@ negRes=NULL
 
 x <- suppressMessages(xcmsRaw(filepath))
 
+# New: generate TIC plots. Prepare txt files with data for plots
+tic_intensity_persample <- cbind(round(x@scantime, 2), x@tic)
+colnames(tic_intensity_persample) <- c("retentionTime", "TIC")
+write.table(tic_intensity_persample, file=paste0(outdir, "/2-pklist/", sampname, "_TIC.txt"))
+
 load(paste(outdir, "breaks.fwhm.RData", sep="/"))
 
 # Create empty placeholders for later use

pipeline/scripts/3-averageTechReplicates.R

@@ -1,7 +1,9 @@
 #!/usr/bin/Rscript
 
-# load required packages 
-# none 
+# load required packages  for TIC plots
+library("ggplot2")
+library("gridExtra")
+pdf(NULL) 
 
 # define parameters 
 cmd_args <- commandArgs(trailingOnly = TRUE)
@@ -122,3 +124,69 @@ repl.pattern.filtered <- retVal$pattern
 save(repl.pattern.filtered, file=paste(outdir, "repl.pattern.negative.RData", sep="/"))
 write.table(remove_neg, file=paste(outdir, "miss_infusions_neg.txt", sep="/"), row.names=FALSE, col.names=FALSE ,sep= "\t")
 
+# New: generate TIC plots
+run_name <- basename(outdir)
+tic_input_dir <- paste(outdir, "2-pklist", sep = "/")
+
+# get replication pattern
+load(paste(outdir, "logs", "init.RData", sep="/"))
+# get misinjections
+bad_pos <- read.table(paste(outdir, "miss_infusions_pos.txt", sep="/"))
+bad_neg <- read.table(paste(outdir, "miss_infusions_neg.txt", sep="/"))
+
+# get all txt files
+tic_files = list.files(tic_input_dir, full.names=TRUE, pattern="*TIC.txt")
+all_samps <- sub('_TIC\\..*$', '', basename(tic_files))
+
+print("\n")
+print(paste0("reading TIC files", tic_files[1], " - till - ", tic_files[length(tic_files)]))
+# determine maximum intensity
+highest_tic_max <- 0
+for (file in tic_files) {
+  tic <- read.table(file)
+  this_tic_max <- max(tic$TIC)
+  if (this_tic_max > highest_tic_max) {
+    highest_tic_max <- this_tic_max
+    max_sample <- sub('_TIC\\..*$', '', basename(file))
+  }
+}
+
+tic_plot_list <- list()
+k = 0
+for (i in c(1:length(repl.pattern))) { 
+  tech_reps <- as.vector(unlist(repl.pattern[i]))
+  sampleName <- names(repl.pattern)[i]
+  for (j in 1:length(tech_reps)) {
+    k = k + 1
+    repl1.nr <- read.table(paste(paste(outdir, "2-pklist/", sep="/"), tech_reps[j], "_TIC.txt", sep=""))
+    bad_color_pos <- tech_reps[j] %in% bad_pos[[1]]
+    bad_color_neg <- tech_reps[j] %in% bad_neg[[1]]
+    if (bad_color_neg & bad_color_pos) {
+	    plotcolor = '#F8766D'
+    } else if (bad_color_pos) {
+	    plotcolor = "#ED8141"
+    } else if (bad_color_neg) {
+	    plotcolor = "#BF80FF"
+    } else {
+	    plotcolor = 'white'
+    }
+    tic_plot <- ggplot(repl1.nr, aes(retentionTime, TIC)) +
+      geom_line(linewidth = 0.3) +
+      geom_hline(yintercept = highest_tic_max, col = "grey", linetype = 2, linewidth = 0.3) +
+      labs(x = 't (s)', y = 'TIC', title = paste0(tech_reps[j], "  ||  ", sampleName)) +
+      theme(plot.background = element_rect(fill = plotcolor), axis.text = element_text(size = 4), axis.title = element_text(size = 4), plot.title = element_text(size = 6))
+    tic_plot_list[[k]] <- tic_plot
+  }
+
+}
+# Create a layout matrix with the size of the number of replicates as number of columns
+layout <- matrix(1:(10 * nrepl), 10, nrepl, TRUE)
+
+tic_plot_pdf <- marrangeGrob(grobs = tic_plot_list, 
+			     nrow = 10, 
+			     ncol = nrepl, 
+			     layout_matrix = layout, 
+			     top = quote(paste("TICs of run", run_name," \n colors: red = both modes misinjection, orange = pos mode misinjection, purple = neg mode misinjection \n ", g, "/", npages)))
+ggsave(filename = paste0(outdir, "/", run_name, "_TICplots.pdf"), tic_plot_pdf, width = 21, height = 29.7, units = "cm")
+
+

pipeline/scripts/AddOnFunctions/create_violin_plots.R

@@ -42,17 +42,25 @@ create_violin_plots <- function(pdf_dir, pt_name, metab_perpage, top_metab_pt=NU
   for (page_index in 1:length(metab_perpage)) {
     # extract list of metabolites to plot on a page
     metab_list_2plot <- metab_perpage[[page_index]]
+    # extract original data for patient of interest (pt_name) before cut-offs
+    pt_list_2plot_orig <- metab_list_2plot[which(metab_list_2plot$variable == pt_name), ]
     # cut off Z-scores higher than 20 or lower than -5 (for nicer plots)
     metab_list_2plot$value[metab_list_2plot$value >  20] <-  20
     metab_list_2plot$value[metab_list_2plot$value <  -5] <-  -5
     # extract data for patient of interest (pt_name)
     pt_list_2plot <- metab_list_2plot[which(metab_list_2plot$variable == pt_name), ]
+    # restore original Z-score before cut-off, for showing Z-scores in PDF
+    pt_list_2plot$value_orig <- pt_list_2plot_orig$value
     # remove patient of interest (pt_name) from list; violins will be made up of controls and other patients
     metab_list_2plot <- metab_list_2plot[-which(metab_list_2plot$variable == pt_name), ]
     # subtitle per page
     sub_perpage <- gsub("_", " ", page_headers[page_index])
     # for IEM plots, put subtitle on two lines
     sub_perpage <- gsub("probability", "\nprobability", sub_perpage)
+    # add size parameter for showing Z-score of patient per metabolite
+    Z_size <- rep(3, nrow(pt_list_2plot))
+    # set size to 0 if row is empty
+    Z_size[is.na(pt_list_2plot$value)] <- 0
 
     # draw violin plot. shape=22 gives square for patient of interest
     ggplot_object <- ggplot(metab_list_2plot, aes(x=value, y=HMDB_name)) +
@@ -61,6 +69,8 @@ create_violin_plots <- function(pdf_dir, pt_name, metab_perpage, top_metab_pt=NU
       geom_violin(scale="width") +
       geom_point(data = pt_list_2plot, aes(color=value), size = 3.5*circlesize, shape=22, fill="white") +
       scale_fill_gradientn(colors = colors_4plot, values = NULL, space = "Lab", na.value = "grey50", guide = "colourbar", aesthetics = "colour") +
+      # add Z-score value for patient of interest at x=16
+      geom_text(data = pt_list_2plot, aes(16, label = paste0("Z=", round(value_orig, 2))), hjust = "left", vjust = +0.2, size = Z_size) +
       # add labels. Use font Courier to get all the plots in the same location.
       labs(x = "Z-scores", y = "Metabolites", subtitle = sub_perpage, color = "z-score") + 
       theme(axis.text.y = element_text(family = "Courier", size=6)) +