# 16 Jan 2011, 06h00
#library(ALL)
library(ConsensusClusterPlus)

#source("/projects/remc_bigdata/TCGA/CORE/consensus_clustering/20110610/ConsensusClusterPlus.CORE.xhost09.20110610.R")

##############################################################################
# data
cat("\nread data\n")
#mydata <- read.table(file="/projects/remc_bigdata/TCGA/CORE/miRNA-seq/20110610/expn_matrix_norm_passed_TCGA.221_COAD-READ.isomiR-thresh.q80_278-miRs.20110610.txt", header = TRUE, row.names = 1, sep="\t")
mydata <- read.table(file="/projects/remc_bigdata/TCGA/CORE/miRNA-seq/20110610/expn_matrix_norm_passed_TCGA.221_COAD-READ.isomiR-thresh.q0_1393-miRs.20110610.txt", header = TRUE, row.names = 1, sep="\t")
#/Users/grobertson/GENEREG/ChIP/CHIP_SEQ/SOLEXA/TCGA/COAD-READ/data/miRNA-Seq/20110610/

dim(mydata)
mydata.matrix <- as.matrix(mydata)


##############################################################################
# output path for results
# localhost
cat("\nset directories\n")

output.folder <- "/projects/remc_bigdata/TCGA/CORE/consensus_clustering/20110610"
#output.folder <- "/Users/grobertson/GENEREG/ChIP/CHIP_SEQ/SOLEXA/TCGA/COAD-READ/data/miRNA-Seq/20110610/consensus_clustering"
setwd(output.folder)

#subfolder.title <- "star_km_spearman_q80Var_reps200"
#subfolder.title <- "COAD-READ-221.q80_278-isomiRs.spearman_km_reps1000"
subfolder.title <- "COAD-READ-221.q0_1393-isomiRs.spearman_km_reps1000"

#icl.string <- "CORE.q80.star.20110121.txt"
icl.string <- "CORE-221.q80.km.txt"

output.subfolder <- paste( output.folder,"/", subfolder.title,"/", sep="" )


##############################################################################
kmax=12
cat("\nrun ConsensusClusterPlus()\n")
results = ConsensusClusterPlus( mydata.matrix, maxK=kmax, reps=1000, pItem=0.8, pFeature=1, title=subfolder.title, clusterAlg="km", distance="spearman", seed=1262118388.71279, plot="pdf", writeTable=T, verbose=T )


# Identify an interesting clustering result
#cT <- results[[17]][["consensusTree"]]
# tree order
#cTorder <- cT$order
#write.table( cTorder, file="cT17order.LAML.DNAmeth.hclust_beta_spearman_q50_noNA_reps100.BRCA.20110113.txt",sep="\t")

# a function to take an integer k value and write out an order file (13 Jan 2010)
writeOrder <- function( k, results ){
	cT <-       results[[k]][["consensusTree"]]
	cTorder <-  cT$order
	pathfile <- paste( output.subfolder, "cTorder.k", k, ".txt", sep="" )
	write.table( cTorder, file=pathfile, sep="\t")
}
cat("\nwriting 'item order' files for each clustering result (K value)\n")
for(k in 2:kmax) writeOrder(k,results)

#results[[2]][["consensusTree"]]



##############################################################################
cat("\nrun calcICL()\n")
icl = calcICL( results, title=subfolder.title, plot="pdf", writeTable=T )
icl.cC <- icl[["clusterConsensus"]]
icl.iC <- icl[["itemConsensus"]]

cat("\nwriting calcICL() files\n")
setwd(output.subfolder)
write.table( icl, file=paste( output.subfolder, "icl.", icl.string, sep="" ), sep="\t" )
write.table( icl.cC, file=paste( output.subfolder, "icl.clusterConsensus.", icl.string, sep="" ), sep="\t" )
write.table( icl.iC, file=paste( output.subfolder, "icl.itemConsensus.", icl.string, sep="" ), sep="\t" )