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