Mercurial > repos > iuc > volcanoplot

comparison test-data/out.rscript @ 6:83c573f2e73c draft

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
"planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/volcanoplot commit e880a5736a65605ae8a38ca9bc6f17c92c798ccd"
author iuc
date Thu, 10 Jun 2021 08:38:12 +0000
parents 44608d0193ed
children ab01e379d29e
comparison
equal deleted inserted replaced
5:44608d0193ed 6:83c573f2e73c
5 options(show.error.messages = F, error = function() {cat(geterrmessage(), file = stderr()); q("no", 1, F)}) 5 options(show.error.messages = F, error = function() {cat(geterrmessage(), file = stderr()); q("no", 1, F)})
6 6
7 # we need that to not crash galaxy with an UTF8 error on German LC settings. 7 # we need that to not crash galaxy with an UTF8 error on German LC settings.
8 loc <- Sys.setlocale("LC_MESSAGES", "en_US.UTF-8") 8 loc <- Sys.setlocale("LC_MESSAGES", "en_US.UTF-8")
9 9
10 # Galaxy settings end -----------------------------------------------------
10 11
11 # Load packages ----------------------------------------------------------- 12 # Load packages -----------------------------------------------------------
12 13
13 suppressPackageStartupMessages({ 14 suppressPackageStartupMessages({
14 library(dplyr) 15 library(dplyr)
17 }) 18 })
18 19
19 20
20 # Import data ------------------------------------------------------------ 21 # Import data ------------------------------------------------------------
21 22
22 # Check if header is present by checking if P value column is numeric or not 23 results <- read.delim('/private/var/folders/zn/m_qvr9zd7tq0wdtsbq255f8xypj_zg/T/tmprh4qip75/files/d/2/2/dataset_d2255b46-f0f6-4900-8b9e-bd352e34f303.dat', header = TRUE)
23
24 first_line <- read.delim('/private/var/folders/zn/m_qvr9zd7tq0wdtsbq255f8xypj_zg/T/tmpfpemuuun/files/8/3/7/dataset_8374ef6b-02c7-46f1-afc9-408a2a6cbde4.dat', header = FALSE, nrow = 1)
25
26 first_pvalue <- first_line[, 3]
27
28 if (is.numeric(first_pvalue)) {
29 print("No header row detected")
30 results <- read.delim('/private/var/folders/zn/m_qvr9zd7tq0wdtsbq255f8xypj_zg/T/tmpfpemuuun/files/8/3/7/dataset_8374ef6b-02c7-46f1-afc9-408a2a6cbde4.dat', header = FALSE)
31 } else {
32 print("Header row detected")
33 results <- read.delim('/private/var/folders/zn/m_qvr9zd7tq0wdtsbq255f8xypj_zg/T/tmpfpemuuun/files/8/3/7/dataset_8374ef6b-02c7-46f1-afc9-408a2a6cbde4.dat', header = TRUE)
34 }
35 24
36 25
37 # Format data ------------------------------------------------------------ 26 # Format data ------------------------------------------------------------
38 27
39 # Create columns from the column numbers specified 28 # Create columns from the column numbers specified
50 # Set colours 39 # Set colours
51 colours <- setNames(c("cornflowerblue", "grey", "firebrick"), c(down, notsig, up)) 40 colours <- setNames(c("cornflowerblue", "grey", "firebrick"), c(down, notsig, up))
52 41
53 # Create significant (sig) column 42 # Create significant (sig) column
54 results <- mutate(results, sig = case_when( 43 results <- mutate(results, sig = case_when(
55 fdr < 0.05 & logfc > 0.0 ~ up, 44 fdr < 0.05 & logfc > 0.0 ~ up,
56 fdr < 0.05 & logfc < -0.0 ~ down, 45 fdr < 0.05 & logfc < -0.0 ~ down,
57 TRUE ~ notsig)) 46 TRUE ~ notsig))
58 47
59 48
60 # Specify genes to label -------------------------------------------------- 49 # Specify genes to label --------------------------------------------------
61 labelfile <- read.delim('/private/var/folders/zn/m_qvr9zd7tq0wdtsbq255f8xypj_zg/T/tmpfpemuuun/files/4/2/f/dataset_42fc8a63-f9cc-435b-9bb3-dd106b708cd9.dat') 50
51 # Import file with genes of interest
52 labelfile <- read.delim('/private/var/folders/zn/m_qvr9zd7tq0wdtsbq255f8xypj_zg/T/tmprh4qip75/files/5/e/5/dataset_5e5b8fb0-bf65-438e-9b5b-03a540d9aa5d.dat', header = TRUE)
53
54 # Label the genes of interest in results table
62 results <- mutate(results, labels = ifelse(labels %in% labelfile[, 1], labels, "")) 55 results <- mutate(results, labels = ifelse(labels %in% labelfile[, 1], labels, ""))
56
63 57
64 58
65 # Create plot ------------------------------------------------------------- 59 # Create plot -------------------------------------------------------------
66 60
67 pdf("out.pdf") 61 # Open file to save plot as PDF
68 p <- ggplot(results, aes(x = logfc, y = -log10(pvalue))) + 62 pdf("volcano_plot.pdf")
63
64 # Set up base plot
65 p <- ggplot(data = results, aes(x = logfc, y = -log10(pvalue))) +
69 geom_point(aes(colour = sig)) + 66 geom_point(aes(colour = sig)) +
70 scale_color_manual(values = colours) + 67 scale_color_manual(values = colours) +
71 scale_fill_manual(values = colours) +
72 theme(panel.grid.major = element_blank(), 68 theme(panel.grid.major = element_blank(),
73 panel.grid.minor = element_blank(), 69 panel.grid.minor = element_blank(),
74 panel.background = element_blank(), 70 panel.background = element_blank(),
75 axis.line = element_line(colour = "black"), 71 axis.line = element_line(colour = "black"),
76 legend.key = element_blank()) 72 legend.key = element_blank())
73
74 # Add gene labels
75 p <- p + geom_text_repel(data = filter(results, labels != ""), aes(label = labels),
76 min.segment.length = 0,
77 max.overlaps = Inf,
78 show.legend = FALSE)
77 79
78 80
79 81
80 82
81 83
82 84
83 # Set legend title 85 # Set legend title
84 p <- p + labs(colour = "") 86 p <- p + theme(legend.title = element_blank())
85 87
86 # Add gene labels in boxes 88 # Print plot
87 p <- p + geom_label_repel(aes(label = labels, fill = sig), 89 print(p)
88 segment.colour = "black",
89 colour = "white",
90 min.segment.length = 0,
91 show.legend = FALSE)
92 90
93 print(p) 91 # Close PDF graphics device
94 dev.off() 92 dev.off()
95
96
97 # Save RData -------------------------------------------------------------
98 save.image(file="volcanoplot.RData")
99 93
100 94
101 # R and Package versions ------------------------------------------------- 95 # R and Package versions -------------------------------------------------
102 sessionInfo() 96 sessionInfo()
103 97