Mercurial > repos > earlhaminst > plotheatmap

comparison script.R @ 1:05fbc5c891c5 draft default tip

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
planemo upload for repository https://github.com/TGAC/earlham-galaxytools/tree/main/tools/plotheatmap commit 6e40bbe92367310e9d3ec69571d08eb49af7c0a6-dirty
author earlhaminst
date Mon, 24 Feb 2025 11:03:02 +0000
parents bd8fd161908b
children
comparison
equal deleted inserted replaced
0:bd8fd161908b 1:05fbc5c891c5
1 #!/usr/bin/env Rscript 1 #!/usr/bin/env Rscript
2 2
3 args = commandArgs(TRUE) 3 args <- commandArgs(TRUE)
4 countDataPath = args[1] 4 countDataPath <- args[1]
5 statsDataPath = args[2] 5 statsDataPath <- args[2]
6 logFC = args[3] 6 logFC <- args[3]
7 logCPM = args[4] 7 logCPM <- args[4]
8 pValue = args[5] 8 pValue <- args[5]
9 fdr = args[6] 9 fdr <- args[6]
10 10
11 clusterRow = args[7] 11 clusterRow <- args[7]
12 clusterCol = args[8] 12 clusterCol <- args[8]
13 hclustMethod = args[9] 13 hclustMethod <- args[9]
14 14
15 mgColumnNm = as.numeric(args[10]) 15 mgColumnNm <- as.numeric(args[10])
16 mgRowNm = as.numeric(args[11]) 16 mgRowNm <- as.numeric(args[11])
17 17
18 pdfWidth = as.numeric(args[12]) 18 pdfWidth <- as.numeric(args[12])
19 pdfHeight = as.numeric(args[13]) 19 pdfHeight <- as.numeric(args[13])
20 20
21 21
22 if(clusterRow == "Yes"){ 22 if (clusterRow == "Yes") {
23 clusterRow = TRUE 23 clusterRow <- TRUE
24 } else { 24 } else {
25 clusterRow = NA 25 clusterRow <- NA
26 } 26 }
27 27
28 if(clusterCol == "Yes"){ 28 if (clusterCol == "Yes") {
29 clusterCol = TRUE 29 clusterCol <- TRUE
30 } else { 30 } else {
31 clusterCol = NA 31 clusterCol <- NA
32 } 32 }
33 33
34 require(preprocessCore) 34 require(preprocessCore)
35 require(gplots) 35 require(gplots)
36 36
37 #prepare counts data -------------------------------------------------------- 37 # prepare counts data --------------------------------------------------------
38 countData = read.table(countDataPath, comment = "", 38 countData <- read.table(countDataPath,
39 sep = "\t") 39 comment = "",
40 sep = "\t"
41 )
40 42
41 groups = sapply(as.character(countData[1, -1]), strsplit, ":") 43 groups <- sapply(as.character(countData[1, -1]), strsplit, ":")
42 groups = as.vector(t(countData[1, -1])) 44 groups <- as.vector(t(countData[1, -1]))
43 45
44 names = as.vector(t(countData[2, -1])) 46 names <- as.vector(t(countData[2, -1]))
45 47
46 countData = countData[-c(1,2), ] 48 countData <- countData[-c(1, 2), ]
47 rownames(countData) = countData[, 1] 49 rownames(countData) <- countData[, 1]
48 countData = countData[, -1] 50 countData <- countData[, -1]
49 colnames(countData) = names 51 colnames(countData) <- names
50 52
51 countData = countData[, order(groups, names)] 53 countData <- countData[, order(groups, names)]
52 54
53 # prepare stats data ------------------------------------------------------ 55 # prepare stats data ------------------------------------------------------
54 56
55 statsData = read.table(statsDataPath, sep = "\t", 57 statsData <- read.table(statsDataPath,
56 header = T) 58 sep = "\t",
59 header = T
60 )
57 61
58 colnames(statsData)[-1] = sapply(colnames(statsData)[-1], function(x){ 62 colnames(statsData)[-1] <- sapply(colnames(statsData)[-1], function(x) {
59 unlist(strsplit(x, ".", fixed = T))[3] 63 unlist(strsplit(x, ".", fixed = T))[3]
60 }) 64 })
61 65
62 wh = which(abs(statsData$logFC) >= logFC & statsData$logCPM >= logCPM & statsData$PValue <= pValue & statsData$FDR <= fdr) 66 wh <- which(abs(statsData$logFC) >= logFC & statsData$logCPM >= logCPM & statsData$PValue <= pValue & statsData$FDR <= fdr)
63 67
64 for(i in 1:ncol(countData)){ 68 for (i in 1:ncol(countData)) {
65 countData[,i] = as.numeric(countData[,i]) 69 countData[, i] <- as.numeric(countData[, i])
66 } 70 }
67 71
68 countDataNorm = normalize.quantiles(as.matrix(countData), copy = T) 72 countDataNorm <- normalize.quantiles(as.matrix(countData), copy = T)
69 countDataNormLog = log(countDataNorm + 1, 2) 73 countDataNormLog <- log(countDataNorm + 1, 2)
70 74
71 colnames(countDataNormLog) = colnames(countData) 75 colnames(countDataNormLog) <- colnames(countData)
72 rownames(countDataNormLog) = rownames(countData) 76 rownames(countDataNormLog) <- rownames(countData)
73 77
74 #svg("heatmap.svg", width = 3+length(names), height = 1/2*length(wh)) 78 # svg("heatmap.svg", width = 3+length(names), height = 1/2*length(wh))
75 pdf("heatmap.pdf", width = pdfWidth, height = pdfHeight) 79 pdf("heatmap.pdf", width = pdfWidth, height = pdfHeight)
76 80
77 heatmap.2( 81 heatmap.2(
78 countDataNormLog[wh, ], 82 countDataNormLog[wh, ],
79 density.info=c("none"), 83 density.info = c("none"),
80 hclustfun = function(x) hclust(x, method = hclustMethod), 84 hclustfun = function(x) hclust(x, method = hclustMethod),
81 distfun = function(x) as.dist(1-cor(t(x))), 85 distfun = function(x) as.dist(1 - cor(t(x))),
82 col = bluered(50), 86 col = bluered(50),
83 scale = "row", 87 scale = "row",
84 trace = "none", 88 trace = "none",
85 Rowv = clusterRow, 89 Rowv = clusterRow,
86 Colv = clusterCol, 90 Colv = clusterCol,
87 margins = c(mgColumnNm, mgRowNm) 91 margins = c(mgColumnNm, mgRowNm)
88 ) 92 )
89 93
90 dev.off() 94 dev.off()
91