Mercurial > repos > mir-bioinf > heatmap_colormanipulation
view heatmap_colormanipulation/heatmap_extra_v2beta_VERSION.R @ 0:3797463c65f8 draft
Initial upload
author | mir-bioinf |
---|---|
date | Mon, 20 Apr 2015 15:26:53 -0400 |
parents | |
children | aa592c2c26cb |
line wrap: on
line source
sink(file="/tmp/none") sink("/dev/null") options(warn=-1) options(echo=F) args <- commandArgs(trailingOnly = T) #title <- args[17] Rowcorr <- args[2] Rowlink <- args[3] Colcorr <- args[4] Collink <- args[5] #Xlab <- args[18] #Ylab <- args[19] inputfile <- args[1] Var_cols <- args[6] Scale_var <- args[7] Remove_na <- args[8] header_yes <- args[9] rowhead_yes <- args[10] color_grad <- args[11] color_min <- args[12] #max_val_mincol <- args[13] color_max <- args[13] #min_val_maxcol <- args[15] out_file <- args[14] #logfile <- args[15] ##Now for title, xlabel, and ylabel (spaces are hard to deal with here): #title <- args[17] #Xlab <- args[18] #Ylab <- args[19] stoptime = 0 argIndex = 16 everything = args[argIndex] debugcounter = 0 suppressMessages(library(gplots)) Rinfo = sessionInfo() Rinfo_pkg = sessionInfo(package="gplots") gplots_info = Rinfo_pkg$otherPkgs #sink(logfile) Rinfo gplots_info sink(file="/tmp/none") sink("/dev/null") #cat(paste("arg value is ",args[argIndex],".\n"),file=logfile,append="TRUE") #while (stoptime < 1){ while ((stoptime < 1)&&(debugcounter<50)){ argIndex=argIndex+1 # cat(paste("in while loop now, arg index is ",argIndex,".\n"),file=logfile,append="TRUE") # cat(paste("arg value is ",args[argIndex],".\n"),file=logfile,append="TRUE") everything = paste(everything,args[argIndex]) if (args[argIndex]=="ZZZZ_END") { stoptime = 1 } debugcounter=debugcounter+1 } argIndex=argIndex+1 #cat(paste("Out of while loop. arg index value is now ",argIndex,".\n"),file=logfile,append="TRUE") splitThese = strsplit(everything,"[@]") title = splitThese[[1]][1] Xlab = splitThese[[1]][2] Ylab = splitThese[[1]][3] ##Now grab the rest of the arguments passed in: colorManip = args[argIndex] argIndex = argIndex+1 #cat(paste("Color manip value is ",colorManip,".\n"),file=logfile,append="TRUE") if ((colorManip == "InnerClip") || (colorManip == "OuterClip")) { LowClipVal = as.numeric(args[argIndex]) HighClipVal = as.numeric(args[argIndex+1]) # cat(paste("Two vals to clip: ",LowClipVal,HighClipVal,".\n"),file=logfile,append="TRUE") } else { ClipVal = as.numeric(args[argIndex]) # cat(paste("One val to clip: ",ClipVal,".\n"),file=logfile,append="TRUE") } if (header_yes == "yes") { inp = read.table(inputfile,stringsAsFactors=F, header=T, sep="\t") } else { inp = read.table(inputfile,stringsAsFactors=F, sep="\t") } these_cols = read.csv(text=Var_cols,header=F) if (ncol(these_cols)<2) { x = data.frame(cbind(inp[, c(as.matrix(these_cols))],inp[, c(as.matrix(these_cols))])) currentColNames=colnames(x) labColVar = c(currentColNames[1],"") } else { x = inp[, c(as.matrix(these_cols))] labColVar = colnames(x) } genemat = do.call(cbind,x) x = apply(genemat,2,as.numeric) scale_value = Scale_var na_rm_value = FALSE if (Remove_na == "yes") { na_rm_value = TRUE } if (rowhead_yes == "yes") { rownames(x)=inp[[1]] } pdf(out_file) if ((Rowcorr=="none") && (Colcorr!="none")) { dendro_val = "column" } else if ((Rowcorr!="none") && (Colcorr=="none")) { dendro_val = "row" } if ((Rowcorr=="none") && (Colcorr=="none")) { dendro_val = "none" } if ((Rowcorr!="none") && (Colcorr!="none")) { dendro_val = "both" } if (Rowcorr == "none") { Rowv_val = FALSE } else { Rcor = cor(t(x),method=Rowcorr) R_clust = hclust(as.dist(1-Rcor),method=Rowlink) R_dendro = as.dendrogram(R_clust) Rowv_val = R_dendro } ##Column clustering (if any) set up: if (Colcorr == "none") { Colv_val = FALSE } else { Ccor = cor(x,method=Colcorr) C_clust = hclust(as.dist(1-Ccor),method=Collink) C_dendro = as.dendrogram(C_clust) Colv_val = C_dendro } par(cex.main=0.8) ##font size for title ##Estimate good guesses for font sizes of rows and columns: font_r1 = 0.2 + 1/log10(nrow(x)) ##default done in heatmap, based on number of rows font_size_r = min(0.8,font_r1) font_c1 = 0.2 + 1/log10(ncol(x)) ##default done in heatmap, based on number of columns font_size_c = min(0.8,font_c1) #min_value = min(x) ##x should be the original data matrix #max_value = max(x) if (colorManip == "InnerClip") { min_value = min(x) max_value = max(x) max_val_mincol = LowClipVal min_val_maxcol = HighClipVal } else if (colorManip == "OuterClip") { min_value = LowClipVal max_value = HighClipVal ##How do we set the other values if 0 isn't included in the range? Probably want central color to be center value: if ((min_value<=0)&&(max_value>=0)) { max_val_mincol = 0 ##will be reset later to account for slight tolerance (so black is included) min_val_maxcol = 0 } else { ##0 is not in range, center around halfway point max_val_mincol = (min_value+max_value)/2 + 0.00005 min_val_maxcol = (min_value+max_value)/2 - 0.00005 } } else if (colorManip == "ClipMax") { min_value = min(x) max_value = ClipVal if ((min_value<=0)&&(max_value>=0)) { max_val_mincol = 0 ##will be reset later to account for slight tolerance (so black is included) min_val_maxcol = 0 } else { ##0 is not in range, center around halfway point max_val_mincol = (min_value+max_value)/2 + 0.00005 min_val_maxcol = (min_value+max_value)/2 - 0.00005 } } else { min_value = ClipVal max_value = max(x) if ((min_value<=0)&&(max_value>=0)) { max_val_mincol = 0 ##will be reset later to account for slight tolerance (so black is included) min_val_maxcol = 0 } else { ##0 is not in range, center around halfway point max_val_mincol = (min_value+max_value)/2 + 0.00005 min_val_maxcol = (min_value+max_value)/2 - 0.00005 } } ##is 0 included in the data range? if so we want it centered if ((min_value <= 0) && (max_value >=0)) { sym_breaks_value = "TRUE" sym_key_value = "TRUE" } else { sym_breaks_value = "FALSE" sym_key_value = "FALSE" } if (color_grad == "double") { my_palette = colorRampPalette(c(color_min,"black",color_max))(n=299) } else { my_palette = colorRampPalette(c(color_min,color_max))(n=299) } ##Need some tolerance otherwise black won't be included for 0 if ((max_val_mincol==0) && (min_val_maxcol==0)) { max_val_mincol = -0.00005 min_val_maxcol = 0.00005 } #cat(paste("max_val_mincol value is ",max_val_mincol,".\n"),file=logfile,append="TRUE") #cat(paste("min_val_maxcol value is ",min_val_maxcol,".\n"),file=logfile,append="TRUE") colors = c(seq(min_value,max_val_mincol,length=100),seq(max_val_mincol,min_val_maxcol,length=100),seq(min_val_maxcol,max_value,length=100)) ##Call heatmap.2. cexCol value is constant to account for long sample names heatmap.2(x, margins=c(9,10), main=title, xlab=Xlab, ylab=Ylab, cexCol=font_size_c, cexRow=font_size_r, scale=scale_value, symbreaks=sym_breaks_value, symm=F, symkey=sym_key_value, na.rm=na_rm_value, trace="none", col=my_palette, breaks=colors, dendrogram=dendro_val, Rowv=Rowv_val, Colv=Colv_val, labCol=labColVar) ## Close the PDF file devname = dev.off()