Mercurial > repos > peter-waltman > ucsc_cluster_tools
comparison cluster.tools/cluster.2.centroid.R @ 2:b442996b66ae draft
Uploaded
author | peter-waltman |
---|---|
date | Wed, 27 Feb 2013 20:17:04 -0500 |
parents | |
children |
comparison
equal
deleted
inserted
replaced
1:e25d2bece0a2 | 2:b442996b66ae |
---|---|
1 #!/usr/bin/env Rscript | |
2 argspec <- c("tab.2.cdt.R converts a data matrix to cdt format | |
3 | |
4 Usage: | |
5 | |
6 Optional: | |
7 | |
8 \n\n") | |
9 args <- commandArgs(TRUE) | |
10 if ( length( args ) == 1 && args =="--help") { | |
11 write(argspec, stderr()) | |
12 q(); | |
13 } | |
14 | |
15 lib.load.quiet <- function( package ) { | |
16 package <- as.character(substitute(package)) | |
17 suppressPackageStartupMessages( do.call( "library", list( package=package ) ) ) | |
18 } | |
19 | |
20 lib.load.quiet( getopt ) | |
21 lib.load.quiet( amap ) | |
22 if ( any( c( 'flashClust', 'fastcluster' ) %in% installed.packages() ) ) { | |
23 if ( 'flashClust' %in% installed.packages() ) { | |
24 lib.load.quiet( flashClust ) | |
25 } else { | |
26 if ( 'fastcluster' %in% installed.packages() ) { | |
27 lib.load.quiet( fastcluster ) | |
28 } | |
29 } | |
30 } | |
31 | |
32 spec <- matrix( c( "dataset", "d", 1, "character", | |
33 "gen.new.dgram", "g", 2, "character", | |
34 "output.fname", "o", 2, "character" | |
35 ), | |
36 nc=4, | |
37 byrow=TRUE | |
38 ) | |
39 | |
40 | |
41 opt <- getopt( spec=spec ) | |
42 if ( is.null( opt$output.report.dir ) ) { opt$output.report.dir <- "report" } | |
43 if ( is.null( opt$output.fname ) ) { opt$output.fname <- file.path( opt$output.report.dir, paste( "data", opt$output.format, sep="." ) ) } | |
44 if ( is.null( opt$gen.new.dgram ) ) { | |
45 opt$gen.new.dgram <- FALSE | |
46 } else { | |
47 if ( ! opt$gen.new.dgram %in% c( "no", "yes" ) ) { | |
48 stop( "invalid input to gen.new.dgram param", opt$gen.new.dgram, "\n" ) | |
49 } | |
50 ## set to TRUE/FALSE | |
51 opt$gen.new.dgram <- ( opt$gen.new.dgram == "yes" ) | |
52 } | |
53 | |
54 | |
55 load( opt$dataset ) ## should load the cl, treecl.res (or partcl.res) and data | |
56 | |
57 if ( ! exists( 'data' ) ) stop( "No data object in the rdata file provided for", opt$output.format, "format!!\n" ) | |
58 if ( inherits( data, "dist" ) ) stop( "data provided is a distance matrix - not a data matrix. Can't generate TreeView or Tab-delimited files w/distance matrices!\n" ) | |
59 | |
60 ## the rest of this is for the remaining output formats | |
61 ## pre-set the cluster results for rows & cols to NULL | |
62 direction <- NULL | |
63 if ( exists( 'treecl.res' ) ) { | |
64 cl.res <- treecl.res | |
65 if ( is.null( treecl.res$dist.method ) ) treecl.res$dist.method <- 'euclidean' # just set it to some stub so that the ctc fn's don't complain | |
66 } else { | |
67 if ( exists( 'partcl.res' ) ) { | |
68 cl.res <- partcl.res | |
69 } | |
70 else { | |
71 stop( 'could not find a valid cluster result to use for primary direction\n' ) | |
72 } | |
73 } | |
74 | |
75 if ( all( names( cl ) %in% rownames( data ) ) ) { | |
76 direction <- "rows" | |
77 } else if ( all( names( cl ) %in% colnames( data ) ) ) { | |
78 direction <- "cols" | |
79 data <- t( data ) | |
80 } else { | |
81 stop( "Specified cluster result does not come from this data set\n" ) | |
82 } | |
83 | |
84 | |
85 centroids <- NULL | |
86 cl <- sort( cl ) | |
87 if ( inherits( cl.res, "kmeans" ) ) { | |
88 ## already comes pre-calculated for us!! | |
89 centroids <- cl.res$centers | |
90 } else { | |
91 data <- data[ names( cl ), ] | |
92 cl.list <- unique( cl ) | |
93 cl.list <- lapply( cl.list, function(i) cl[ cl %in% i ] ) | |
94 centroids <- sapply( cl.list, | |
95 function(x) { | |
96 return( apply( data[ names(x), ], 2, mean, na.rm=T ) ) | |
97 } | |
98 ) | |
99 centroids <- t( centroids ) ## get them back to the same number of columns that data has now | |
100 } | |
101 | |
102 data <- centroids | |
103 colnames( data ) <- sapply( 1:max( cl ), function(i) sprintf( "cluster-%02d", i ) ) | |
104 | |
105 if ( opt$gen.new.dgram ) { | |
106 distance <- 'euclidean' | |
107 if ( inherits( cl.res, 'hclust' ) ) { | |
108 distance <- cl.res$dist.method | |
109 } | |
110 dmat <- Dist( data, distance ) | |
111 treecl.res <- hclust( dmat ) | |
112 cl <- cutree( treecl.res ) | |
113 } | |
114 | |
115 if ( direction == "cols" ) { | |
116 data <- t( data ) | |
117 } | |
118 | |
119 save( file=opt$output.name, treecl.res, cl, data ) |