view cluster.tools/partition.R @ 8:a58527c632b7 draft

Uploaded

#!/usr/bin/env Rscript

argspec <- c("partition.R help TBD
                \n\n")
args <- commandArgs(TRUE)
if ( length( args ) == 1 && args =="--help") { 
  write(argspec, stderr())
  q();
}

lib.load.quiet <- function( package ) {
   package <- as.character(substitute(package))
   suppressPackageStartupMessages( do.call( "library", list( package=package ) ) )
}
lib.load.quiet(getopt)
lib.load.quiet( amap )
lib.load.quiet( cluster )

##  we're going to use the amap Dist function, but they misname their correlation
##  functions, so re-name them correctly
amap.distance <- c( "euclidean", "maximum", "manhattan", "canberra", "binary",
                    "pearson", "abspearson", "correlation", "abscorrelation", "spearman", "kendall" )
names( amap.distance ) <- c( "euclidean", "maximum", "manhattan", "canberra", "binary",
                             "cosine", "abscosine", "pearson", "abspearson", "spearman", "kendall" )

spec <- matrix( c( "data.fname",         "d", 1, "character",
                   "algorithm",        "a", 2, "character",
                   "distance.metric",    "m", 2, "character", ## must be one supported by R's dist function
                   "dist.obj",           "D", 2, "logical",
                   "direction",          "n", 2, "character",
                   "num.k",              "k", 2, "integer",
                   "output.name",        "o", 2, "character"
                   ),
                nc=4,
                byrow=TRUE
               )

opt <- getopt( spec=spec )
data <- as.matrix( read.delim( opt$data.fname, header=T, row.names=1 , check.names=FALSE ) )

if ( is.null( opt$distance.metric ) ) { opt$distance.metric <- "euclidean" }
if ( is.null( opt$algorithm ) ) { opt$algorithm <- "km" }
if ( is.null( opt$dist.obj ) ) { opt$dist.obj <- FALSE }
if ( is.null( opt$direction ) ) { opt$direction <- "cols"  }
if ( is.null( opt$output.name ) ) { opt$output.name <- "partition.result" }
if ( is.null( opt$num.k ) || ( opt$num.k == -1 )) {
  if ( opt$direction == 'cols' ) {
    opt$num.k <- 5
  } else if ( opt$direction == 'rows' ) {
    opt$num.k <- nrow( data ) / 30  ## we use an estimated average size of gene clusters to be 30
    if ( opt$num.k > 1000 ) {
      opt$num.k <- ( opt$num.k %/% 10 ) * 10
    } else {
      opt$num.k <- ( opt$num.k %/% 5 ) * 5
    }
  }
}


if ( opt$direction == "cols" ) {
  ## need to transpose b/c both kmeans & pam cluster the rows
  ## this shouldn't have an effect upon a distance matrix
  data <- t( data )
}
if ( opt$num.k > nrow( data ) ) {
  err.msg <- paste( "K specified is greater than the number of elements (", opt$direction, ") in data matrix to be clustereed\n", sep="" )
  stop( err.msg )
}

mat.2.b.clustered <- data
if ( opt$dist.obj ) {
  ## To be updated

  mat.2.b.clustered <- as.dist( data )

  if ( opt$algorithm=="km" ) {
    ##clusterAlg is kmeans
    if (verbose) warning()
  }
} else {
  ## this is a data matrix -- we always generate a dist.mat object (b/c we need it
  ##  in case this result is used with a heatmap
  
  ## PAM clustering
  if ( opt$algorithm != "km" ) {

    if ( ! opt$distance.metric %in% names( amap.distance ) ) stop("unsupported distance.")
    mat.2.b.clustered <- Dist( data, method=as.character( amap.distance[ opt$distance.metric ] ) )
    attr( mat.2.b.clustered, "method" ) <- opt$distance.metric

  } else {
    mat.2.b.clustered <- data
  }
}

## now run the clustering
partcl.res <- cl <- NA
if ( opt$algorithm=="pam" ) {
  partcl.res <- pam( x=mat.2.b.clustered,
                     k=opt$num.k,
                     metric=( ifelse( inherits(data, "dist" ), 
                                      attr( data, "method" ),  ## this is ok if data is a dist object (b/c pam will ignore it)
                                      opt$distance.metric) ) ) ##,
                     ##cluster.only=TRUE )
  cl <- partcl.res$clustering

  if ( is.character( partcl.res$medoids ) ) {
    medoids <- data[ partcl.res$medoids, ]
  } else {
    ##partcl.res$medoids is a matrix -- we shouldn't get this (only if mat.2.b.clustered is a data matrix)
    medoids <- partcl.res$medoids
  }
  med.names <- rownames( medoids )
  med.hc <- hclust( as.dist( as.matrix( mat.2.b.clustered )[ med.names, med.names ] ) )
  med.cls <- as.numeric( cl[ med.names[ med.hc$order ] ] )

  cl.list <- lapply( med.cls, function(i) names( cl[ cl %in% i ] ) )
  names( cl.list ) <- med.cls 

  cl.list <- lapply( cl.list,
                     function( elts ) {
                       if ( length( elts ) == 1 ) {
                         retval <- 1
                         names( retval ) <- elts
                       } else {
                         subdist <- as.dist( as.matrix( mat.2.b.clustered )[ elts, elts ] )
                         sub.hc <- hclust( subdist )
                         retval <- sub.hc$order
                         names( retval ) <- sub.hc$labels
                         retval <- sort( retval )
                       }
                       return( retval )
                     }
                    )
  
  fnl.ord <- as.character( unlist( lapply( cl.list, names ) ) )
  cl <- cl[ fnl.ord ]
} else {
  partcl.res <- kmeans( x=mat.2.b.clustered,
                        centers=opt$num.k )
  cl <- partcl.res$cluster
  centroids <- partcl.res$centers
  cent.hc <- hclust( Dist( centroids, method=as.character( amap.distance[ opt$distance.metric ] ) ) )
  cent.cls <- as.numeric( cent.hc$labels[ cent.hc$order ] )

  cl.list <- lapply( cent.cls, function(i) names( cl[ cl %in% i ] ) )
  names( cl.list ) <- cent.cls 

  cl.list <- lapply( cl.list,
                     function( elts ) {
                       if ( length( elts ) == 1 ) {
                         retval <- 1
                         names( retval ) <- elts
                       } else {
                         if ( all( elts %in% colnames( mat.2.b.clustered ) ) ) {
                           submat <- mat.2.b.clustered[ elts, elts ]
                         } else {
                           submat <- mat.2.b.clustered[ elts, ]
                         }
                         subdist <- Dist( submat, method=as.character( amap.distance[ opt$distance.metric ] ) )
                         sub.hc <- hclust( subdist )
                         retval <- sub.hc$order
                         names( retval ) <- sub.hc$labels
                         retval <- sort( retval )
                       }
                       return( retval )
                     }
                    )
  
  fnl.ord <- as.character( unlist( lapply( cl.list, names ) ) )
  cl <- cl[ fnl.ord ]  
}

if ( opt$direction == "cols" ) {
  ## need to transpose back
  data <- t( data )
}
save( file=opt$output.name, partcl.res, cl, data )