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