-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathheatmap_plain.R
74 lines (55 loc) · 2.08 KB
/
heatmap_plain.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
#!/usr/bin/Rscript
#source("https://bioconductor.org/biocLite.R")
#biocLite("Heatplus")
#biocLite("vegan")
#biocLite("ape")
#install.packages("phangorn")
#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,check.names=FALSE)
row.names(all.data) <- all.data$otu
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),"red")
#group1 <- replace(group1,which(group1==2),"blue")
#group1 <-t(group1)
#cbind(names(data.prop), 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.prop[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, leaflab = "none")
pdf(file=args[2],width=20,height=20)
xm = 20 #as.numeric(args[5])
ym = 20 #as.numeric(args[6])
pdf(file=args[2],width=as.numeric(args[5]),height=as.numeric(args[6]))
font_size=as.numeric(args[7])
#mode(xm)
#mode(ym)
heatmap.2(as.matrix(x),Rowv = as.dendrogram(otu.clus), Colv = as.dendrogram(samples.clus), col = colscale,margins=c(xm,ym),trace=c("none"),srtCol=90,key=FALSE, cexRow=font_size, cexCol=font_size)
dev.off()