|
1
|
1 suppressMessages(library(BoutrosLab.plotting.general));
|
|
|
2 suppressMessages(library(dplyr));
|
|
|
3
|
|
|
4 ### GET ARGUMENTS ##################################################################################
|
|
|
5
|
|
|
6 suppressMessages(library(argparse));
|
|
|
7 parser <- ArgumentParser(description="Create a Gene by Patient Heatmap BPGs multiplot fucntion");
|
|
|
8
|
|
|
9 ### VARIANT INPUT ARGUMENTS ########################################################################
|
|
|
10
|
|
|
11 parser$add_argument(
|
|
|
12 "--input_snv", "-snv",
|
|
|
13 required="True",
|
|
|
14 help="Input SNV files in MAF format"
|
|
|
15 );
|
|
|
16
|
|
|
17 parser$add_argument(
|
|
|
18 "--input_cnv", "-cnv",
|
|
|
19 help="Input CNV files in MAF format"
|
|
|
20 );
|
|
|
21
|
|
|
22
|
|
|
23 ### GENE INPUT ARGUMENTS ###########################################################################
|
|
|
24
|
|
|
25 parser$add_argument(
|
|
|
26 "--gene_list", "-gl",
|
|
|
27 help=""
|
|
|
28 );
|
|
|
29
|
|
|
30 parser$add_argument(
|
|
|
31 "--gene_oncodrive", "-go",
|
|
|
32 help=""
|
|
|
33 );
|
|
|
34
|
|
|
35 parser$add_argument(
|
|
|
36 "--gene_mutsig", "-gc",
|
|
|
37 help=""
|
|
|
38 )
|
|
|
39
|
|
|
40 ### GENERAL ARGUMENTS ##############################################################################
|
|
|
41
|
|
|
42 parser$add_argument(
|
|
|
43 "--center_plot",
|
|
|
44 choices=c('impact', 'variant_classification'),
|
|
|
45 required=TRUE,
|
|
|
46 help="What should be displayed in the Center Heatmap?"
|
|
|
47 );
|
|
|
48
|
|
|
49 parser$add_argument(
|
|
|
50 "--cancer_genomics_functions_path", "-cgfp",
|
|
|
51 required='True',
|
|
|
52 help="Path to cancer genomics functions in R"
|
|
|
53 )
|
|
|
54
|
|
|
55 parser$add_argument(
|
|
|
56 "--output", "-o",
|
|
|
57 required="True",
|
|
|
58 help="Output Image File"
|
|
|
59 )
|
|
|
60
|
|
|
61 ### PATIENT ARGUMENTS ##############################################################################
|
|
|
62
|
|
|
63 parser$add_argument(
|
|
|
64 "--patient_covariate_data",
|
|
|
65 help="Input Patient Covariate Matrix"
|
|
|
66 );
|
|
|
67
|
|
|
68 parser$add_argument(
|
|
|
69 "--plot_patient_covariate",
|
|
|
70 action="append",
|
|
|
71 help="Should We Plot the Input Covariate Data?"
|
|
|
72 );
|
|
|
73
|
|
|
74 parser$add_argument(
|
|
|
75 "--plot_patient_snv_counts",
|
|
|
76 action='store_true',
|
|
|
77 help="Should We Plot Patient SNV Counts?"
|
|
|
78 );
|
|
|
79
|
|
|
80 parser$add_argument(
|
|
|
81 "--plot_patient_snv_distribution",
|
|
|
82 action='store_true',
|
|
|
83 help="Should We Plot Patient SNV Distribution?"
|
|
|
84 );
|
|
|
85
|
|
|
86 parser$add_argument(
|
|
|
87 "--patient_order",
|
|
|
88 action="append",
|
|
|
89 help="How Should We Determine Patient Order?"
|
|
|
90 );
|
|
|
91
|
|
|
92
|
|
|
93 ### GENE ARGUMENTS #################################################################################
|
|
|
94
|
|
|
95 parser$add_argument(
|
|
|
96 "--plot_gene_snv_counts",
|
|
|
97 action="store_true",
|
|
|
98 help="Should We Plot Gene SNV Counts?"
|
|
|
99 );
|
|
|
100
|
|
|
101 parser$add_argument(
|
|
|
102 "--plot_gene_mutsig",
|
|
|
103 action="store_true",
|
|
|
104 help="Should We Plot Gene Significance?"
|
|
|
105 );
|
|
|
106
|
|
|
107 parser$add_argument(
|
|
|
108 "--plot_gene_oncodrive",
|
|
|
109 action="store_true",
|
|
|
110 help="How Should We Determine Gene Order?"
|
|
|
111 );
|
|
|
112
|
|
|
113 parser$add_argument(
|
|
|
114 "--gene_order",
|
|
|
115 choices=c('mutsig', 'oncodrive', 'snv_counts', 'gene_list'),
|
|
|
116 help="How Should We Determine Gene Order?"
|
|
|
117 );
|
|
|
118
|
|
|
119 ####################################################################################################
|
|
|
120
|
|
|
121 args <- parser$parse_args();
|
|
|
122
|
|
|
123 source(paste(args$cancer_genomics_functions_path, "functions.R", sep="/"));
|
|
|
124
|
|
|
125 ### ORGANIZE GENE DATA #############################################################################
|
|
|
126
|
|
|
127 snv_data <- read_vcf2maf_maf(args$input_snv);
|
|
|
128
|
|
|
129 snv_data$patient
|
|
|
130
|
|
|
131 snv_counts <- calculate_gene_snv_counts(snv_data);
|
|
|
132
|
|
|
133 patient_covariate_data <- calculate_patient_snv_counts(snv_data);
|
|
|
134
|
|
|
135 patient_snv_distribution <- calculate_patient_snv_frequency(snv_data);
|
|
|
136
|
|
|
137 if (!is.null(args$patient_covariate_data)) {
|
|
|
138 tmp_patient_covariate_data <- read_patient_covariate_data(args$patient_covariate_data);
|
|
|
139
|
|
|
140 if (!all(rownames(patient_covariate_data) %in% rownames(tmp_patient_covariate_data))) {
|
|
|
141 stop("Some patients identified in the MAF file are not present in patient covariate data file");
|
|
|
142 }
|
|
|
143 indexes <- which( rownames(patient_covariate_data) %in% rownames(tmp_patient_covariate_data))
|
|
|
144 current_names <- colnames(patient_covariate_data);
|
|
|
145 patient_covariate_data <- cbind(patient_covariate_data, tmp_patient_covariate_data[rownames(patient_covariate_data)[indexes],-1]);
|
|
|
146 colnames(patient_covariate_data) <- c(current_names, colnames(tmp_patient_covariate_data)[-1]);
|
|
|
147 }
|
|
|
148
|
|
|
149
|
|
|
150 ##################################
|
|
|
151
|
|
|
152 gene_list <- NULL;
|
|
|
153 mutsig <- NULL;
|
|
|
154 mutsig_sig <- NULL;
|
|
|
155 oncodrive <- NULL;
|
|
|
156 oncodrive_sig <- NULL;
|
|
|
157
|
|
|
158 if (!is.null(args$gene_list)) {
|
|
|
159 gene_list <- read_gene_list(args$gene_list);
|
|
|
160 }
|
|
|
161 if (!is.null(args$gene_mutsig)) {
|
|
|
162 mutsig <- read_mutsigCV_genes(args$gene_mutsig);
|
|
|
163 mutsig_sig <- filter_gene_list(mutsig);
|
|
|
164 };
|
|
|
165 if (!is.null(args$gene_oncodrive)) {
|
|
|
166 oncodrive <- read_oncodriveFM_genes(args$gene_oncodrive);
|
|
|
167 oncodrive_sig <- filter_gene_list(oncodrive);
|
|
|
168 };
|
|
|
169
|
|
|
170 genes_of_interest <- c(
|
|
|
171 gene_list$gene,
|
|
|
172 mutsig_sig$gene,
|
|
|
173 oncodrive_sig$gene
|
|
|
174 );
|
|
|
175
|
|
|
176 if (!is.null(args$gene_list)) {
|
|
|
177 genes_of_interest <- gene_list$gene
|
|
|
178 }
|
|
|
179
|
|
|
180 genes_of_interest <- unique(genes_of_interest);
|
|
|
181
|
|
|
182 if (is.null(genes_of_interest)) {
|
|
|
183 genes_of_interest <- calculate_gene_snv_counts(snv_data, best=TRUE)$gene
|
|
|
184 }
|
|
|
185
|
|
|
186 gene_data <- data.frame(
|
|
|
187 gene = genes_of_interest,
|
|
|
188 user = genes_of_interest %in% gene_list,
|
|
|
189 mutsig = genes_of_interest %in% mutsig_sig$gene,
|
|
|
190 oncodrive = genes_of_interest %in% oncodrive_sig$gene,
|
|
|
191 stringsAsFactors = F
|
|
|
192 );
|
|
|
193
|
|
|
194 gene_order <- NULL;
|
|
|
195
|
|
|
196 if (is.null(args$gene_order)) {
|
|
|
197 gene_order <- gene_data$gene;
|
|
|
198
|
|
|
199 } else if (args$gene_order == 'mutsig') {
|
|
|
200 add_to_end <- NULL;
|
|
|
201 if (!all(gene_data$gene %in% mutsig$gene)) {
|
|
|
202 warning("mutsigCV input file is missing some genes, assuming a p value of 1");
|
|
|
203 add_to_end <- gene_data$gene[!gene_data$gene %in% mutsig$gene]
|
|
|
204 }
|
|
|
205
|
|
|
206 gene_order <- c(mutsig$gene[mutsig$gene %in% gene_data$gene], add_to_end);
|
|
|
207
|
|
|
208 } else if (args$gene_order == 'oncodrive') {
|
|
|
209 add_to_end <- NULL;
|
|
|
210 if (!all(gene_data$gene %in% oncodrive$gene)) {
|
|
|
211 warning("oncodriveFM input file is missing some genes, assuming a p value of 1");
|
|
|
212 add_to_end <- gene_data$gene[!gene_data$gene %in% oncodrive$gene]
|
|
|
213 }
|
|
|
214
|
|
|
215 gene_order <- c(oncodrive$gene[oncodrive$gene %in% gene_data$gene], add_to_end);
|
|
|
216
|
|
|
217 } else if (args$gene_order == 'snv_counts') {
|
|
|
218 add_to_end <- NULL;
|
|
|
219 if ( !all(gene_data$gene %in% snv_counts$gene)) {
|
|
|
220 warning("snv input file is missing some genes, assuming snv count of 0");
|
|
|
221 add_to_end <- gene_data$gene[!gene_data$gene %in% snv_counts$gene]
|
|
|
222 }
|
|
|
223
|
|
|
224 gene_order <- c( snv_counts[snv_counts$gene %in% gene_data$gene,]$gene, add_to_end);
|
|
|
225 }
|
|
|
226
|
|
|
227 patient_order <- NULL;
|
|
|
228
|
|
|
229
|
|
|
230
|
|
|
231 if ( is.null(args$patient_order) ) {
|
|
|
232 patient_order <- unique(snv_data$patient);
|
|
|
233 patient_covariate_data <- patient_covariate_data[patient_order,];
|
|
|
234
|
|
|
235 } else {
|
|
|
236 for (covariate in args$patient_order) {
|
|
|
237 print(covariate)
|
|
|
238 print(colnames(patient_covariate_data))
|
|
|
239 current_order <- order(patient_covariate_data[,covariate]);
|
|
|
240 patient_covariate_data <- patient_covariate_data[current_order,];
|
|
|
241 }
|
|
|
242
|
|
|
243 patient_order <- rownames(patient_covariate_data);
|
|
|
244
|
|
|
245 }
|
|
|
246
|
|
|
247 ### GENE PLOTS #####################################################################################
|
|
|
248
|
|
|
249 gene_plots <- list();
|
|
|
250
|
|
|
251
|
|
|
252 if (args$plot_gene_snv_counts) {
|
|
|
253
|
|
|
254 plot_data <- data.frame(
|
|
|
255 gene = gene_order,
|
|
|
256 order = length(gene_order):1,
|
|
|
257 counts = log10(snv_counts[gene_order,]$counts)
|
|
|
258 );
|
|
|
259
|
|
|
260 max <- ceiling(max(plot_data$counts))
|
|
|
261 middle <- round(median(plot_data$counts))
|
|
|
262
|
|
|
263 plot <- create.barplot(
|
|
|
264 data = plot_data,
|
|
|
265 formula = order ~ counts,
|
|
|
266 plot.horizontal = TRUE,
|
|
|
267 xlimits = c(-0.1, max+0.1),
|
|
|
268 xat = c(0, middle, max),
|
|
|
269 xaxis.lab = c( expression(10^0), bquote(10^.(middle)), bquote(10^.(max)) ),
|
|
|
270 yat = 0
|
|
|
271 );
|
|
|
272
|
|
|
273 gene_plots <- c(gene_plots, list(plot));
|
|
|
274
|
|
|
275 }
|
|
|
276
|
|
|
277
|
|
|
278 if (args$plot_gene_mutsig) {
|
|
|
279
|
|
|
280 plot_data <- data.frame(
|
|
|
281 gene = gene_order,
|
|
|
282 order = length(gene_order):1,
|
|
|
283 p.value = -1*log10(mutsig[gene_order,]$p.value),
|
|
|
284 q.value = mutsig[gene_order,]$q.value,
|
|
|
285 stringsAsFactors=F
|
|
|
286 );
|
|
|
287
|
|
|
288 if (any( is.infinite(plot_data$p.value) ) ) {
|
|
|
289 plot_data$p.value[is.infinite(plot_data$p.value)] <- max( c( plot_data$p.value[!is.infinite(plot_data$p.value)], 10) );
|
|
|
290 }
|
|
|
291
|
|
|
292 plot_data$p.value[is.na(plot_data$p.value)] <- 1
|
|
|
293 plot_data$q.value[is.na(plot_data$q.value)] <- 1
|
|
|
294 plot_data$colour = c("black", "grey80")[(plot_data$q.value <= 0.05) + 1]
|
|
|
295
|
|
|
296 max <- ceiling(max(plot_data$p.value))
|
|
|
297 middle <- round(median(plot_data$p.value))
|
|
|
298
|
|
|
299 plot <- create.barplot(
|
|
|
300 data = plot_data,
|
|
|
301 formula = order ~ p.value,
|
|
|
302 col = rev(plot_data$colour),
|
|
|
303 plot.horizontal = TRUE,
|
|
|
304 xlimits = c(-1, max+1),
|
|
|
305 xat = c(0, middle, max),
|
|
|
306 xaxis.lab = c( expression(10^-0), bquote(10^-.(middle)), bquote(10^-.(max)) ) ,
|
|
|
307 yat = 0,
|
|
|
308 );
|
|
|
309
|
|
|
310 gene_plots <- c(gene_plots, list(plot))
|
|
|
311
|
|
|
312 }
|
|
|
313
|
|
|
314
|
|
|
315 if (args$plot_gene_oncodrive) {
|
|
|
316
|
|
|
317 plot_data <- data.frame(
|
|
|
318 gene = gene_order,
|
|
|
319 order = length(gene_order):1,
|
|
|
320 p.value = -1 * log10(oncodrive[gene_order,]$p.value),
|
|
|
321 q.value = oncodrive[gene_order,]$q.value,
|
|
|
322 stringsAsFactors=F
|
|
|
323 );
|
|
|
324
|
|
|
325 plot_data$p.value[is.na(plot_data$p.value)] <- 1
|
|
|
326 plot_data$q.value[is.na(plot_data$q.value)] <- 1
|
|
|
327 plot_data$colour = c("black", "grey80")[(plot_data$q.value <= 0.05) + 1]
|
|
|
328
|
|
|
329 max <- ceiling(max(plot_data$p.value))
|
|
|
330 middle <- round(median(plot_data$p.value))
|
|
|
331
|
|
|
332 plot <- create.barplot(
|
|
|
333 data = plot_data,
|
|
|
334 formula = order ~ p.value,
|
|
|
335 col = rev(plot_data$colour),
|
|
|
336 plot.horizontal = TRUE,
|
|
|
337 xlimits = c(-1, max),
|
|
|
338 xat = c(0, middle, max),
|
|
|
339 xaxis.lab = c( expression(10^-0), bquote(10^-.(middle)), bquote(10^-.(max)) ),
|
|
|
340 yat = 0
|
|
|
341 );
|
|
|
342
|
|
|
343 gene_plots <- c(gene_plots, list(plot))
|
|
|
344
|
|
|
345 }
|
|
|
346
|
|
|
347
|
|
|
348 ### PATIENT PLOTS ##################################################################################
|
|
|
349
|
|
|
350 patient_plots <- list();
|
|
|
351
|
|
|
352 if (args$plot_patient_snv_counts) {
|
|
|
353
|
|
|
354 plot_data <- data.frame(
|
|
|
355 patient = patient_order,
|
|
|
356 order = 1:length(patient_order),
|
|
|
357 counts = log10(patient_covariate_data[patient_order,]$snv_counts)
|
|
|
358 )
|
|
|
359
|
|
|
360 middle <- round(median(plot_data$counts))
|
|
|
361 max <- ceiling(max(plot_data$counts))
|
|
|
362
|
|
|
363 plot <- create.barplot(
|
|
|
364 data = plot_data,
|
|
|
365 formula = counts ~ order,
|
|
|
366 ylimits = c(-0.1, max+0.1),
|
|
|
367 yaxis.lab = c(expression(10^0), bquote(10^.(middle)), bquote(10^.(max)) ),
|
|
|
368 yat = c(0,middle,max),
|
|
|
369 xat=0
|
|
|
370 );
|
|
|
371
|
|
|
372 patient_plots <- c(patient_plots, list(plot));
|
|
|
373
|
|
|
374 }
|
|
|
375
|
|
|
376 if (args$plot_patient_snv_distribution) {
|
|
|
377
|
|
|
378 patient_to_order <- data.frame(
|
|
|
379 order = 1:length(patient_order)
|
|
|
380 );
|
|
|
381
|
|
|
382 rownames(patient_to_order) <- patient_order;
|
|
|
383
|
|
|
384 plot_data <- data.frame(
|
|
|
385 patient = patient_snv_distribution$patient,
|
|
|
386 category = patient_snv_distribution$category,
|
|
|
387 counts = patient_snv_distribution$counts,
|
|
|
388 order = patient_to_order[patient_snv_distribution$patient,],
|
|
|
389 colours = rep(c("#0e8efa", "#010101","#fe0000", "#bfbfbf", "#00f100", "#ffbfcb"), by=6),
|
|
|
390 stringsAsFactors=F
|
|
|
391 );
|
|
|
392
|
|
|
393 plot <- create.barplot(
|
|
|
394 data = plot_data,
|
|
|
395 formula = counts ~ order,
|
|
|
396 groups = category,
|
|
|
397 stack = TRUE,
|
|
|
398 xat = 0,
|
|
|
399 yaxis.lab = c("0.0", "0.2", "0.4", "0.6", "0.8", "1.0"),
|
|
|
400 ylimits = c(-0.1, 1.1),
|
|
|
401 yat = c(0, .2, .4, .6, .8, 1.0),
|
|
|
402 col = plot_data$colours
|
|
|
403 );
|
|
|
404
|
|
|
405 patient_plots <- c(patient_plots, list(plot));
|
|
|
406
|
|
|
407 }
|
|
|
408
|
|
|
409 ### PATIENT COVARIATE PLOT #########################################################################
|
|
|
410
|
|
|
411 patient_covariate_plots <- list();
|
|
|
412 patient_covariate_legend <- list();
|
|
|
413
|
|
|
414 if (!is.null(args$plot_patient_covariate)) {
|
|
|
415
|
|
|
416 cov.data <- patient_covariate_data[patient_order,c('snv_counts', args$plot_patient_covariate)];
|
|
|
417
|
|
|
418 cov.colors <- c();
|
|
|
419 for (cov in args$plot_patient_covariate) {
|
|
|
420 tmp <- paste0(cov, "_col");
|
|
|
421 if ( any(tmp == colnames(patient_covariate_data)) ) {
|
|
|
422 tmp.colors <- unique(patient_covariate_data[,tmp]);
|
|
|
423 for (col_index in (length(cov.colors)+1):((length(cov.colors)+1)+length(tmp.colors))) {
|
|
|
424 col <- tmp.colors[col_index - length(cov.colors)];
|
|
|
425 cov.data[which(col == patient_covariate_data[,tmp]), cov] <- col_index;
|
|
|
426 }
|
|
|
427 }
|
|
|
428 cov.colors <- c(cov.colors, tmp.colors);
|
|
|
429 }
|
|
|
430 cov.data <- cov.data[,-1];
|
|
|
431
|
|
|
432 covariates <- create.heatmap(
|
|
|
433 x = matrix(as.numeric(as.matrix(cov.data)), ncol=length(args$plot_patient_covariate)),
|
|
|
434 clustering.method = 'none',
|
|
|
435 colour.scheme = as.vector(cov.colors),
|
|
|
436 total.colours = length(cov.colors)+1,
|
|
|
437 row.colour = 'white',
|
|
|
438 col.colour = 'white',
|
|
|
439 grid.row = TRUE,
|
|
|
440 grid.col = TRUE,
|
|
|
441 yaxis.lab = args$plot_patient_covariate,
|
|
|
442 yat = 1:length(args$plot_patient_covariate),
|
|
|
443 xaxis.lab = patient_order,
|
|
|
444 xat = 1:length(patient_order),
|
|
|
445 yaxis.tck = 0,
|
|
|
446 print.colour.key = FALSE
|
|
|
447 );
|
|
|
448
|
|
|
449 patient_covariate_plots <- c(patient_covariate_plots, list(covariates));
|
|
|
450
|
|
|
451 for ( cov in args$plot_patient_covariate) {
|
|
|
452
|
|
|
453 legend <- list(
|
|
|
454 legend = list(
|
|
|
455 colours = rev(unique(patient_covariate_data[, paste0(cov, "_col")])),
|
|
|
456 labels = rev(unique(patient_covariate_data[, cov])),
|
|
|
457 border = 'white',
|
|
|
458 title = cov,
|
|
|
459 pch = 15
|
|
|
460 )
|
|
|
461 );
|
|
|
462
|
|
|
463 patient_covariate_legend <- c(patient_covariate_legend, legend);
|
|
|
464
|
|
|
465 }
|
|
|
466
|
|
|
467 }
|
|
|
468
|
|
|
469 ### CENTRAL PLOT ###################################################################################
|
|
|
470
|
|
|
471 central_plot <- NULL;
|
|
|
472
|
|
|
473 if ( args$center_plot == 'impact' ) {
|
|
|
474
|
|
|
475 plot_data <- matrix(
|
|
|
476 data = rep(0, length(patient_order) * length(gene_order)),
|
|
|
477 nrow = length(gene_order)
|
|
|
478 )
|
|
|
479
|
|
|
480 rownames(plot_data) <- rev(gene_order);
|
|
|
481 colnames(plot_data) <- patient_order;
|
|
|
482
|
|
|
483 for (patient in patient_order) {
|
|
|
484 for (gene in gene_order) {
|
|
|
485 truth <- snv_data$patient == patient & snv_data$gene == gene;
|
|
|
486 if (any(truth)) {
|
|
|
487 if (any(snv_data$impact[truth] == "HIGH")) {
|
|
|
488 plot_data[gene, patient] <- 3;
|
|
|
489 } else if (any(snv_data$impact[truth] == "MODERATE")) {
|
|
|
490 plot_data[gene, patient] <- 2;
|
|
|
491 } else if (any(snv_data$impact[truth] == "LOW")) {
|
|
|
492 plot_data[gene, patient] <- 1;
|
|
|
493 }
|
|
|
494 }
|
|
|
495 }
|
|
|
496 }
|
|
|
497
|
|
|
498 gene.labs <- rownames(plot_data);
|
|
|
499 patient.labs <- colnames(plot_data);
|
|
|
500 if (!is.null(args$plot_patient_covariate)) {
|
|
|
501 patient.labs <- rep("", length(patient.labs));
|
|
|
502 }
|
|
|
503
|
|
|
504 plot <- create.heatmap(
|
|
|
505 t(plot_data),
|
|
|
506 clustering.method = "none",
|
|
|
507 colour.scheme = c('grey95',c('#ffeda0', '#feb24c','#f03b20')),
|
|
|
508 total.colours = 5,
|
|
|
509 xaxis.tck = 0,
|
|
|
510 yaxis.tck = 0,
|
|
|
511 print.colour.key = F,
|
|
|
512 xaxis.lab = patient.labs,
|
|
|
513 yaxis.lab = gene.labs,
|
|
|
514 row.colour = 'white',
|
|
|
515 col.colour = 'white',
|
|
|
516 grid.row = TRUE,
|
|
|
517 grid.col = TRUE
|
|
|
518 );
|
|
|
519
|
|
|
520 central_plot <- list(plot);
|
|
|
521
|
|
|
522 }
|
|
|
523
|
|
|
524 which_present <- rep(FALSE, 9)
|
|
|
525
|
|
|
526 if (args$center_plot == 'variant_classification') {
|
|
|
527
|
|
|
528 plot_data <- matrix(
|
|
|
529 data = rep(0, length(patient_order) * length(gene_order)),
|
|
|
530 nrow = length(gene_order)
|
|
|
531 )
|
|
|
532
|
|
|
533 rownames(plot_data) <- rev(gene_order);
|
|
|
534 colnames(plot_data) <- patient_order;
|
|
|
535
|
|
|
536 order.of.variants <- data.frame(
|
|
|
537 Truncation = 'Truncation',
|
|
|
538 Splicing = 'Splicing',
|
|
|
539 Missense = 'Missense',
|
|
|
540 UTR = 'UTR',
|
|
|
541 Non_coding = "Non_coding",
|
|
|
542 Regulatory = "Regulatory",
|
|
|
543 miRNA = "miRNA",
|
|
|
544 Synonymous = "Synonymous",
|
|
|
545 stringsAsFactors = F
|
|
|
546 )
|
|
|
547
|
|
|
548 categories <- recode(
|
|
|
549 snv_data[,]$ensembl_variant_class,
|
|
|
550 # Truncating variants
|
|
|
551 stop_gained = "Truncation",
|
|
|
552 frameshift_variant = "Truncation",
|
|
|
553 stop_lost = "Truncation",
|
|
|
554 start_lost = "Truncation",
|
|
|
555 incomplete_terminal_codon_variant = "Truncation",
|
|
|
556 transcript_ablation = "Truncation",
|
|
|
557 # Missense(-ish) variants
|
|
|
558 missense_variant = "Missense",
|
|
|
559 inframe_deletion = "Missense",
|
|
|
560 inframe_insertion = "Missense",
|
|
|
561 protein_altering_variant = "Missense",
|
|
|
562 # UTR Variants
|
|
|
563 `5_prime_UTR_variant` = "UTR",
|
|
|
564 `3_prime_UTR_variant` = "UTR",
|
|
|
565 # Splicing variants
|
|
|
566 splice_acceptor_variant = "Splicing",
|
|
|
567 splice_donor_variant = "Splicing",
|
|
|
568 splice_region_variant = "Splicing",
|
|
|
569 # Non-coding variants
|
|
|
570 non_coding_transcript_exon_variant = "Non_coding",
|
|
|
571 intron_variant = "Non_coding",
|
|
|
572 upstream_gene_variant = "Non_coding",
|
|
|
573 downstream_gene_variant = "Non_coding",
|
|
|
574 intergenic_variant = "Non_coding",
|
|
|
575 # Regulatory variants
|
|
|
576 TF_binding_site_variant = "Regulatory",
|
|
|
577 regulatory_region_variant = "Regulatory",
|
|
|
578 # Synonymous variants
|
|
|
579 synonymous_variant = "Synonymous",
|
|
|
580 coding_sequence_variant = "Synonymous",
|
|
|
581 stop_retained_variant = "Synonymous",
|
|
|
582 # miRNA variants
|
|
|
583 mature_miRNA_variant = "miRNA",
|
|
|
584 # Default
|
|
|
585 .default = NA_character_)
|
|
|
586
|
|
|
587 for( gene in gene_order ) {
|
|
|
588 for( patient in patient_order ) {
|
|
|
589 truth <- snv_data$patient == patient & snv_data$gene == gene;
|
|
|
590 if (gene == "SRSF7" && any(truth)) {
|
|
|
591 }
|
|
|
592 if (any(truth)) {
|
|
|
593 vars <- colnames(order.of.variants) %in% categories[truth]
|
|
|
594 if(any(vars)) {
|
|
|
595 plot_data[gene, patient] <- which(vars)[1]
|
|
|
596 which_present[which(vars)[1]] <- TRUE
|
|
|
597 }
|
|
|
598 else {
|
|
|
599 which_present[9] <- TRUE
|
|
|
600 }
|
|
|
601 }
|
|
|
602 }
|
|
|
603 }
|
|
|
604
|
|
|
605 col_one <- "#512C6F";
|
|
|
606 col_two <- "#0F6A99";
|
|
|
607 col_thr <- "#46823C";
|
|
|
608 col_fou <- "#B367A7";
|
|
|
609 col_fiv <- "#64B4D5";
|
|
|
610 col_six <- "#F7D72E";
|
|
|
611 col_sev <- "#EF922A";
|
|
|
612 col_eig <- "#B12025";
|
|
|
613
|
|
|
614 gene.labs <- rownames(plot_data);
|
|
|
615 patient.labs <- colnames(plot_data);
|
|
|
616 if (!is.null(args$plot_patient_covariate)) {
|
|
|
617 patient.labs <- rep("", length(patient.labs));
|
|
|
618 }
|
|
|
619
|
|
|
620 plot <- create.heatmap(
|
|
|
621 t(plot_data),
|
|
|
622 clustering.method = "none",
|
|
|
623 colour.scheme = c('grey95', c(col_one, col_two, col_thr, col_fou, col_fiv, col_six, col_sev, col_eig) ),
|
|
|
624 total.colours = 10,
|
|
|
625 xaxis.tck = 0,
|
|
|
626 yaxis.tck = 0,
|
|
|
627 print.colour.key = F,
|
|
|
628 xaxis.lab = patient.labs,
|
|
|
629 yaxis.lab = gene.labs,
|
|
|
630 row.colour = 'white',
|
|
|
631 col.colour = 'white',
|
|
|
632 grid.row = TRUE,
|
|
|
633 grid.col = TRUE
|
|
|
634 );
|
|
|
635
|
|
|
636 central_plot <- list(plot);
|
|
|
637
|
|
|
638 }
|
|
|
639
|
|
|
640 ### LEGENDS ########################################################################################
|
|
|
641
|
|
|
642 legends <- list();
|
|
|
643
|
|
|
644 if (args$plot_patient_snv_distribution ) {
|
|
|
645
|
|
|
646 legend <- list(
|
|
|
647 legend = list(
|
|
|
648 colours = rev(c("#0e8efa", "#010101","#fe0000", "#bfbfbf", "#00f100", "#ffbfcb" )),
|
|
|
649 labels = rev(expression(C > A,C > G,C > T, T > A, T > C, T > G)),
|
|
|
650 border = 'white',
|
|
|
651 title = 'Mutation Frequencies',
|
|
|
652 pch = 15
|
|
|
653 )
|
|
|
654 );
|
|
|
655
|
|
|
656 legends <- c(legends, legend);
|
|
|
657
|
|
|
658 }
|
|
|
659
|
|
|
660 if (args$plot_gene_mutsig || args$plot_gene_oncodrive) {
|
|
|
661
|
|
|
662 legend <- list(
|
|
|
663 legend = list(
|
|
|
664 colours = c('grey80', 'black'),
|
|
|
665 labels = expression(Q <= 10^-1, Q > 10^-1),
|
|
|
666 border = 'white',
|
|
|
667 title = 'Q-value Groupings',
|
|
|
668 pch=15
|
|
|
669 )
|
|
|
670 )
|
|
|
671
|
|
|
672 legends <- c(legends, legend);
|
|
|
673
|
|
|
674 }
|
|
|
675
|
|
|
676 if ( args$center_plot == 'impact' ) {
|
|
|
677
|
|
|
678 legend <- list(
|
|
|
679 legend = list(
|
|
|
680 colours = c('grey95',c('#ffeda0', '#feb24c','#f03b20')),
|
|
|
681 labels = c("None", "Low", "Moderate", "High"),
|
|
|
682 border = 'white',
|
|
|
683 title = 'SNV Impact',
|
|
|
684 pch = 15
|
|
|
685 )
|
|
|
686 )
|
|
|
687
|
|
|
688 legends <- c(legends, legend);
|
|
|
689
|
|
|
690 }
|
|
|
691
|
|
|
692 if ( args$center_plot == 'variant_classification' ) {
|
|
|
693
|
|
|
694 col_one <- "#512C6F";
|
|
|
695 col_two <- "#0F6A99";
|
|
|
696 col_thr <- "#46823C";
|
|
|
697 col_fou <- "#B367A7";
|
|
|
698 col_fiv <- "#64B4D5";
|
|
|
699 col_six <- "#F7D72E";
|
|
|
700 col_sev <- "#EF922A";
|
|
|
701 col_eig <- "#B12025";
|
|
|
702
|
|
|
703 legend <- list(
|
|
|
704 legend = list(
|
|
|
705 colours = c(c(col_one, col_two, col_thr, col_fou, col_fiv, col_six, col_sev, col_eig), 'grey95' )[which_present],
|
|
|
706 labels = c("Truncation", "Splicing", "Missense", "UTR", "Non-coding", "Regulatory", "miRNA", "Synonymous", "None")[which_present],
|
|
|
707 border = 'white',
|
|
|
708 title = 'SNV Category',
|
|
|
709 pch = 15
|
|
|
710 )
|
|
|
711 )
|
|
|
712
|
|
|
713 legends <- c(legends, legend);
|
|
|
714
|
|
|
715 }
|
|
|
716
|
|
|
717 if (length(patient_covariate_legend) > 0) {
|
|
|
718
|
|
|
719 legends <- c(legends, patient_covariate_legend);
|
|
|
720
|
|
|
721 }
|
|
|
722
|
|
|
723 legend_grob <- legend.grob(
|
|
|
724 legends = legends,
|
|
|
725 title.just = 'left',
|
|
|
726 label.cex = 1.0,
|
|
|
727 title.cex = 1.0
|
|
|
728 );
|
|
|
729
|
|
|
730 ### MULTIPLOT ######################################################################################
|
|
|
731
|
|
|
732 plots <- c(
|
|
|
733 patient_covariate_plots,
|
|
|
734 central_plot,
|
|
|
735 gene_plots,
|
|
|
736 patient_plots
|
|
|
737 );
|
|
|
738
|
|
|
739 plot.layout <- c(1 + length(gene_plots), 1 + length(patient_plots));
|
|
|
740 layout.skip <- c(F, rep(F, length(gene_plots)), rep(c(F, rep(T, length(gene_plots))), length(patient_plots)));
|
|
|
741 panel.heights <- c(rep(3, length(patient_plots)),20);
|
|
|
742 panel.widths <- c(20,rep(3, length(gene_plots)));
|
|
|
743 if (length(patient_covariate_plots) > 0) {
|
|
|
744 plot.layout[2] <- plot.layout[2] + 1;
|
|
|
745 layout.skip <- c(F, rep(T, length(gene_plots)), layout.skip);
|
|
|
746 panel.heights <- c(panel.heights, 0.5 * length(args$plot_patient_covariate) );
|
|
|
747 }
|
|
|
748
|
|
|
749 create.multiplot(
|
|
|
750 filename = args$output,
|
|
|
751 res = 300,
|
|
|
752 plot.objects = plots,
|
|
|
753 plot.layout = plot.layout,
|
|
|
754 layout.skip = layout.skip,
|
|
|
755 panel.heights = panel.heights,
|
|
|
756 panel.widths = panel.widths,
|
|
|
757 height = 16,
|
|
|
758 width = 16,
|
|
|
759 x.spacing = -.8,
|
|
|
760 y.spacing = -.33 * max(nchar(patient_order)),
|
|
|
761 left.padding = 4,
|
|
|
762 xaxis.cex = .65,
|
|
|
763 yaxis.cex = .9,
|
|
|
764 xaxis.rot = 90,
|
|
|
765 retrieve.plot.labels=T,
|
|
|
766 legend = list(
|
|
|
767 right = list(
|
|
|
768 x = 0.10,
|
|
|
769 y = 0.50,
|
|
|
770 fun = legend_grob
|
|
|
771 )
|
|
|
772 ),
|
|
|
773
|
|
|
774 print.new.legend= TRUE
|
|
|
775
|
|
|
776 );
|
