Mercurial > repos > marpiech > rnaseq_pro_workflow_tools

comparison tools/script.R @ 0:c5a812cdf478 draft

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
planemo upload
author marpiech
date Fri, 09 Dec 2016 10:52:35 -0500
parents
children db5d126bf8d0
comparison
equal deleted inserted replaced
-1:000000000000 0:c5a812cdf478
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
11 require(preprocessCore)
12 require(gplots)
13
14 #prepare counts data --------------------------------------------------------
15 countData = read.table(countDataPath, comment = "",
16 sep = "\t")
17
18 groups = sapply(as.character(countData[1, -1]), strsplit, ":")
19 groups = as.vector(t(countData[1, -1]))
20
21 names = as.vector(t(countData[2, -1]))
22
23 countData = countData[-c(1,2), ]
24 rownames(countData) = countData[, 1]
25 countData = countData[, -1]
26 colnames(countData) = names
27
28 countData = countData[, order(groups, names)]
29
30 # prepare stats data ------------------------------------------------------
31
32 statsData = read.table(statsDataPath, sep = "\t",
33 header = T)
34
35 colnames(statsData)[-1] = sapply(colnames(statsData)[-1], function(x){
36 unlist(strsplit(x, ".", fixed = T))[3]
37 })
38
39 wh = which(abs(statsData$logFC) >= logFC & statsData$logCPM >= logCPM & statsData$PValue <= pValue & statsData$FDR <= fdr)
40
41 for(i in 1:ncol(countData)){
42 countData[,i] = as.numeric(countData[,i])
43 }
44
45 countDataNorm = normalize.quantiles(as.matrix(countData), copy = T)
46 countDataNormLog = log(countDataNorm + 1, 2)
47
48 colnames(countDataNormLog) = colnames(countData)
49 rownames(countDataNormLog) = rownames(countData)
50
51 #svg("heatmap.svg", width = 3+length(names), height = 1/2*length(wh))
52 pdf("heatmap.pdf")
53
54 heatmap.2(
55 countDataNormLog[wh, ],
56 density.info=c("none"),
57 hclustfun = function(x) hclust(x, method = "average"),
58 distfun = function(x) as.dist(1-cor(t(x))),
59 col = bluered(50),
60 scale = 'row',
61 trace = "none", #lwid = c(1, length(names)), lhei = c(1,1/3*length(wh)),
62 # Rowv=NA,
63 Colv = NA,
64 margins = c(7, 8)
65 )
66
67 dev.off()
68