Mercurial > repos > earlhaminst > plotheatmap
diff script.R @ 0:bd8fd161908b draft default tip
planemo upload for repository https://github.com/TGAC/earlham-galaxytools/tree/master/tools/plotheatmap commit 7e9d6b283879b45d5fe673767668c06104d31609-dirty
author | earlhaminst |
---|---|
date | Wed, 12 Jul 2017 14:29:51 -0400 |
parents | |
children |
line wrap: on
line diff
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/script.R Wed Jul 12 14:29:51 2017 -0400 @@ -0,0 +1,91 @@ +#!/usr/bin/env Rscript + +args = commandArgs(TRUE) +countDataPath = args[1] +statsDataPath = args[2] +logFC = args[3] +logCPM = args[4] +pValue = args[5] +fdr = args[6] + +clusterRow = args[7] +clusterCol = args[8] +hclustMethod = args[9] + +mgColumnNm = as.numeric(args[10]) +mgRowNm = as.numeric(args[11]) + +pdfWidth = as.numeric(args[12]) +pdfHeight = as.numeric(args[13]) + + +if(clusterRow == "Yes"){ + clusterRow = TRUE +} else { + clusterRow = NA +} + +if(clusterCol == "Yes"){ + clusterCol = TRUE +} else { + clusterCol = NA +} + +require(preprocessCore) +require(gplots) + +#prepare counts data -------------------------------------------------------- +countData = read.table(countDataPath, comment = "", + sep = "\t") + +groups = sapply(as.character(countData[1, -1]), strsplit, ":") +groups = as.vector(t(countData[1, -1])) + +names = as.vector(t(countData[2, -1])) + +countData = countData[-c(1,2), ] +rownames(countData) = countData[, 1] +countData = countData[, -1] +colnames(countData) = names + +countData = countData[, order(groups, names)] + +# prepare stats data ------------------------------------------------------ + +statsData = read.table(statsDataPath, sep = "\t", + header = T) + +colnames(statsData)[-1] = sapply(colnames(statsData)[-1], function(x){ + unlist(strsplit(x, ".", fixed = T))[3] +}) + +wh = which(abs(statsData$logFC) >= logFC & statsData$logCPM >= logCPM & statsData$PValue <= pValue & statsData$FDR <= fdr) + +for(i in 1:ncol(countData)){ + countData[,i] = as.numeric(countData[,i]) +} + +countDataNorm = normalize.quantiles(as.matrix(countData), copy = T) +countDataNormLog = log(countDataNorm + 1, 2) + +colnames(countDataNormLog) = colnames(countData) +rownames(countDataNormLog) = rownames(countData) + +#svg("heatmap.svg", width = 3+length(names), height = 1/2*length(wh)) +pdf("heatmap.pdf", width = pdfWidth, height = pdfHeight) + +heatmap.2( + countDataNormLog[wh, ], + density.info=c("none"), + hclustfun = function(x) hclust(x, method = hclustMethod), + distfun = function(x) as.dist(1-cor(t(x))), + col = bluered(50), + scale = "row", + trace = "none", + Rowv = clusterRow, + Colv = clusterCol, + margins = c(mgColumnNm, mgRowNm) +) + +dev.off() +