#!/usr/bin/Rscript
library(multtest)
data <- read.table("ex29_hierarchical_clustering.txt", header=T, row.names=1)
data <- log2(data)
dev.new()
plot(density(data[,1]))
for(p in 2:dim(data)[2]) lines(density(data[,p]), col=p)
dev.new()
genes <- sample(1:dim(data)[1], 20)
data.sample <- data[genes,]
plot(x=c(1,20), y=c(min(data.sample),max(data.sample)), col="white", xlab="Patient", ylab="Expression value")
for(i in 1:20)
{
	lines(unlist(data.sample[i,]), col=rainbow(20)[i])
	points(unlist(data.sample[i,]), col=rainbow(20)[i], pch=i)
}
legend(17, max(data.sample), rownames(data.sample), fill=rainbow(20), pch=1:20)

tmp <- mt.teststat(data, rep(0:1, each=8), test="t")
p <- pt(tmp, 15)
dev.new()
plot(density(p))

for(method in c("Bonferroni", "Holm", "FDR"))
{
bon <- mt.rawp2adjp(p, method)
dev.new()
plot(density(bon$adjp))

top20idx <- bon$index[1:20]
top20genes <- rownames(data)[top20idx]
data.top20 <- data[top20idx,]
dev.new()
plot(x=c(1,20), y=c(min(data.top20),max(data.top20)), col="white", xlab="Patient", ylab="Expression value")
for(i in 1:20)
{
	lines(unlist(data.top20[i,]), col=rainbow(20)[i])
	points(unlist(data.top20[i,]), col=rainbow(20)[i], pch=i)
}
legend(17, max(data.top20), rownames(data.top20), fill=rainbow(20), pch=1:20)
}