annotate R/mutationSpectra_Galaxy.r @ 7:eda59b985b1c draft default tip

Uploaded
author iarc
date Mon, 13 Mar 2017 08:21:19 -0400
parents 46a10309dfe2
children
Ignore whitespace changes - Everywhere: Within whitespace: At end of lines:
rev   line source
7
eda59b985b1c Uploaded
iarc
parents: 6
diff changeset
1 #!/usr/bin/Rscript
eda59b985b1c Uploaded
iarc
parents: 6
diff changeset
2
eda59b985b1c Uploaded
iarc
parents: 6
diff changeset
3 #-----------------------------------#
eda59b985b1c Uploaded
iarc
parents: 6
diff changeset
4 # Author: Maude #
eda59b985b1c Uploaded
iarc
parents: 6
diff changeset
5 # Script: mutationSpectra_Galaxy.r #
eda59b985b1c Uploaded
iarc
parents: 6
diff changeset
6 # Last update: 23/07/15 #
eda59b985b1c Uploaded
iarc
parents: 6
diff changeset
7 #-----------------------------------#
eda59b985b1c Uploaded
iarc
parents: 6
diff changeset
8
eda59b985b1c Uploaded
iarc
parents: 6
diff changeset
9 #########################################################################################################################################
eda59b985b1c Uploaded
iarc
parents: 6
diff changeset
10 # Represent the mutation spectra with a bar graph #
eda59b985b1c Uploaded
iarc
parents: 6
diff changeset
11 #########################################################################################################################################
eda59b985b1c Uploaded
iarc
parents: 6
diff changeset
12
eda59b985b1c Uploaded
iarc
parents: 6
diff changeset
13 #-------------------------------------------------------------------------------
eda59b985b1c Uploaded
iarc
parents: 6
diff changeset
14 # Print a usage message if there is no argument pass to the command line
eda59b985b1c Uploaded
iarc
parents: 6
diff changeset
15 #-------------------------------------------------------------------------------
eda59b985b1c Uploaded
iarc
parents: 6
diff changeset
16 args <- commandArgs(TRUE)
eda59b985b1c Uploaded
iarc
parents: 6
diff changeset
17 usage <- function()
eda59b985b1c Uploaded
iarc
parents: 6
diff changeset
18 {
eda59b985b1c Uploaded
iarc
parents: 6
diff changeset
19 msg <- paste0('Usage:\n',
eda59b985b1c Uploaded
iarc
parents: 6
diff changeset
20 ' mutationSpectra_Galaxy.r input_Mutation_Spectra Sample_Name Output_Folder_High_Resolution Output_Folder_Low_Resolution Count_ca Count_cg Count_ta Count_tc Count_tg\n',
eda59b985b1c Uploaded
iarc
parents: 6
diff changeset
21 '\ninput_Mutation_Spectra should be tab-separated: alteration context value\n',
eda59b985b1c Uploaded
iarc
parents: 6
diff changeset
22 '\nOutput_Folder_High_Resolution: Folder for saving the high resolution image (display on the HTML page)\n',
eda59b985b1c Uploaded
iarc
parents: 6
diff changeset
23 '\nOutput_Folder_Low_Resolution: Folder for saving the low resolution image (display on the Excel report)\n'
eda59b985b1c Uploaded
iarc
parents: 6
diff changeset
24 )
eda59b985b1c Uploaded
iarc
parents: 6
diff changeset
25 cat(msg, '\n', file="/dev/stderr")
eda59b985b1c Uploaded
iarc
parents: 6
diff changeset
26 quit(status=1)
eda59b985b1c Uploaded
iarc
parents: 6
diff changeset
27 }
eda59b985b1c Uploaded
iarc
parents: 6
diff changeset
28
eda59b985b1c Uploaded
iarc
parents: 6
diff changeset
29 input = args[length(args)]
eda59b985b1c Uploaded
iarc
parents: 6
diff changeset
30
eda59b985b1c Uploaded
iarc
parents: 6
diff changeset
31 if (length(args) == 0) { usage() }
eda59b985b1c Uploaded
iarc
parents: 6
diff changeset
32
eda59b985b1c Uploaded
iarc
parents: 6
diff changeset
33
eda59b985b1c Uploaded
iarc
parents: 6
diff changeset
34 #-------------------------------------------------------------------------------
eda59b985b1c Uploaded
iarc
parents: 6
diff changeset
35 # Load library
eda59b985b1c Uploaded
iarc
parents: 6
diff changeset
36 #-------------------------------------------------------------------------------
eda59b985b1c Uploaded
iarc
parents: 6
diff changeset
37 suppressMessages(suppressWarnings(library(ggplot2)))
eda59b985b1c Uploaded
iarc
parents: 6
diff changeset
38 suppressMessages(suppressWarnings(library(reshape)))
eda59b985b1c Uploaded
iarc
parents: 6
diff changeset
39 suppressMessages(suppressWarnings(library(grid)))
eda59b985b1c Uploaded
iarc
parents: 6
diff changeset
40 suppressMessages(suppressWarnings(library(scales)))
eda59b985b1c Uploaded
iarc
parents: 6
diff changeset
41 suppressMessages(suppressWarnings(library(gridExtra)))
eda59b985b1c Uploaded
iarc
parents: 6
diff changeset
42
eda59b985b1c Uploaded
iarc
parents: 6
diff changeset
43
eda59b985b1c Uploaded
iarc
parents: 6
diff changeset
44
eda59b985b1c Uploaded
iarc
parents: 6
diff changeset
45 #-------------------------------------------------------------------------------
eda59b985b1c Uploaded
iarc
parents: 6
diff changeset
46 # Recover the arguments
eda59b985b1c Uploaded
iarc
parents: 6
diff changeset
47 #-------------------------------------------------------------------------------
eda59b985b1c Uploaded
iarc
parents: 6
diff changeset
48 input <- args[1]
eda59b985b1c Uploaded
iarc
parents: 6
diff changeset
49 sampleName <- args[2]
eda59b985b1c Uploaded
iarc
parents: 6
diff changeset
50 output_html <- args[3]
eda59b985b1c Uploaded
iarc
parents: 6
diff changeset
51 output_report <- args[4]
eda59b985b1c Uploaded
iarc
parents: 6
diff changeset
52 count_ca <- as.numeric(args[5])
eda59b985b1c Uploaded
iarc
parents: 6
diff changeset
53 count_cg <- as.numeric(args[6])
eda59b985b1c Uploaded
iarc
parents: 6
diff changeset
54 count_ct <- as.numeric(args[7])
eda59b985b1c Uploaded
iarc
parents: 6
diff changeset
55 count_ta <- as.numeric(args[8])
eda59b985b1c Uploaded
iarc
parents: 6
diff changeset
56 count_tc <- as.numeric(args[9])
eda59b985b1c Uploaded
iarc
parents: 6
diff changeset
57 count_tg <- as.numeric(args[10])
eda59b985b1c Uploaded
iarc
parents: 6
diff changeset
58
eda59b985b1c Uploaded
iarc
parents: 6
diff changeset
59 count_ca <- paste("C>A (", count_ca, ")")
eda59b985b1c Uploaded
iarc
parents: 6
diff changeset
60 count_cg <- paste("C>G (", count_cg, ")")
eda59b985b1c Uploaded
iarc
parents: 6
diff changeset
61 count_ct <- paste("C>T (", count_ct, ")")
eda59b985b1c Uploaded
iarc
parents: 6
diff changeset
62 count_ta <- paste("T>A (", count_ta, ")")
eda59b985b1c Uploaded
iarc
parents: 6
diff changeset
63 count_tc <- paste("T>C (", count_tc, ")")
eda59b985b1c Uploaded
iarc
parents: 6
diff changeset
64 count_tg <- paste("T>G (", count_tg, ")")
eda59b985b1c Uploaded
iarc
parents: 6
diff changeset
65
eda59b985b1c Uploaded
iarc
parents: 6
diff changeset
66
eda59b985b1c Uploaded
iarc
parents: 6
diff changeset
67
eda59b985b1c Uploaded
iarc
parents: 6
diff changeset
68 ###############################################################################
eda59b985b1c Uploaded
iarc
parents: 6
diff changeset
69 # Load the functions #
eda59b985b1c Uploaded
iarc
parents: 6
diff changeset
70 ###############################################################################
eda59b985b1c Uploaded
iarc
parents: 6
diff changeset
71
eda59b985b1c Uploaded
iarc
parents: 6
diff changeset
72 #-------------------------------------------------------------------------------
eda59b985b1c Uploaded
iarc
parents: 6
diff changeset
73 # Set the font depending on X11 availability
eda59b985b1c Uploaded
iarc
parents: 6
diff changeset
74 #-------------------------------------------------------------------------------
eda59b985b1c Uploaded
iarc
parents: 6
diff changeset
75 font <- ""
eda59b985b1c Uploaded
iarc
parents: 6
diff changeset
76 # Check the device available
eda59b985b1c Uploaded
iarc
parents: 6
diff changeset
77 device <- capabilities()
eda59b985b1c Uploaded
iarc
parents: 6
diff changeset
78 # X11 is available
eda59b985b1c Uploaded
iarc
parents: 6
diff changeset
79 if(device[5]) { font <- "Helvetica" } else { font <- "mono" }
eda59b985b1c Uploaded
iarc
parents: 6
diff changeset
80
eda59b985b1c Uploaded
iarc
parents: 6
diff changeset
81 #-------------------------------------------------------------------------------
eda59b985b1c Uploaded
iarc
parents: 6
diff changeset
82 # My own thme
eda59b985b1c Uploaded
iarc
parents: 6
diff changeset
83 #-------------------------------------------------------------------------------
eda59b985b1c Uploaded
iarc
parents: 6
diff changeset
84 theme_custom <- function(base_size = 12, base_family = "")
eda59b985b1c Uploaded
iarc
parents: 6
diff changeset
85 {
eda59b985b1c Uploaded
iarc
parents: 6
diff changeset
86 # Starts with theme_grey and then modify some parts
eda59b985b1c Uploaded
iarc
parents: 6
diff changeset
87 theme_grey(base_size = base_size, base_family = base_family) %+replace%
eda59b985b1c Uploaded
iarc
parents: 6
diff changeset
88 theme(
eda59b985b1c Uploaded
iarc
parents: 6
diff changeset
89 axis.text = element_text(size = rel(1), family=font),
eda59b985b1c Uploaded
iarc
parents: 6
diff changeset
90 axis.ticks = element_line(colour = "black"),
eda59b985b1c Uploaded
iarc
parents: 6
diff changeset
91 axis.line = element_line(colour = "black", size = .5),
eda59b985b1c Uploaded
iarc
parents: 6
diff changeset
92 legend.key = element_blank(),
eda59b985b1c Uploaded
iarc
parents: 6
diff changeset
93 panel.background = element_blank(),
eda59b985b1c Uploaded
iarc
parents: 6
diff changeset
94 panel.border = element_blank(),
eda59b985b1c Uploaded
iarc
parents: 6
diff changeset
95 panel.grid.major = element_blank(),
eda59b985b1c Uploaded
iarc
parents: 6
diff changeset
96 panel.grid.minor = element_blank()
eda59b985b1c Uploaded
iarc
parents: 6
diff changeset
97 )
eda59b985b1c Uploaded
iarc
parents: 6
diff changeset
98 }
eda59b985b1c Uploaded
iarc
parents: 6
diff changeset
99
eda59b985b1c Uploaded
iarc
parents: 6
diff changeset
100
eda59b985b1c Uploaded
iarc
parents: 6
diff changeset
101 #-------------------------------------------------------------------------------
eda59b985b1c Uploaded
iarc
parents: 6
diff changeset
102 # Customize the theme for adding a y axis
eda59b985b1c Uploaded
iarc
parents: 6
diff changeset
103 #-------------------------------------------------------------------------------
eda59b985b1c Uploaded
iarc
parents: 6
diff changeset
104 mytheme <- theme_custom()
eda59b985b1c Uploaded
iarc
parents: 6
diff changeset
105 mytheme$axis.line.x <- mytheme$axis.line.y <- mytheme$axis.line
eda59b985b1c Uploaded
iarc
parents: 6
diff changeset
106 mytheme$axis.line.x$colour <- 'white'
eda59b985b1c Uploaded
iarc
parents: 6
diff changeset
107
eda59b985b1c Uploaded
iarc
parents: 6
diff changeset
108 #-------------------------------------------------------------------------------
eda59b985b1c Uploaded
iarc
parents: 6
diff changeset
109 # Set the decimal precision to 0.0
eda59b985b1c Uploaded
iarc
parents: 6
diff changeset
110 #------------------------------------------------------------------------------
eda59b985b1c Uploaded
iarc
parents: 6
diff changeset
111 fmt <- function()
eda59b985b1c Uploaded
iarc
parents: 6
diff changeset
112 {
eda59b985b1c Uploaded
iarc
parents: 6
diff changeset
113 function(x) format(x,nsmall = 1,scientific = FALSE, digits=1)
eda59b985b1c Uploaded
iarc
parents: 6
diff changeset
114 }
eda59b985b1c Uploaded
iarc
parents: 6
diff changeset
115
eda59b985b1c Uploaded
iarc
parents: 6
diff changeset
116
eda59b985b1c Uploaded
iarc
parents: 6
diff changeset
117
eda59b985b1c Uploaded
iarc
parents: 6
diff changeset
118 ###############################################################################
eda59b985b1c Uploaded
iarc
parents: 6
diff changeset
119 # MAIN #
eda59b985b1c Uploaded
iarc
parents: 6
diff changeset
120 ###############################################################################
eda59b985b1c Uploaded
iarc
parents: 6
diff changeset
121
eda59b985b1c Uploaded
iarc
parents: 6
diff changeset
122 matrixW_inputggplot2 <- read.table(input, header=T)
eda59b985b1c Uploaded
iarc
parents: 6
diff changeset
123 matrixW_melt <- melt(matrixW_inputggplot2)
eda59b985b1c Uploaded
iarc
parents: 6
diff changeset
124 max_matrixW <- max(matrixW_inputggplot2[,3:ncol(matrixW_inputggplot2)])
eda59b985b1c Uploaded
iarc
parents: 6
diff changeset
125
eda59b985b1c Uploaded
iarc
parents: 6
diff changeset
126
eda59b985b1c Uploaded
iarc
parents: 6
diff changeset
127 p <- ggplot(matrixW_melt, aes(x=context, y=value, fill=alteration)) + geom_bar(stat="identity", width=0.5) + facet_grid(variable ~ alteration, scales="free_y")
eda59b985b1c Uploaded
iarc
parents: 6
diff changeset
128 # Color the mutation like Alexandrov et al.
eda59b985b1c Uploaded
iarc
parents: 6
diff changeset
129 p <- p + scale_fill_manual(values=c("blue", "black", "red", "#828282", "#00CC33", "pink"))
eda59b985b1c Uploaded
iarc
parents: 6
diff changeset
130 # Remove the legend
eda59b985b1c Uploaded
iarc
parents: 6
diff changeset
131 p <- p + guides(fill=FALSE)
eda59b985b1c Uploaded
iarc
parents: 6
diff changeset
132 # customized theme (no background, no facet grid and strip, y axis only)
eda59b985b1c Uploaded
iarc
parents: 6
diff changeset
133 p <- p + mytheme
eda59b985b1c Uploaded
iarc
parents: 6
diff changeset
134 # Remove the title of the x facet strip
eda59b985b1c Uploaded
iarc
parents: 6
diff changeset
135 p <- p + theme(strip.text.x=element_blank(), strip.text.y=element_blank())
eda59b985b1c Uploaded
iarc
parents: 6
diff changeset
136 # Remove the x axis label, thicks and title
eda59b985b1c Uploaded
iarc
parents: 6
diff changeset
137 p <- p + theme(axis.title.x=element_blank(), axis.ticks.x=element_blank(), axis.title.y=element_text(size=15))
eda59b985b1c Uploaded
iarc
parents: 6
diff changeset
138 # Scale the y axis to the maximum value
eda59b985b1c Uploaded
iarc
parents: 6
diff changeset
139 p <- p + scale_y_continuous(limits=c(0,max_matrixW), oob=squish, breaks=c(0,max_matrixW), labels=fmt())
eda59b985b1c Uploaded
iarc
parents: 6
diff changeset
140 # Rename the y axis
eda59b985b1c Uploaded
iarc
parents: 6
diff changeset
141 p <- p + ylab("percent")
eda59b985b1c Uploaded
iarc
parents: 6
diff changeset
142 # Add a title to the plot
eda59b985b1c Uploaded
iarc
parents: 6
diff changeset
143 p <- p + ggtitle(sampleName) + theme(plot.title = element_text(vjust = 3.4, family=font))
eda59b985b1c Uploaded
iarc
parents: 6
diff changeset
144 # Add a top margin for writing the title of the plot
eda59b985b1c Uploaded
iarc
parents: 6
diff changeset
145 p <- p + theme(plot.margin=unit(c(.7,0,0,0), "cm"))
eda59b985b1c Uploaded
iarc
parents: 6
diff changeset
146 p <- p + scale_x_discrete(breaks = c("A_A","A_C","A_G","A_T", "C_A","C_C","C_G","C_T", "G_A","G_C","G_G","G_T", "T_A","T_C","T_G","T_T"),
eda59b985b1c Uploaded
iarc
parents: 6
diff changeset
147 labels =c('A\nA',"\nC","\nG","\nT", 'C\nA',"\nC","\nG","\nT",
eda59b985b1c Uploaded
iarc
parents: 6
diff changeset
148 'G\nA',"\nC","\nG","\nT", 'T\nA',"\nC","\nG","\nT"
eda59b985b1c Uploaded
iarc
parents: 6
diff changeset
149 )
eda59b985b1c Uploaded
iarc
parents: 6
diff changeset
150 )
eda59b985b1c Uploaded
iarc
parents: 6
diff changeset
151
eda59b985b1c Uploaded
iarc
parents: 6
diff changeset
152 #------------------------------------------------------------------------------------------------------------------------------
eda59b985b1c Uploaded
iarc
parents: 6
diff changeset
153 # Change the color of the facets for the mutation type
eda59b985b1c Uploaded
iarc
parents: 6
diff changeset
154 #------------------------------------------------------------------------------------------------------------------------------
eda59b985b1c Uploaded
iarc
parents: 6
diff changeset
155 cols <- rep( c("blue", "black", "red", "#828282", "#00CC33", "pink")) # Facet strip colours
eda59b985b1c Uploaded
iarc
parents: 6
diff changeset
156
eda59b985b1c Uploaded
iarc
parents: 6
diff changeset
157 # Make a grob object
eda59b985b1c Uploaded
iarc
parents: 6
diff changeset
158 Pg <- ggplotGrob(p)
eda59b985b1c Uploaded
iarc
parents: 6
diff changeset
159 # To keep track of strip.background grobs
eda59b985b1c Uploaded
iarc
parents: 6
diff changeset
160 idx <- 0
eda59b985b1c Uploaded
iarc
parents: 6
diff changeset
161 # Find each strip.background and alter its backround colour
eda59b985b1c Uploaded
iarc
parents: 6
diff changeset
162 for( g in 1:length(Pg$grobs) )
eda59b985b1c Uploaded
iarc
parents: 6
diff changeset
163 {
eda59b985b1c Uploaded
iarc
parents: 6
diff changeset
164 if( grepl( "strip.absoluteGrob" , Pg$grobs[[g]]$name ) )
eda59b985b1c Uploaded
iarc
parents: 6
diff changeset
165 {
eda59b985b1c Uploaded
iarc
parents: 6
diff changeset
166 idx <- idx + 1
eda59b985b1c Uploaded
iarc
parents: 6
diff changeset
167 sb <- which( grepl( "strip\\.background" , names( Pg$grobs[[g]]$children ) ) )
eda59b985b1c Uploaded
iarc
parents: 6
diff changeset
168 Pg$grobs[[g]]$children[[sb]][]$gp$fill <- cols[idx]
eda59b985b1c Uploaded
iarc
parents: 6
diff changeset
169 }
eda59b985b1c Uploaded
iarc
parents: 6
diff changeset
170 }
eda59b985b1c Uploaded
iarc
parents: 6
diff changeset
171
eda59b985b1c Uploaded
iarc
parents: 6
diff changeset
172 # Reduce the size of the facet strip
eda59b985b1c Uploaded
iarc
parents: 6
diff changeset
173 Pg$heights[[3]] = unit(.1,"cm")
eda59b985b1c Uploaded
iarc
parents: 6
diff changeset
174
eda59b985b1c Uploaded
iarc
parents: 6
diff changeset
175
eda59b985b1c Uploaded
iarc
parents: 6
diff changeset
176 # Save the plot for the HTML page (higher resolution)
eda59b985b1c Uploaded
iarc
parents: 6
diff changeset
177 graphics.off() # close graphics windows
eda59b985b1c Uploaded
iarc
parents: 6
diff changeset
178 # Use cairo device as isn't possible to install X11 on the server...
eda59b985b1c Uploaded
iarc
parents: 6
diff changeset
179 png(paste0(output_html, "/", sampleName, "-MutationSpectraPercent-Genomic.png"), width=3500, heigh=500, res=300, type=c("cairo-png"))
eda59b985b1c Uploaded
iarc
parents: 6
diff changeset
180 plot(Pg)
eda59b985b1c Uploaded
iarc
parents: 6
diff changeset
181 ## Add label for the mutation type above the strip facet
eda59b985b1c Uploaded
iarc
parents: 6
diff changeset
182 grid.text(0.13, unit(0.90,"npc") - unit(1,"line"), label=count_ca)
eda59b985b1c Uploaded
iarc
parents: 6
diff changeset
183 grid.text(0.29, unit(0.90,"npc") - unit(1,"line"), label=count_cg)
eda59b985b1c Uploaded
iarc
parents: 6
diff changeset
184 grid.text(0.45, unit(0.90,"npc") - unit(1,"line"), label=count_ct)
eda59b985b1c Uploaded
iarc
parents: 6
diff changeset
185 grid.text(0.6, unit(0.90,"npc") - unit(1,"line"), label=count_ta)
eda59b985b1c Uploaded
iarc
parents: 6
diff changeset
186 grid.text(0.76, unit(0.90,"npc") - unit(1,"line"), label=count_tc)
eda59b985b1c Uploaded
iarc
parents: 6
diff changeset
187 grid.text(0.92, unit(0.90,"npc") - unit(1,"line"), label=count_tg)
eda59b985b1c Uploaded
iarc
parents: 6
diff changeset
188 invisible( dev.off() )
eda59b985b1c Uploaded
iarc
parents: 6
diff changeset
189
eda59b985b1c Uploaded
iarc
parents: 6
diff changeset
190 # Save the plot for the report
eda59b985b1c Uploaded
iarc
parents: 6
diff changeset
191 png(paste0(output_report, "/", sampleName, "-MutationSpectraPercent-Genomic-Report.png"), width=1000, heigh=150, type=c("cairo-png"))
eda59b985b1c Uploaded
iarc
parents: 6
diff changeset
192 plot(Pg)
eda59b985b1c Uploaded
iarc
parents: 6
diff changeset
193 ## Add label for the mutation type above the strip facet
eda59b985b1c Uploaded
iarc
parents: 6
diff changeset
194 grid.text(0.13, unit(0.90,"npc") - unit(1,"line"), label=count_ca)
eda59b985b1c Uploaded
iarc
parents: 6
diff changeset
195 grid.text(0.29, unit(0.90,"npc") - unit(1,"line"), label=count_cg)
eda59b985b1c Uploaded
iarc
parents: 6
diff changeset
196 grid.text(0.45, unit(0.90,"npc") - unit(1,"line"), label=count_ct)
eda59b985b1c Uploaded
iarc
parents: 6
diff changeset
197 grid.text(0.6, unit(0.90,"npc") - unit(1,"line"), label=count_ta)
eda59b985b1c Uploaded
iarc
parents: 6
diff changeset
198 grid.text(0.76, unit(0.90,"npc") - unit(1,"line"), label=count_tc)
eda59b985b1c Uploaded
iarc
parents: 6
diff changeset
199 grid.text(0.92, unit(0.90,"npc") - unit(1,"line"), label=count_tg)
eda59b985b1c Uploaded
iarc
parents: 6
diff changeset
200 invisible( dev.off() )
eda59b985b1c Uploaded
iarc
parents: 6
diff changeset
201
eda59b985b1c Uploaded
iarc
parents: 6
diff changeset
202 # Delete the empty plot created by the script
eda59b985b1c Uploaded
iarc
parents: 6
diff changeset
203 if (file.exists("Rplots.pdf")) invisible( file.remove("Rplots.pdf") )