Original::HEATMAP from mail

R codes for producing the heatmap clustering for the item 23 at http://www.microbiota.org/cgi-bin/mum/pilot/analysis.cgi.


The inputs for these R code are BASED on the "genera count" and "phyla count" at http://www.microbiota.org/cgi-bin/mum/pilot/analysis.cgi. Specially, the "genera count" and "phyla count" are the absolute count for each type of bacteria. The R code uses the ratio (i.e., the "absolute count"/"TotalSeq") as the input (might as well write a R code to convert the genera/phyla count data into a genera/phyla ratio data).



library(ClassComparison)
library(ClassDiscovery)


################## read in data

Genus.ratio <- read.table("summary.genus.ratio.txt", header=TRUE)

Genus.ratio.data <- Genus.ratio[5:dim(Genus.ratio)[1], 2:dim(Genus.ratio)[2]]
colnames(Genus.ratio.data) <- colnames(Genus.ratio)[2:dim(Genus.ratio)[2]]
rownames(Genus.ratio.data) <- Genus.ratio[5:dim(Genus.ratio)[1],1]
Infection <- as.vector(t(Genus.ratio[3, 2:dim(Genus.ratio)[2]]))
Infection <- ifelse(Infection=="None", "None", "Yes")

########### hierarchical clustering
plotdata <- matrix(0, dim(Genus.ratio.data)[1], dim(Genus.ratio.data)[2], byrow=TRUE)
rownames(plotdata)<- rownames(Genus.ratio.data)
colnames(plotdata)<- colnames(Genus.ratio.data)
for(i1 in 1:dim(plotdata)[1])
{
plotdata[i1,] <- as.numeric(as.vector(t(Genus.ratio.data[i1,])))
}



labRowSymbol <- rownames(Genus.ratio.data)
sc <- hclust(distanceMatrix((plotdata),"spearman"), "ward")
ddc <- as.dendrogram(sc)
colInd <- order.dendrogram(ddc)

gc <- hclust(distanceMatrix(t(plotdata),"spearman"), "ward")
ddr <- as.dendrogram(gc)
rowInd <- order.dendrogram(ddr)

########### define color
col.set <- c( "purple","yellow")


png(file="Genus-ratio-2.png", width=860, height=1000, res=90)

margins = c(5, 18,13)
colt<-redgreen(75)
col1 <- col.set
keysize <- 1.5
lmat<-rbind(c(0,4,4), c(0,1,1), c(3,2,2), c(0,5, 6))
lwid <- c( keysize, 2,3)
lhei <- c(keysize, 0.4, 16, 1.5)


layout(lmat, widths = lwid, heights = lhei, respect = FALSE)



########## color bar NormTumor######################
par(mar = c(0.5, 0, 0, margins[2]))
image(matrix(as.numeric(as.factor(Infection))[colInd], ncol=1), col = col1,
axes = FALSE, xaxt='n', yaxt='n', xlab='', ylab='')
mtext(side = 2, "Infection", line = 0.5, cex=1.3, las=1)


############# main image ################################
par(mar = c(margins[1],0, 0, margins[2]))

image(1:dim(plotdata)[2], 1:dim(plotdata)[1], t(plotdata[rowInd,colInd]), xlim = 0.5 + c(0, dim(plotdata)[2]), ylim = 0.5 +
c(0, dim(plotdata)[1]), axes = FALSE, xlab = "", ylab = "", col = colt)
axis(1,at=(1:dim(plotdata)[2])+0.3, colnames(plotdata)[colInd],las=2 , cex.axis=2)
axis(4,at=(1:dim(plotdata)[1])+0.3, as.vector(labRowSymbol)[rowInd],las=1 , cex.axis=1.3)



########## color bar DMRest######################
par(mar = c(margins[1], 0, 0, 0))
plot(ddr, horiz = TRUE, axes = FALSE, yaxs = "i", leaflab = "none")

############# dendrogram ##########################
par(mar = c(0, 0, 0, margins[2]))
plot(ddc, axes = FALSE, xaxs = "i", leaflab = "none")


############# dendrogram ##########################


par(mar = c(0, 0, 0, 1.5))
plot(c(-5,5),c(-5,5), axes = FALSE, xaxt = "n", yaxt = "n", main="", xlab="", ylab="", type="n")


rect(-3, -4.5, 2, 4)


legend(-3, 4, attributes(as.factor(Infection))$levels,
col=col1, pch=15 , cex=2, bty="n")

############# color##########################


par(mar=c(2.1, 0.2, 1, 4), cex=1)
dummy.x <- seq(min(plotdata), max(plotdata),
length = length(colt))
dummy.z <- matrix(dummy.x, ncol = 1)
image(x = dummy.x, y = 1, z = dummy.z, xlab = "",
ylab="", yaxt = "n", col = colt )


dev.off()

par(mar=c(5, 4, 4, 2) + 0.1,mfrow=c(1,1))