annotate Dotplot_Release/pheatmap_j.R @ 6:b4bb92033d59 draft default tip

Uploaded

1 lo = function(rown, coln, nrow, ncol, cellheight = NA, cellwidth = NA, treeheight_col, treeheight_row, legend, annotation, annotation_colors, annotation_legend, main, fontsize, fontsize_row, fontsize_col, ...){
2 	# Get height of colnames and length of rownames
3 	if(!is.null(coln[1])){
4 		longest_coln = which.max(strwidth(coln, units = 'in'))
5 		gp = list(fontsize = fontsize_col, ...)
6 		coln_height = unit(1, "grobheight", textGrob(coln[longest_coln], rot = 90, gp = do.call(gpar, gp))) + unit(5, "bigpts")
7 	}
8 	else{
9 		coln_height = unit(5, "bigpts")
10 	}
11 	
12 	if(!is.null(rown[1])){
13 		longest_rown = which.max(strwidth(rown, units = 'in'))
14 		gp = list(fontsize = fontsize_row, ...)
15 		rown_width = unit(1, "grobwidth", textGrob(rown[longest_rown], gp = do.call(gpar, gp))) + unit(10, "bigpts")
16 	}
17 	else{
18 		rown_width = unit(5, "bigpts")
19 	}
20 	
21 	gp = list(fontsize = fontsize, ...)
22 	# Legend position
23 	if(!is.na(legend[1])){
24 		longest_break = which.max(nchar(names(legend)))
25 		longest_break = unit(1.1, "grobwidth", textGrob(as.character(names(legend))[longest_break], gp = do.call(gpar, gp)))
26 		title_length = unit(1.1, "grobwidth", textGrob("Scale", gp = gpar(fontface = "bold", ...)))
27 		legend_width = unit(12, "bigpts") + longest_break * 1.2
28 		legend_width = max(title_length, legend_width)
29 	}
30 	else{
31 		legend_width = unit(0, "bigpts")
32 	}
33 	
34 	# Set main title height
35 	if(is.na(main)){
36 		main_height = unit(0, "npc")
37 	}
38 	else{
39 		main_height = unit(1.5, "grobheight", textGrob(main, gp = gpar(fontsize = 1.3 * fontsize, ...)))
40 	}
41 	
42 	# Column annotations
43 	if(!is.na(annotation[[1]][1])){
44 		# Column annotation height 
45 		annot_height = unit(ncol(annotation) * (8 + 2) + 2, "bigpts")
46 		# Width of the correponding legend
47 		longest_ann = which.max(nchar(as.matrix(annotation)))
48 		annot_legend_width = unit(1.2, "grobwidth", textGrob(as.matrix(annotation)[longest_ann], gp = gpar(...))) + unit(12, "bigpts")
49 		if(!annotation_legend){
50 			annot_legend_width = unit(0, "npc")
51 		}
52 	}
53 	else{
54 		annot_height = unit(0, "bigpts")
55 		annot_legend_width = unit(0, "bigpts")
56 	}
57 	
58 	# Tree height
59 	treeheight_col = unit(treeheight_col, "bigpts") + unit(5, "bigpts")
60 	treeheight_row = unit(treeheight_row, "bigpts") + unit(5, "bigpts") 
61 	
62 	# Set cell sizes
63 	if(is.na(cellwidth)){
64 		matwidth = unit(1, "npc") - rown_width - legend_width - treeheight_row - annot_legend_width
65 	}
66 	else{
67 		matwidth = unit(cellwidth * ncol, "bigpts")
68 	}
69 	
70 	if(is.na(cellheight)){
71 		matheight = unit(1, "npc") - main_height - coln_height - treeheight_col - annot_height
72 	}
73 	else{
74 		matheight = unit(cellheight * nrow, "bigpts")
75 	}	
76 	
77 	
78 	# Produce layout()
79 	pushViewport(viewport(layout = grid.layout(nrow = 5, ncol = 5, widths = unit.c(treeheight_row, matwidth, rown_width, legend_width, annot_legend_width), heights = unit.c(main_height, treeheight_col, annot_height, matheight, coln_height)), gp = do.call(gpar, gp)))
80 	
81 	# Get cell dimensions
82 	pushViewport(vplayout(4, 2))
83 	cellwidth = convertWidth(unit(0:1, "npc"), "bigpts", valueOnly = T)[2] / ncol
84 	cellheight = convertHeight(unit(0:1, "npc"), "bigpts", valueOnly = T)[2] / nrow
85 	upViewport()
86 	
87 	# Return minimal cell dimension in bigpts to decide if borders are drawn
88 	mindim = min(cellwidth, cellheight) 
89 	return(mindim)
90 }
91 
92 draw_dendrogram = function(hc, horizontal = T){
93 	h = hc$height / max(hc$height) / 1.05
94 	m = hc$merge
95 	o = hc$order
96 	n = length(o)
97 
98 	m[m > 0] = n + m[m > 0] 
99 	m[m < 0] = abs(m[m < 0])
00 
01 	dist = matrix(0, nrow = 2 * n - 1, ncol = 2, dimnames = list(NULL, c("x", "y"))) 
02 	dist[1:n, 1] = 1 / n / 2 + (1 / n) * (match(1:n, o) - 1)
03 
04 	for(i in 1:nrow(m)){
05 		dist[n + i, 1] = (dist[m[i, 1], 1] + dist[m[i, 2], 1]) / 2
06 		dist[n + i, 2] = h[i]
07 	}
08 	
09 	draw_connection = function(x1, x2, y1, y2, y){
10 		grid.lines(x = c(x1, x1), y = c(y1, y))
11 		grid.lines(x = c(x2, x2), y = c(y2, y))
12 		grid.lines(x = c(x1, x2), y = c(y, y))
13 	}
14 	
15 	if(horizontal){
16 		for(i in 1:nrow(m)){
17 			draw_connection(dist[m[i, 1], 1], dist[m[i, 2], 1], dist[m[i, 1], 2], dist[m[i, 2], 2], h[i])
18 		}
19 	}
20 	
21 	else{
22 		gr = rectGrob()
23 		pushViewport(viewport(height = unit(1, "grobwidth", gr), width = unit(1, "grobheight", gr), angle = 90))
24 		dist[, 1] = 1 - dist[, 1] 
25 		for(i in 1:nrow(m)){
26 			draw_connection(dist[m[i, 1], 1], dist[m[i, 2], 1], dist[m[i, 1], 2], dist[m[i, 2], 2], h[i])
27 		}
28 		upViewport()
29 	}
30 }
31 
32 draw_matrix = function(matrix, border_color, border_width, fmat, fontsize_number){
33 	n = nrow(matrix)
34 	m = ncol(matrix)
35 	x = (1:m)/m - 1/2/m
36 	y = 1 - ((1:n)/n - 1/2/n)
37 	for(i in 1:m){
38 		grid.rect(x = x[i], y = y[1:n], width = 1/m, height = 1/n, gp = gpar(fill = matrix[,i], col = border_color, lwd = border_width))
39 		if(attr(fmat, "draw")){
40 			grid.text(x = x[i], y = y[1:n], label = fmat[, i], gp = gpar(col = "grey30", fontsize = fontsize_number))
41 		}
42 	}
43 }
44 
45 draw_colnames = function(coln, ...){
46 	m = length(coln)
47 	x = (1:m)/m - 1/2/m
48 	grid.text(coln, x = x, y = unit(0.96, "npc"), just="right", rot = 90, gp = gpar(...))
49 }
50 
51 draw_rownames = function(rown, ...){
52 	n = length(rown)
53 	y = 1 - ((1:n)/n - 1/2/n)
54 	grid.text(rown, x = unit(0.04, "npc"), y = y, vjust = 0.5, hjust = 0, gp = gpar(...))	
55 }
56 
57 draw_legend = function(color, breaks, legend, ...){
58 	height = min(unit(1, "npc"), unit(150, "bigpts"))
59 	pushViewport(viewport(x = 0, y = unit(1, "npc"), just = c(0, 1), height = height))
60 	legend_pos = (legend - min(breaks)) / (max(breaks) - min(breaks))
61 	breaks = (breaks - min(breaks)) / (max(breaks) - min(breaks))
62 	h = breaks[-1] - breaks[-length(breaks)]
63 	grid.rect(x = 0, y = breaks[-length(breaks)], width = unit(10, "bigpts"), height = h, hjust = 0, vjust = 0, gp = gpar(fill = color, col = "#FFFFFF00"))
64 	grid.text(names(legend), x = unit(12, "bigpts"), y = legend_pos, hjust = 0, gp = gpar(...))
65 	upViewport()
66 }
67 
68 convert_annotations = function(annotation, annotation_colors){
69 	new = annotation
70 	for(i in 1:ncol(annotation)){
71 		a = annotation[, i]
72 		b = annotation_colors[[colnames(annotation)[i]]]
73 		if(is.character(a) | is.factor(a)){
74 			a = as.character(a)
75 			if(length(setdiff(a, names(b))) > 0){
76 				stop(sprintf("Factor levels on variable %s do not match with annotation_colors", colnames(annotation)[i]))
77 			}
78 			new[, i] = b[a]
79 		}
80 		else{
81 			a = cut(a, breaks = 100)
82 			new[, i] = colorRampPalette(b)(100)[a]
83 		}
84 	}
85 	return(as.matrix(new))
86 }
87 
88 draw_annotations = function(converted_annotations, border_color, border_width){
89 	n = ncol(converted_annotations)
90 	m = nrow(converted_annotations)
91 	x = (1:m)/m - 1/2/m
92 	y = cumsum(rep(8, n)) - 4 + cumsum(rep(2, n))
93 	for(i in 1:m){
94 		grid.rect(x = x[i], unit(y[1:n], "bigpts"), width = 1/m, height = unit(8, "bigpts"), gp = gpar(fill = converted_annotations[i, ], col = border_color, lwd = border_width))
95 	}
96 }
97 
98 draw_annotation_legend = function(annotation, annotation_colors, border_color, border_width, ...){
99 	y = unit(1, "npc")
00 	text_height = unit(1, "grobheight", textGrob("FGH", gp = gpar(...)))
01 	for(i in names(annotation_colors)){
02 		grid.text(i, x = 0, y = y, vjust = 1, hjust = 0, gp = gpar(fontface = "bold", ...))
03 		y = y - 1.5 * text_height
04 		if(is.character(annotation[, i]) | is.factor(annotation[, i])){
05 			for(j in 1:length(annotation_colors[[i]])){
06 				grid.rect(x = unit(0, "npc"), y = y, hjust = 0, vjust = 1, height = text_height, width = text_height, gp = gpar(col = border_color, lwd = border_width, fill = annotation_colors[[i]][j]))
07 				grid.text(names(annotation_colors[[i]])[j], x = text_height * 1.3, y = y, hjust = 0, vjust = 1, gp = gpar(...))
08 				y = y - 1.5 * text_height
09 			}
10 		}
11 		else{
12 			yy = y - 4 * text_height + seq(0, 1, 0.02) * 4 * text_height
13 			h = 4 * text_height * 0.02
14 			grid.rect(x = unit(0, "npc"), y = yy, hjust = 0, vjust = 1, height = h, width = text_height, gp = gpar(col = "#FFFFFF00", fill = colorRampPalette(annotation_colors[[i]])(50)))
15 			txt = rev(range(grid.pretty(range(annotation[, i], na.rm = TRUE))))
16 			yy = y - c(0, 3) * text_height
17 			grid.text(txt, x = text_height * 1.3, y = yy, hjust = 0, vjust = 1, gp = gpar(...))
18 			y = y - 4.5 * text_height
19 		}
20 		y = y - 1.5 * text_height
21 	}
22 }
23 
24 draw_main = function(text, ...){
25 	grid.text(text, gp = gpar(fontface = "bold", ...))
26 }
27 
28 vplayout = function(x, y){
29 	return(viewport(layout.pos.row = x, layout.pos.col = y))
30 }
31 
32 heatmap_motor = function(matrix, border_color, border_width, cellwidth, cellheight, tree_col, tree_row, treeheight_col, treeheight_row, filename, width, height, breaks, color, legend, annotation, annotation_colors, annotation_legend, main, fontsize, fontsize_row, fontsize_col, fmat, fontsize_number, ...){
33 	grid.newpage()
34 	
35 	# Set layout
36 	mindim = lo(coln = colnames(matrix), rown = rownames(matrix), nrow = nrow(matrix), ncol = ncol(matrix), cellwidth = cellwidth, cellheight = cellheight, treeheight_col = treeheight_col, treeheight_row = treeheight_row, legend = legend, annotation = annotation, annotation_colors = annotation_colors, annotation_legend = annotation_legend, main = main, fontsize = fontsize, fontsize_row = fontsize_row, fontsize_col = fontsize_col,  ...)
37 	
38 	if(!is.na(filename)){
39 		pushViewport(vplayout(1:5, 1:5))
40 		
41 		if(is.na(height)){
42 			height = convertHeight(unit(0:1, "npc"), "inches", valueOnly = T)[2]
43 		}
44 		if(is.na(width)){
45 			width = convertWidth(unit(0:1, "npc"), "inches", valueOnly = T)[2]
46 		}
47 		
48 		# Get file type
49 		r = regexpr("\\.[a-zA-Z]*$", filename)
50 		if(r == -1) stop("Improper filename")
51 		ending = substr(filename, r + 1, r + attr(r, "match.length"))
52 
53 		f = switch(ending,
54 			pdf = function(x, ...) pdf(x, ...),
55 			png = function(x, ...) png(x, units = "in", res = 300, ...),
56 			jpeg = function(x, ...) jpeg(x, units = "in", res = 300, ...),
57 			jpg = function(x, ...) jpeg(x, units = "in", res = 300, ...),
58 			tiff = function(x, ...) tiff(x, units = "in", res = 300, compression = "lzw", ...),
59 			bmp = function(x, ...) bmp(x, units = "in", res = 300, ...),
60 			stop("File type should be: pdf, png, bmp, jpg, tiff")
61 		)
62 		
63 		# print(sprintf("height:%f width:%f", height, width))
64 		f(filename, height = height, width = width)
65 		heatmap_motor(matrix, cellwidth = cellwidth, cellheight = cellheight, border_color = border_color, border_width = border_width, tree_col = tree_col, tree_row = tree_row, treeheight_col = treeheight_col, treeheight_row = treeheight_row, breaks = breaks, color = color, legend = legend, annotation = annotation, annotation_colors = annotation_colors, annotation_legend = annotation_legend, filename = NA, main = main, fontsize = fontsize, fontsize_row = fontsize_row, fontsize_col = fontsize_col, fmat = fmat, fontsize_number =  fontsize_number, ...)
66 		dev.off()
67 		upViewport()
68 		return()
69 	}
70 	
71 	# Omit border color if cell size is too small 
72 	if(mindim < 3) border_color = NA
73 	
74 	# Draw title
75 	if(!is.na(main)){
76 		pushViewport(vplayout(1, 2))
77 		draw_main(main, fontsize = 1.3 * fontsize, ...)
78 		upViewport()
79 	}
80 	
81 	# Draw tree for the columns
82 	if(!is.na(tree_col[[1]][1]) & treeheight_col != 0){
83 		pushViewport(vplayout(2, 2))
84 		draw_dendrogram(tree_col, horizontal = T)
85 		upViewport()
86 	}
87 	
88 	# Draw tree for the rows
89 	if(!is.na(tree_row[[1]][1]) & treeheight_row != 0){
90 		pushViewport(vplayout(4, 1))
91 		draw_dendrogram(tree_row, horizontal = F)
92 		upViewport()
93 	}
94 	
95 	# Draw matrix
96 	pushViewport(vplayout(4, 2))
97 	draw_matrix(matrix, border_color, border_width, fmat, fontsize_number)
98 	upViewport()
99 	
00 	# Draw colnames
01 	if(length(colnames(matrix)) != 0){
02 		pushViewport(vplayout(5, 2))
03 		pars = list(colnames(matrix), fontsize = fontsize_col, ...)
04 		do.call(draw_colnames, pars)
05 		upViewport()
06 	}
07 	
08 	# Draw rownames
09 	if(length(rownames(matrix)) != 0){
10 		pushViewport(vplayout(4, 3))
11 		pars = list(rownames(matrix), fontsize = fontsize_row, ...)
12 		do.call(draw_rownames, pars)
13 		upViewport()
14 	}
15 	
16 	# Draw annotation tracks
17 	if(!is.na(annotation[[1]][1])){
18 		pushViewport(vplayout(3, 2))
19 		converted_annotation = convert_annotations(annotation, annotation_colors)
20 		draw_annotations(converted_annotation, border_color, border_width)
21 		upViewport()
22 	}
23 	
24 	# Draw annotation legend
25 	if(!is.na(annotation[[1]][1]) & annotation_legend){
26 		if(length(rownames(matrix)) != 0){
27 			pushViewport(vplayout(4:5, 5))
28 		}
29 		else{
30 			pushViewport(vplayout(3:5, 5))
31 		}
32 		draw_annotation_legend(annotation, annotation_colors, border_color, border_width, fontsize = fontsize, ...)
33 		upViewport()
34 	}
35 	
36 	# Draw legend
37 	if(!is.na(legend[1])){
38 		length(colnames(matrix))
39 		if(length(rownames(matrix)) != 0){
40 			pushViewport(vplayout(4:5, 4))
41 		}
42 		else{
43 			pushViewport(vplayout(3:5, 4))
44 		}
45 		draw_legend(color, breaks, legend, fontsize = fontsize, ...)
46 		upViewport()
47 	}
48 	
49 	
50 }
51 
52 generate_breaks = function(x, n, center = F){
53 	if(center){
54 		m = max(abs(c(min(x, na.rm = T), max(x, na.rm = T))))
55 		res = seq(-m, m, length.out = n + 1)
56 	}
57 	else{
58 		res = seq(min(x, na.rm = T), max(x, na.rm = T), length.out = n + 1)
59 	}
60 	
61 	return(res)
62 }
63 
64 scale_vec_colours = function(x, col = rainbow(10), breaks = NA){
65 	return(col[as.numeric(cut(x, breaks = breaks, include.lowest = T))])
66 }
67 
68 scale_colours = function(mat, col = rainbow(10), breaks = NA){
69 	mat = as.matrix(mat)
70 	return(matrix(scale_vec_colours(as.vector(mat), col = col, breaks = breaks), nrow(mat), ncol(mat), dimnames = list(rownames(mat), colnames(mat))))
71 }
72 
73 cluster_mat = function(mat, distance, method){
74 	if(!(method %in% c("ward", "single", "complete", "average", "mcquitty", "median", "centroid"))){
75 		stop("clustering method has to one form the list: 'ward', 'single', 'complete', 'average', 'mcquitty', 'median' or 'centroid'.")
76 	}
77 	if(!(distance[1] %in% c("correlation", "euclidean", "maximum", "manhattan", "canberra", "binary", "minkowski")) & class(distance) != "dist"){
78 		print(!(distance[1] %in% c("correlation", "euclidean", "maximum", "manhattan", "canberra", "binary", "minkowski")) | class(distance) != "dist")
79 		stop("distance has to be a dissimilarity structure as produced by dist or one measure  form the list: 'correlation', 'euclidean', 'maximum', 'manhattan', 'canberra', 'binary', 'minkowski'")
80 	}
81 	if(distance[1] == "correlation"){
82 		d = as.dist(1 - cor(t(mat)))
83 	}
84 	else{
85 		if(class(distance) == "dist"){
86 			d = distance
87 		}
88 		else{
89 			d = dist(mat, method = distance)
90 		}
91 	}
92 	
93 	return(hclust(d, method = method))
94 }
95 
96 scale_rows = function(x){
97 	m = apply(x, 1, mean, na.rm = T)
98 	s = apply(x, 1, sd, na.rm = T)
99 	return((x - m) / s)
00 }
01 
02 scale_mat = function(mat, scale){
03 	if(!(scale %in% c("none", "row", "column"))){
04 		stop("scale argument shoud take values: 'none', 'row' or 'column'")
05 	}
06 	mat = switch(scale, none = mat, row = scale_rows(mat), column = t(scale_rows(t(mat))))
07 	return(mat)
08 }
09 
10 generate_annotation_colours = function(annotation, annotation_colors, drop){
11 	if(is.na(annotation_colors)[[1]][1]){
12 		annotation_colors = list()
13 	}
14 	count = 0
15 	for(i in 1:ncol(annotation)){
16 		if(is.character(annotation[, i]) | is.factor(annotation[, i])){
17 			if (is.factor(annotation[, i]) & !drop){
18 				count = count + length(levels(annotation[, i]))
19 			}
20 			else{
21 				count = count + length(unique(annotation[, i]))
22 			}
23 		}
24 	}
25 	
26 	factor_colors = hsv((seq(0, 1, length.out = count + 1)[-1] + 
27       0.2)%%1, 0.7, 0.95)
28 	
29 	set.seed(3453)
30 	
31 	for(i in 1:ncol(annotation)){
32 		if(!(colnames(annotation)[i] %in% names(annotation_colors))){
33 			if(is.character(annotation[, i]) | is.factor(annotation[, i])){
34 				n = length(unique(annotation[, i]))
35 				if (is.factor(annotation[, i]) & !drop){
36 					n = length(levels(annotation[, i]))
37 				}
38 				ind = sample(1:length(factor_colors), n)
39 				annotation_colors[[colnames(annotation)[i]]] = factor_colors[ind]
40 				l = levels(as.factor(annotation[, i]))
41 				l = l[l %in% unique(annotation[, i])]
42 				if (is.factor(annotation[, i]) & !drop){
43 					l = levels(annotation[, i])
44 				}
45 				names(annotation_colors[[colnames(annotation)[i]]]) = l
46 				factor_colors = factor_colors[-ind]
47 			}
48 			else{
49 				r = runif(1)
50 				annotation_colors[[colnames(annotation)[i]]] = hsv(r, c(0.1, 1), 1)
51 			}
52 		}
53 	}
54 	return(annotation_colors)
55 }
56 
57 kmeans_pheatmap = function(mat, k = min(nrow(mat), 150), sd_limit = NA, ...){
58 	# Filter data
59 	if(!is.na(sd_limit)){
60 		s = apply(mat, 1, sd)
61 		mat = mat[s > sd_limit, ]	
62 	}
63 	
64 	# Cluster data
65 	set.seed(1245678)
66 	km = kmeans(mat, k, iter.max = 100)
67 	mat2 = km$centers
68 	
69 	# Compose rownames
70 	t = table(km$cluster)
71 	rownames(mat2) = sprintf("cl%s_size_%d", names(t), t)
72 	
73 	# Draw heatmap
74 	pheatmap(mat2, ...)
75 }
76  
77 #' A function to draw clustered heatmaps.
78 #' 
79 #' A function to draw clustered heatmaps where one has better control over some graphical 
80 #' parameters such as cell size, etc. 
81 #' 
82 #' The function also allows to aggregate the rows using kmeans clustering. This is 
83 #' advisable if number of rows is so big that R cannot handle their hierarchical 
84 #' clustering anymore, roughly more than 1000. Instead of showing all the rows 
85 #' separately one can cluster the rows in advance and show only the cluster centers. 
86 #' The number of clusters can be tuned with parameter kmeans_k.
87 #'
88 #' @param mat numeric matrix of the values to be plotted.
89 #' @param color vector of colors used in heatmap.
90 #' @param kmeans_k the number of kmeans clusters to make, if we want to agggregate the 
91 #' rows before drawing heatmap. If NA then the rows are not aggregated.
92 #' @param breaks a sequence of numbers that covers the range of values in mat and is one 
93 #' element longer than color vector. Used for mapping values to colors. Useful, if needed 
94 #' to map certain values to certain colors, to certain values. If value is NA then the 
95 #' breaks are calculated automatically.
96 #' @param border_color color of cell borders on heatmap, use NA if no border should be 
97 #' drawn.
98 #' @param cellwidth individual cell width in points. If left as NA, then the values 
99 #' depend on the size of plotting window.
00 #' @param cellheight individual cell height in points. If left as NA, 
01 #' then the values depend on the size of plotting window.
02 #' @param scale character indicating if the values should be centered and scaled in 
03 #' either the row direction or the column direction, or none. Corresponding values are 
04 #' \code{"row"}, \code{"column"} and \code{"none"}
05 #' @param cluster_rows boolean values determining if rows should be clustered,
06 #' @param cluster_cols boolean values determining if columns should be clustered.
07 #' @param clustering_distance_rows distance measure used in clustering rows. Possible 
08 #' values are \code{"correlation"} for Pearson correlation and all the distances 
09 #' supported by \code{\link{dist}}, such as \code{"euclidean"}, etc. If the value is none 
10 #' of the above it is assumed that a distance matrix is provided.
11 #' @param clustering_distance_cols distance measure used in clustering columns. Possible 
12 #' values the same as for clustering_distance_rows.
13 #' @param clustering_method clustering method used. Accepts the same values as 
14 #' \code{\link{hclust}}.
15 #' @param treeheight_row the height of a tree for rows, if these are clustered. 
16 #' Default value 50 points.
17 #' @param treeheight_col the height of a tree for columns, if these are clustered. 
18 #' Default value 50 points.
19 #' @param legend logical to determine if legend should be drawn or not.
20 #' @param legend_breaks vector of breakpoints for the legend.
21 #' @param legend_labels vector of labels for the \code{legend_breaks}.
22 #' @param annotation data frame that specifies the annotations shown on top of the 
23 #' columns. Each row defines the features for a specific column. The columns in the data 
24 #' and rows in the annotation are matched using corresponding row and column names. Note 
25 #' that color schemes takes into account if variable is continuous or discrete.
26 #' @param annotation_colors list for specifying annotation track colors manually. It is 
27 #' possible to define the colors for only some of the features. Check examples for 
28 #' details.
29 #' @param annotation_legend boolean value showing if the legend for annotation tracks 
30 #' should be drawn. 
31 #' @param drop_levels logical to determine if unused levels are also shown in the legend
32 #' @param show_rownames boolean specifying if column names are be shown.
33 #' @param show_colnames boolean specifying if column names are be shown.
34 #' @param main the title of the plot
35 #' @param fontsize base fontsize for the plot 
36 #' @param fontsize_row fontsize for rownames (Default: fontsize) 
37 #' @param fontsize_col fontsize for colnames (Default: fontsize) 
38 #' @param display_numbers logical determining if the numeric values are also printed to 
39 #' the cells. 
40 #' @param number_format format strings (C printf style) of the numbers shown in cells. 
41 #' For example "\code{\%.2f}" shows 2 decimal places and "\code{\%.1e}" shows exponential 
42 #' notation (see more in \code{\link{sprintf}}).    
43 #' @param fontsize_number fontsize of the numbers displayed in cells
44 #' @param filename file path where to save the picture. Filetype is decided by 
45 #' the extension in the path. Currently following formats are supported: png, pdf, tiff,
46 #'  bmp, jpeg. Even if the plot does not fit into the plotting window, the file size is 
47 #' calculated so that the plot would fit there, unless specified otherwise.
48 #' @param width manual option for determining the output file width in inches.
49 #' @param height manual option for determining the output file height in inches.
50 #' @param \dots graphical parameters for the text used in plot. Parameters passed to 
51 #' \code{\link{grid.text}}, see \code{\link{gpar}}. 
52 #' 
53 #' @return 
54 #' Invisibly a list of components 
55 #' \itemize{
56 #' 	\item \code{tree_row} the clustering of rows as \code{\link{hclust}} object 
57 #' 	\item \code{tree_col} the clustering of columns as \code{\link{hclust}} object
58 #' 	\item \code{kmeans} the kmeans clustering of rows if parameter \code{kmeans_k} was 
59 #' specified 
60 #' }
61 #' 
62 #' @author  Raivo Kolde <rkolde@@gmail.com>
63 #' @examples
64 #'  # Generate some data
65 #' test = matrix(rnorm(200), 20, 10)
66 #' test[1:10, seq(1, 10, 2)] = test[1:10, seq(1, 10, 2)] + 3
67 #' test[11:20, seq(2, 10, 2)] = test[11:20, seq(2, 10, 2)] + 2
68 #' test[15:20, seq(2, 10, 2)] = test[15:20, seq(2, 10, 2)] + 4
69 #' colnames(test) = paste("Test", 1:10, sep = "")
70 #' rownames(test) = paste("Gene", 1:20, sep = "")
71 #' 
72 #' # Draw heatmaps
73 #' pheatmap(test)
74 #' pheatmap(test, kmeans_k = 2)
75 #' pheatmap(test, scale = "row", clustering_distance_rows = "correlation")
76 #' pheatmap(test, color = colorRampPalette(c("navy", "white", "firebrick3"))(50))
77 #' pheatmap(test, cluster_row = FALSE)
78 #' pheatmap(test, legend = FALSE)
79 #' pheatmap(test, display_numbers = TRUE)
80 #' pheatmap(test, display_numbers = TRUE, number_format = "%.1e")
81 #' pheatmap(test, cluster_row = FALSE, legend_breaks = -1:4, legend_labels = c("0", 
82 #' "1e-4", "1e-3", "1e-2", "1e-1", "1"))
83 #' pheatmap(test, cellwidth = 15, cellheight = 12, main = "Example heatmap")
84 #' pheatmap(test, cellwidth = 15, cellheight = 12, fontsize = 8, filename = "test.pdf")
85 #' 
86 #' 
87 #' # Generate column annotations
88 #' annotation = data.frame(Var1 = factor(1:10 %% 2 == 0, 
89 #' 				labels = c("Class1", "Class2")), Var2 = 1:10)
90 #' annotation$Var1 = factor(annotation$Var1, levels = c("Class1", "Class2", "Class3"))
91 #' rownames(annotation) = paste("Test", 1:10, sep = "")
92 #' 
93 #' pheatmap(test, annotation = annotation)
94 #' pheatmap(test, annotation = annotation, annotation_legend = FALSE)
95 #' pheatmap(test, annotation = annotation, annotation_legend = FALSE, drop_levels = FALSE)
96 #' 
97 #' # Specify colors
98 #' Var1 = c("navy", "darkgreen")
99 #' names(Var1) = c("Class1", "Class2")
00 #' Var2 = c("lightgreen", "navy")
01 #' 
02 #' ann_colors = list(Var1 = Var1, Var2 = Var2)
03 #' 
04 #' pheatmap(test, annotation = annotation, annotation_colors = ann_colors, main = "Example")
05 #' 
06 #' # Specifying clustering from distance matrix
07 #' drows = dist(test, method = "minkowski")
08 #' dcols = dist(t(test), method = "minkowski")
09 #' pheatmap(test, clustering_distance_rows = drows, clustering_distance_cols = dcols)
10 #'
11 #' @export
12 pheatmap_j = function(mat, color = colorRampPalette(rev(brewer.pal(n = 7, name = "RdYlBu")))(100), kmeans_k = NA, breaks = NA, border_color = "grey60", border_width = 1, cellwidth = NA, cellheight = NA, scale = "none", cluster_rows = TRUE, cluster_cols = TRUE, clustering_distance_rows = "euclidean", clustering_distance_cols = "euclidean", clustering_method = "complete",  treeheight_row = ifelse(cluster_rows, 50, 0), treeheight_col = ifelse(cluster_cols, 50, 0), legend = TRUE, legend_breaks = NA, legend_labels = NA, annotation = NA, annotation_colors = NA, annotation_legend = TRUE, drop_levels = TRUE, show_rownames = T, show_colnames = T, main = NA, fontsize = 10, fontsize_row = fontsize, fontsize_col = fontsize, display_numbers = F, number_format = "%.2f", fontsize_number = 0.8 * fontsize, filename = NA, width = NA, height = NA, ...){
13 	
14 	# Preprocess matrix
15 	mat = as.matrix(mat)
16 	if(scale != "none"){
17 		mat = scale_mat(mat, scale)
18 		if(is.na(breaks)){
19 			breaks = generate_breaks(mat, length(color), center = T)
20 		}
21 	}
22 	
23 	
24 	# Kmeans
25 	if(!is.na(kmeans_k)){
26 		# Cluster data
27 		km = kmeans(mat, kmeans_k, iter.max = 100)
28 		mat = km$centers
29 
30 		# Compose rownames
31 		t = table(km$cluster)
32 		rownames(mat) = sprintf("cl%s_size_%d", names(t), t)
33 	}
34 	else{
35 		km = NA
36 	}
37 	
38 	# Do clustering
39 	if(cluster_rows){
40 		tree_row = cluster_mat(mat, distance = clustering_distance_rows, method = clustering_method)
41 		mat = mat[tree_row$order, , drop = FALSE]
42 	}
43 	else{
44 		tree_row = NA
45 		treeheight_row = 0
46 	}
47 	
48 	if(cluster_cols){
49 		tree_col = cluster_mat(t(mat), distance = clustering_distance_cols, method = clustering_method)
50 		mat = mat[, tree_col$order, drop = FALSE]
51 	}
52 	else{
53 		tree_col = NA
54 		treeheight_col = 0
55 	}
56 	
57 	# Format numbers to be displayed in cells 
58 	if(display_numbers){
59 		fmat = matrix(sprintf(number_format, mat), nrow = nrow(mat), ncol = ncol(mat))
60 		attr(fmat, "draw") = TRUE
61 	}
62 	else{
63 		fmat = matrix(NA, nrow = nrow(mat), ncol = ncol(mat))
64 		attr(fmat, "draw") = FALSE
65 	}
66 	
67 	
68 	# Colors and scales
69 	if(!is.na(legend_breaks[1]) & !is.na(legend_labels[1])){
70 		if(length(legend_breaks) != length(legend_labels)){
71 			stop("Lengths of legend_breaks and legend_labels must be the same")
72 		}
73 	}
74 	
75 	
76 	if(is.na(breaks[1])){
77       breaks = generate_breaks(as.vector(mat), length(color))
78   }
79   if (legend & is.na(legend_breaks[1])) {
80       legend = grid.pretty(range(as.vector(breaks)))
81 			names(legend) = legend
82   }
83 	else if(legend & !is.na(legend_breaks[1])){
84 		legend = legend_breaks[legend_breaks >= min(breaks) & legend_breaks <= max(breaks)]
85 		
86 		if(!is.na(legend_labels[1])){
87 			legend_labels = legend_labels[legend_breaks >= min(breaks) & legend_breaks <= max(breaks)]
88 			names(legend) = legend_labels
89 		}
90 		else{
91 			names(legend) = legend
92 		}
93 	}
94   else {
95       legend = NA
96   }
97 	mat = scale_colours(mat, col = color, breaks = breaks)
98 	
99 	# Preparing annotation colors
00 	if(!is.na(annotation[[1]][1])){
01 		annotation = annotation[colnames(mat), , drop = F]
02 		annotation_colors = generate_annotation_colours(annotation, annotation_colors, drop = drop_levels)
03 	}
04 	
05 	if(!show_rownames){
06 		rownames(mat) = NULL
07 	}
08 	
09 	if(!show_colnames){
10 		colnames(mat) = NULL
11 	}
12 	
13 	# Draw heatmap
14 	heatmap_motor(mat, border_color = border_color, border_width = border_width, cellwidth = cellwidth, cellheight = cellheight, treeheight_col = treeheight_col, treeheight_row = treeheight_row, tree_col = tree_col, tree_row = tree_row, filename = filename, width = width, height = height, breaks = breaks, color = color, legend = legend, annotation = annotation, annotation_colors = annotation_colors, annotation_legend = annotation_legend, main = main, fontsize = fontsize, fontsize_row = fontsize_row, fontsize_col = fontsize_col, fmat = fmat, fontsize_number = fontsize_number, ...)
15 	
16 	invisible(list(tree_row = tree_row, tree_col = tree_col, kmeans = km))
17 }
18 
19