Mercurial > repos > iuc > snpfreqplot

comparison heatmap_for_variants.R @ 0:1062d6ad6503 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/snpfreqplot/ commit 1f35303af979c16d9a3126dbc882a59f686ace5d"
author iuc
date Wed, 02 Dec 2020 21:23:06 +0000
parents
children e362b3143cde
comparison
equal deleted inserted replaced
-1:000000000000 0:1062d6ad6503
1 #!/usr/bin/env R
2
3 suppressPackageStartupMessages(library(pheatmap))
4 suppressPackageStartupMessages(library(RColorBrewer))
5 suppressPackageStartupMessages(library(tidyverse))
6
7 fapply <- function(vect_ids, func) {
8 #' List apply but preserve the names
9 res <- lapply(vect_ids, func)
10 names(res) <- vect_ids
11 return(res)
12 }
13
14 # M A I N
15 stopifnot(exists("samples"))
16 variant_files <- fapply(samples$ids, read_and_process) # nolint
17
18 extractall_data <- function(id) {
19 variants <- variant_files[[id]]
20 tmp <- variants %>%
21 mutate(posalt = uni_select) %>%
22 select(posalt, AF)
23 colnames(tmp) <- c("Mutation", id)
24 return(tmp)
25 }
26
27 extractall_annots <- function(id) {
28 variants <- variant_files[[id]]
29 tmp <- variants %>%
30 mutate(posalt = uni_select,
31 effect = EFF....EFFECT, gene = EFF....GENE) %>%
32 select(posalt, effect, gene)
33 return(tmp)
34 }
35
36 # process allele frequencies
37 processed_files <- fapply(samples$ids, extractall_data)
38 final <- as_tibble(
39 processed_files %>%
40 reduce(full_join, by = "Mutation", copy = T))
41
42 final <- final[str_order(final$Mutation, numeric = T), ] %>%
43 column_to_rownames("Mutation") ## sort and set rownames
44 final[final < variant_frequency] <- NA ## adjust the variant frequency:
45 final <- final[rowSums(is.na(final)) != ncol(final), ]
46 final <- t(final)
47 final[is.na(final)] <- 0
48 class(final) <- "numeric"
49
50 # add annotations
51 ## readout annotations
52 processed_annots <- fapply(samples$ids, extractall_annots)
53 ann_final <- processed_annots %>%
54 reduce(function(x, y) {
55 unique(rbind(x, y))}) %>%
56 filter(posalt %in% colnames(final)) ## apply frequency filter
57 ann_final <- as_tibble(ann_final[str_order(
58 ann_final$posalt, numeric = T), ]) %>%
59 column_to_rownames("posalt") ## sort
60
61 # rename annotations
62 trans <- function(x, mapping, replace_missing=NULL) {
63 # helper function for translating effects
64 mapped <- mapping[[x]]
65 if (is.null(mapped)) {
66 if (is.null(replace_missing)) x else replace_missing
67 } else {
68 mapped
69 }
70 }
71
72 # handle translation of classic SnpEff effects to sequence ontology terms
73 # The following list defines the complete mapping between classic and So effect
74 # terms even if not all of these are likely to appear in viral variant data.
75 classic_snpeff_effects_to_so <- list(
76 "coding_sequence_variant", "coding_sequence_variant", "disruptive_inframe_deletion", "disruptive_inframe_insertion", "inframe_deletion", "inframe_insertion", "downstream_gene_variant", "exon_variant", "exon_loss_variant", "frameshift_variant", "gene_variant", "intergenic_variant", "intergenic_region", "conserved_intergenic_variant", "intragenic_variant", "intron_variant", "conserved_intron_variant", "missense_variant", "rare_amino_acid_variant", "splice_acceptor_variant", "splice_donor_variant", "splice_region_variant", "5_prime_UTR_premature_start_codon_variant", "start_lost", "stop_gained", "stop_lost", "synonymous_variant", "start_retained_variant", "stop_retained_variant", "transcript_variant", "upstream_gene_variant", "3_prime_UTR_truncation_+_exon_loss_variant", "3_prime_UTR_variant", "5_prime_UTR_truncation_+_exon_loss_variant", "5_prime_UTR_variant", "initiator_codon_variant", "None", "chromosomal_deletion"
77 )
78 names(classic_snpeff_effects_to_so) <- c(
79 "CDS", "CODON_CHANGE", "CODON_CHANGE_PLUS_CODON_DELETION", "CODON_CHANGE_PLUS_CODON_INSERTION", "CODON_DELETION", "CODON_INSERTION", "DOWNSTREAM", "EXON", "EXON_DELETED", "FRAME_SHIFT", "GENE", "INTERGENIC", "INTERGENIC_REGION", "INTERGENIC_CONSERVED", "INTRAGENIC", "INTRON", "INTRON_CONSERVED", "NON_SYNONYMOUS_CODING", "RARE_AMINO_ACID", "SPLICE_SITE_ACCEPTOR", "SPLICE_SITE_DONOR", "SPLICE_SITE_REGION", "START_GAINED", "START_LOST", "STOP_GAINED", "STOP_LOST", "SYNONYMOUS_CODING", "SYNONYMOUS_START", "SYNONYMOUS_STOP", "TRANSCRIPT", "UPSTREAM", "UTR_3_DELETED", "UTR_3_PRIME", "UTR_5_DELETED", "UTR_5_PRIME", "NON_SYNONYMOUS_START", "NONE", "CHROMOSOME_LARGE_DELETION"
80 )
81 # translate classic effects into SO terms leaving unknown terms (possibly SO already) as is
82 so_effects <- sapply(ann_final$effect, function(x) trans(x, classic_snpeff_effects_to_so))
83
84 # handle further translation of effects we care about
85 so_effects_translation <- list(
86 "non-syn", "syn",
87 "deletion", "deletion", "deletion",
88 "insertion", "insertion", "frame shift",
89 "stop gained", "stop lost"
90 )
91 names(so_effects_translation) <- c(
92 "missense_variant", "synonymous_variant",
93 "disruptive_inframe_deletion", "inframe_deletion", "chromosomal_deletion",
94 "disruptive_inframe_insertion", "inframe_insertion", "frameshift_variant",
95 "stop_gained", "stop_lost"
96 )
97 # translate to our simple terms turning undefined terms into '?'
98 simple_effects <- sapply(so_effects, function(x) trans(x, so_effects_translation, replace_missing = "?"))
99 # complex variant effects (those that do more than one thing) are concatenated
100 # with either '+' (for classic terms) or '&' (for SO terms)
101 simple_effects[grepl("+", so_effects, fixed = TRUE)] <- "complex"
102 simple_effects[grepl("&", so_effects, fixed = TRUE)] <- "complex"
103 simple_effects[so_effects == ""] <- "non-coding"
104
105 ann_final$effect <- simple_effects
106 ann_final$gene <- sub("^$", "NCR", ann_final$gene)
107
108 ## automatically determine gaps for the heatmap
109 gap_vector <- which(!(ann_final$gene[1:length(ann_final$gene) - 1] == # nolint
110 ann_final$gene[2:length(ann_final$gene)]))
111
112 # colormanagement
113 my_colors <- colorRampPalette(c("grey93", "brown", "black")) #heatmap
114 count <- length(unique(ann_final$gene)) #annotations (genes)
115 gene_color <- c(brewer.pal(brewer_color_gene_annotation, n = count))
116 names(gene_color) <- unique(ann_final$gene)
117
118 # colormanagement annotations (effect)
119 ## Define the full set of colors for each effect that we can encounter
120 ## This is not bulletproof. The effect names given here were swapped into the
121 ## data (see above substitutions in ann_final$effect) and so are hard-coded,
122 ## as well as their preferred colors.
123
124 all_colors <- data.frame(
125 color = c("white", "green", "orange", "red",
126 "black", "grey", "yellow", "blue", "purple", "brown"),
127 name = c("non-coding", "syn", "non-syn", "deletion",
128 "frame shift", "stop gained", "stop lost", "insertion",
129 "complex", "?"))
130 ## Reduce the full set to just those that we want
131 detected_effects <- unique(ann_final$effect)
132 subset_colors <- subset(all_colors, name %in% detected_effects)
133 effect_color <- subset_colors$color
134 names(effect_color) <- subset_colors$name
135 color_list <- list(gene_color = gene_color, effect_color = effect_color)
136 names(color_list) <- c("gene", "effect")
137
138 # visualize heatmap
139 if (pheat_number_of_clusters > length(samples$ids)) {
140 print(paste0("[INFO] Number of clusters: User-specified clusters (",
141 pheat_number_of_clusters,
142 ") is greater than the number of samples (",
143 length(samples$ids), ")"))
144 pheat_number_of_clusters <- length(samples$ids)
145 print(paste0("[INFO] Number of clusters: now set to ",
146 pheat_number_of_clusters))
147 }
148
149 get_plot_dims <- function(heat_map) {
150 ## get the dimensions of a pheatmap object
151 ## useful for plot formats that can't be written to a file directly, but
152 ## for which we need to set up a plotting device
153 ## source: https://stackoverflow.com/a/61876386
154 plot_height <- sum(sapply(heat_map$gtable$heights,
155 grid::convertHeight, "in"))
156 plot_width <- sum(sapply(heat_map$gtable$widths,
157 grid::convertWidth, "in"))
158 return(list(height = plot_height, width = plot_width))
159 }
160
161 height <- round(max(c(max(c(
162 16 * (length(unique(ann_final$effect)) +
163 length(unique(ann_final$gene))), 160)) /
164 nrow(final), 15)))
165 width <- round(ratio * height)
166
167
168 if (!(out_ext %in% c("svg", "jpeg", "png", "pdf"))) {
169 stop("Unknown extension: ", ext, ", aborting.")
170 }
171 plot_device <- get(out_ext)
172
173
174 ## A constant scaling factor based on the calculated dimensions
175 ## above does not work for PNG, so we resort to feeding pheatmap
176 ## with a direct filename
177 plot_filename <- NA
178 if (out_ext %in% c("jpeg", "png")) {
179 plot_filename <- out_file
180 }
181
182 ## SVG is not a format pheatmap knows how to write to a file directly.
183 ## As a workaround we
184 ## 1. create the plot object
185 ## 2. get its dimensions
186 ## 3. set up a svg plotting device with these dimensions
187 ## 4. print the heatmap object to the device
188 hm <- pheatmap(
189 final,
190 color = my_colors(100),
191 cellwidth = width,
192 cellheight = height,
193 fontsize_col = round(1 / 3 * width),
194 fontsize_row = round(1 / 3 * min(c(height, width))),
195 clustering_method = pheat_clustering_method,
196 cluster_rows = pheat_clustering,
197 cluster_cols = F,
198 cutree_rows = pheat_number_of_clusters,
199 annotation_col = ann_final,
200 annotation_colors = color_list,
201 filename = plot_filename,
202 gaps_col = gap_vector
203 )
204
205 if (out_ext %in% c("pdf", "svg")) {
206 plot_dims <- get_plot_dims(hm)
207 plot_device(out_file,
208 width = plot_dims$width,
209 height = plot_dims$height)
210 print(hm)
211 dev.off()
212 }