view tools/script.R @ 3:db5d126bf8d0 draft

planemo upload

#!/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()