|
0
|
1 #!/usr/bin/env Rscript
|
|
|
2
|
|
|
3 argspec <- c("hclust.R help TBD
|
|
|
4 \n\n")
|
|
|
5 args <- commandArgs(TRUE)
|
|
|
6 if ( length( args ) == 1 && args =="--help") {
|
|
|
7 write(argspec, stderr())
|
|
|
8 q();
|
|
|
9 }
|
|
|
10
|
|
|
11 lib.load.quiet <- function( package ) {
|
|
|
12 package <- as.character(substitute(package))
|
|
|
13 suppressPackageStartupMessages( do.call( "library", list( package=package ) ) )
|
|
|
14 }
|
|
|
15 lib.load.quiet(getopt)
|
|
|
16 lib.load.quiet( amap )
|
|
|
17 ## if any of the faster clustering methods are available on this system, load them
|
|
|
18 if ( any( c( 'flashClust', 'fastcluster' ) %in% installed.packages() ) ) {
|
|
|
19 if ( 'flashClust' %in% installed.packages() ) {
|
|
|
20 lib.load.quiet( flashClust )
|
|
|
21 } else {
|
|
|
22 if ( 'fastcluster' %in% installed.packages() ) {
|
|
|
23 lib.load.quiet( fastcluster )
|
|
|
24 }
|
|
|
25 }
|
|
|
26 }
|
|
|
27
|
|
|
28 spec <- matrix( c( "data.fname", "d", 1, "character",
|
|
|
29 "num.k", "k", 1, "integer",
|
|
|
30 "distance.metric", "m", 2, "character",
|
|
|
31 "dist.obj", "D", 2, "logical",
|
|
|
32 "direction", "n", 2, "character",
|
|
|
33 "linkage", "l", 2, "character",
|
|
|
34 "output.name", "o", 2, "character"
|
|
|
35 ),
|
|
|
36 nc=4,
|
|
|
37 byrow=TRUE
|
|
|
38 )
|
|
|
39
|
|
|
40 opt <- getopt( spec=spec )
|
|
|
41
|
|
8
|
42 data <- as.matrix( read.delim( opt$data.fname, header=T, row.names=1 , check.names=FALSE ) )
|
|
0
|
43 if ( is.null( opt$distance.metric ) ) { opt$distance.metric <- "euclidean" }
|
|
|
44 if ( is.null( opt$dist.obj ) ) { opt$dist.obj <- FALSE }
|
|
|
45 if ( is.null( opt$direction ) ) { opt$direction <- "cols" }
|
|
|
46 if ( is.null( opt$linkage ) ) { opt$linkage <- "average" }
|
|
|
47 if ( is.null( opt$output.name ) ) { opt$output.name <- "hclust.result.rda" }
|
|
8
|
48 if ( is.null( opt$num.k ) || ( opt$num.k == -1 )) {
|
|
|
49 if ( opt$direction == 'cols' ) {
|
|
|
50 opt$num.k <- 5
|
|
|
51 } else if ( opt$direction == 'rows' ) {
|
|
|
52 opt$num.k <- nrow( data ) / 30 ## we use an estimated average size of gene clusters to be 30
|
|
|
53 if ( opt$num.k > 1000 ) {
|
|
|
54 opt$num.k <- ( opt$num.k %/% 10 ) * 10
|
|
|
55 } else {
|
|
|
56 opt$num.k <- ( opt$num.k %/% 5 ) * 5
|
|
|
57 }
|
|
|
58 }
|
|
|
59 }
|
|
0
|
60
|
|
|
61 if ( opt$direction == "cols" ) {
|
|
|
62 ## need to transpose b/c both kmeans & pam cluster the rows
|
|
|
63 ## this shouldn't have an effect upon a distance matrix
|
|
|
64 data <- t( data )
|
|
|
65 }
|
|
|
66 if ( opt$num.k > nrow( data ) ) {
|
|
|
67 err.msg <- paste( "K specified is greater than the number of elements (", opt$direction, ") in data matrix to be clustereed\n", sep="" )
|
|
|
68 stop( err.msg )
|
|
|
69 }
|
|
|
70
|
|
|
71 if ( opt$dist.obj ) {
|
|
|
72 dist.mat <- as.dist( data )
|
|
|
73 } else {
|
|
|
74 ## we're going to use the amap Dist function, but they misname their correlation
|
|
|
75 ## functions, so re-name them correctly
|
|
|
76 amap.distance <- c( "euclidean", "maximum", "manhattan", "canberra", "binary",
|
|
|
77 "pearson", "abspearson", "correlation", "abscorrelation", "spearman", "kendall" )
|
|
|
78 names( amap.distance ) <- c( "euclidean", "maximum", "manhattan", "canberra", "binary",
|
|
|
79 "cosine", "abscosine", "pearson", "abspearson", "spearman", "kendall" )
|
|
|
80
|
|
|
81 if ( ! opt$distance.metric %in% names( amap.distance ) ) stop("unsupported distance.")
|
|
|
82 dist.mat <- Dist( data, method=as.character( amap.distance[ opt$distance.metric ] ) )
|
|
|
83 attr( dist.mat, "method" ) <- opt$distance.metric
|
|
|
84 }
|
|
|
85
|
|
|
86 ## now, do the clustering
|
|
|
87 treecl.res <- hclust( dist.mat, method=opt$linkage )
|
|
|
88 cutree.res <- cutree( treecl.res, k=opt$num.k )
|
|
|
89 ##cl <- cbind( names( cutree.res ), as.numeric( cutree.res ) )
|
|
|
90 ##colnames( cl ) <- c( "ID", "class" )
|
|
|
91
|
|
|
92 if ( opt$direction == "cols" ) {
|
|
|
93 ## need to re-transpose the data back to it's original dimensionality
|
|
|
94 data <- t( data )
|
|
|
95 }
|
|
|
96
|
|
|
97 cl <- cutree.res
|
|
|
98 save( file=opt$output.name, treecl.res, cl, data )
|
