Wednesday, August 11, 2010

R::Direct Heatmap.R Without Ratiofile

Refer LIXIA folder for Correlation :


Check and remove any row with all "0" and any column with all "0", before Heatmap clustering.

Call this file using CCratio file



library(ClassComparison)
library(ClassDiscovery)



#################### read in data
Genus.ratio <- read.table(filename[x], header=TRUE)

Genus.ratio.data <- Genus.ratio[10: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[10:dim(Genus.ratio)[1],1]
Method <- as.vector(t(Genus.ratio[3, 2:dim(Genus.ratio)[2]]))
Method <- ifelse(Method=="U", "Urine", "Swab")
# HIV <- as.vector(t(Genus.ratio[7, 2:dim(Genus.ratio)[2]]))
# Smoker <- as.vector(t(Genus.ratio[8, 2:dim(Genus.ratio)[2]]))

#################### remove two patients with total count as NA


NAIndex <- which(Genus.ratio[9, 2:dim(Genus.ratio)[2]] == "N/A")

Genus.ratio.data <- Genus.ratio.data[, -NAIndex]
HIV <- HIV[-NAIndex]
Smoker <- Smoker[-NAIndex]
Genus.ratio.dataTotal <- as.numeric(as.vector(t(Genus.ratio[9, 2:dim(Genus.ratio)[2]]))[-NAIndex])






########### 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,])))
}


plotdata <- sweep(plotdata, 2, Genus.ratio.dataTotal, "/")


labRowSymbol <- rownames(plotdata)
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")
col.set1 <- c( "white")



png(heatmapname[x],width=1360, height=3000, res=300)

margins = c(5, 15,11)
colt<-redgreen(75)
col1 <- col.set
keysize <- 3
lmat<-rbind(c(0,5,5), c(0,1,1), c(0,2,2), c(4,3,3), c(0,6, 7))
lwid <- c(keysize, 17,17)
lhei <- c(keysize,1, 1,58, 7)
col2 <- col.set1

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



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

########## color bar NormTumor######################
par(mar = c(0.5, 0, 0, margins[2]))
image(matrix(as.numeric(as.factor(Method))[colInd], ncol=1), col = col1,
axes = FALSE, xaxt='n', yaxt='n', xlab='', ylab='')
mtext(side = 2, "", line = 0.5, cex=0.5, 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=0.8)
axis(4,at=(1:dim(plotdata)[1])+0.3, as.vector(labRowSymbol)[rowInd],las=1 , cex.axis=0.7)



########## 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, 0))
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(Method))$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))

No comments:

Post a Comment