-
Notifications
You must be signed in to change notification settings - Fork 0
/
heatmap47.2.R
79 lines (63 loc) · 2.69 KB
/
heatmap47.2.R
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
#!/usr/bin/Rscript
#biocLite("Heatplus")
#biocLite("vegan")
#biocLite("ape")
#makes a heatmap that has rows ordered by upgma of variable
library(Heatplus)
library(vegan)
library(RColorBrewer)
library(gplots)
library(ape)
library(phangorn)
library(seqinr)
args = commandArgs(trailingOnly=TRUE)
all.data <- read.table(args[1],sep="\t",header=TRUE)
row.names(all.data) <- all.data$otu
all.data <- all.data[,-1]
group1 = all.data[1,]
all.data <- all.data[-1,]
head(all.data)
#make data proportions, not abundances, normalise by sample total
#log transfrom if required.
data.prop <- (sweep(all.data, 2, colSums(all.data), FUN="/"))
if (args[4]=="10"){data.s <- log10(data.prop+1)}
if (args[4]=="e"){data.s <- log(data.prop+1)}
if (args[4]=="2"){data.s <- log2(data.prop+1)}
if (args[4]=="a"){data.s <- all.data}
if (args[4]=="p"){data.s <- data.prop}
if (args[4]=="a10"){data.s <- log10(all.data+1)}
if (args[4]=="ae"){data.s <- log(all.data+1)}
if (args[4]=="a2"){data.s <- log2(all.data+1)}
#set group data colours
group1 <- replace(group1,which(group1==1),"black")
group1 <- replace(group1,which(group1==2),"grey") #add a line here to add groups and colours
group1 <- replace(group1,which(group1==3),"blue")
group1 <- replace(group1,which(group1==4),"red")
group1 <- replace(group1,which(group1==5),"yellow")
group1 <- replace(group1,which(group1==6),"orange")
group1 <- replace(group1,which(group1==7),"purple")
group1 <- replace(group1,which(group1==8),"green")
group1 <-t(group1)
s1=args[3]
colscale <- colorRampPalette(c("chartreuse4","yellow", "red"), space = "rgb")(100)
#simple map
#heatmap(as.matrix(data.prop[1:s1,]), Rowv = NA, Colv = NA, col = scaleyellowred, margins = c(10, 2))
#make trees
x <-data.s[1:s1,]
otu.dist <- vegdist(x, method = "bray")
otu.clus <- hclust(otu.dist, "aver")
samples.dist <- vegdist(t(data.prop), method = "bray")
samples.clus <- hclust(samples.dist, "aver")
#plot(as.dendrogram(otu.clus), horiz=TRUE)
pdf(file=args[2],width=20,height=12)
xm = as.numeric(args[5])
ym = as.numeric(args[6])
lwid = c(1,2.8)
#mode(xm)
#mode(ym)
heatmap.2(as.matrix(x), lwid=lwid, Rowv = as.dendrogram(otu.clus), Colv = as.dendrogram(samples.clus), col = colscale,margins=c(xm,ym), ColSideColors = group1,trace=c("none"),srtCol=90,key=FALSE)
dev.off()
#heatmap.2(as.matrix(x),Rowv = as.dendrogram(otu.clus), Colv = as.dendrogram(samples.clus), col = colscale,margins=c(xm,ym), ColSideColors = group1,trace=c("none"),srtCol=0,key=FALSE)
hm<-heatmap.2(as.matrix(x),Rowv = as.dendrogram(otu.clus), Colv = as.dendrogram(samples.clus), col = colscale,margins=c(xm,ym), ColSideColors = group1,trace=c("none"),srtCol=0,key=FALSE)
sorted <- x[match(rev(labels(hm$rowDendrogram)), rownames(x)), ]
write.table(sorted,"table.txt")