comparison heatmap_colormanipulation/heatmap_extra_v2beta_VERSION.R @ 0:3797463c65f8 draft

Initial upload

0:3797463c65f8

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()