annotate cluster.tools/dichotomize.sample.clusters.R @ 9:a3c03541fe6f draft default tip

Uploaded
author peter-waltman
date Mon, 11 Mar 2013 17:30:48 -0400
parents 2efa1a284546
children
Ignore whitespace changes - Everywhere: Within whitespace: At end of lines:
rev   line source
7
2efa1a284546 Uploaded
peter-waltman
parents:
diff changeset
1 #!/usr/bin/env Rscript
2efa1a284546 Uploaded
peter-waltman
parents:
diff changeset
2 argspec <- c("tab.2.cdt.R converts a data matrix to cdt format
2efa1a284546 Uploaded
peter-waltman
parents:
diff changeset
3
2efa1a284546 Uploaded
peter-waltman
parents:
diff changeset
4 Usage:
2efa1a284546 Uploaded
peter-waltman
parents:
diff changeset
5 tab.2.cdt.R -d <data.file>
2efa1a284546 Uploaded
peter-waltman
parents:
diff changeset
6 Optional:
2efa1a284546 Uploaded
peter-waltman
parents:
diff changeset
7 -o <output_file>
2efa1a284546 Uploaded
peter-waltman
parents:
diff changeset
8 \n\n")
2efa1a284546 Uploaded
peter-waltman
parents:
diff changeset
9 args <- commandArgs(TRUE)
2efa1a284546 Uploaded
peter-waltman
parents:
diff changeset
10 if ( length( args ) == 1 && args =="--help") {
2efa1a284546 Uploaded
peter-waltman
parents:
diff changeset
11 write(argspec, stderr())
2efa1a284546 Uploaded
peter-waltman
parents:
diff changeset
12 q();
2efa1a284546 Uploaded
peter-waltman
parents:
diff changeset
13 }
2efa1a284546 Uploaded
peter-waltman
parents:
diff changeset
14
2efa1a284546 Uploaded
peter-waltman
parents:
diff changeset
15 ## some helper fn's
2efa1a284546 Uploaded
peter-waltman
parents:
diff changeset
16 write.2.tab <- function( mat,
2efa1a284546 Uploaded
peter-waltman
parents:
diff changeset
17 fname ) {
2efa1a284546 Uploaded
peter-waltman
parents:
diff changeset
18 mat <- rbind( colnames( mat ), mat )
2efa1a284546 Uploaded
peter-waltman
parents:
diff changeset
19 mat <- cbind( c( "ID", rownames( mat )[-1] ),
2efa1a284546 Uploaded
peter-waltman
parents:
diff changeset
20 mat )
2efa1a284546 Uploaded
peter-waltman
parents:
diff changeset
21 write.table( mat, fname, sep="\t", row.names=FALSE, col.names=FALSE, quote=FALSE )
2efa1a284546 Uploaded
peter-waltman
parents:
diff changeset
22 }
2efa1a284546 Uploaded
peter-waltman
parents:
diff changeset
23
2efa1a284546 Uploaded
peter-waltman
parents:
diff changeset
24 lib.load.quiet <- function( package ) {
2efa1a284546 Uploaded
peter-waltman
parents:
diff changeset
25 package <- as.character(substitute(package))
2efa1a284546 Uploaded
peter-waltman
parents:
diff changeset
26 suppressPackageStartupMessages( do.call( "library", list( package=package ) ) )
2efa1a284546 Uploaded
peter-waltman
parents:
diff changeset
27 }
2efa1a284546 Uploaded
peter-waltman
parents:
diff changeset
28
2efa1a284546 Uploaded
peter-waltman
parents:
diff changeset
29 lib.load.quiet( getopt )
2efa1a284546 Uploaded
peter-waltman
parents:
diff changeset
30 lib.load.quiet( ctc )
2efa1a284546 Uploaded
peter-waltman
parents:
diff changeset
31 if ( any( c( 'flashClust', 'fastcluster' ) %in% installed.packages() ) ) {
2efa1a284546 Uploaded
peter-waltman
parents:
diff changeset
32 if ( 'flashClust' %in% installed.packages() ) {
2efa1a284546 Uploaded
peter-waltman
parents:
diff changeset
33 lib.load.quiet( flashClust )
2efa1a284546 Uploaded
peter-waltman
parents:
diff changeset
34 } else {
2efa1a284546 Uploaded
peter-waltman
parents:
diff changeset
35 if ( 'fastcluster' %in% installed.packages() ) {
2efa1a284546 Uploaded
peter-waltman
parents:
diff changeset
36 lib.load.quiet( fastcluster )
2efa1a284546 Uploaded
peter-waltman
parents:
diff changeset
37 }
2efa1a284546 Uploaded
peter-waltman
parents:
diff changeset
38 }
2efa1a284546 Uploaded
peter-waltman
parents:
diff changeset
39 }
2efa1a284546 Uploaded
peter-waltman
parents:
diff changeset
40
2efa1a284546 Uploaded
peter-waltman
parents:
diff changeset
41 spec <- matrix( c( "dataset", "d", 1, "character",
2efa1a284546 Uploaded
peter-waltman
parents:
diff changeset
42 "num.k", "k", 2, "character",
2efa1a284546 Uploaded
peter-waltman
parents:
diff changeset
43 "output.fname", "o", 2, "character"
2efa1a284546 Uploaded
peter-waltman
parents:
diff changeset
44 ),
2efa1a284546 Uploaded
peter-waltman
parents:
diff changeset
45 nc=4,
2efa1a284546 Uploaded
peter-waltman
parents:
diff changeset
46 byrow=TRUE
2efa1a284546 Uploaded
peter-waltman
parents:
diff changeset
47 )
2efa1a284546 Uploaded
peter-waltman
parents:
diff changeset
48
2efa1a284546 Uploaded
peter-waltman
parents:
diff changeset
49
2efa1a284546 Uploaded
peter-waltman
parents:
diff changeset
50 opt <- getopt( spec=spec )
2efa1a284546 Uploaded
peter-waltman
parents:
diff changeset
51 if ( is.null( opt$output.fname ) ) { opt$output.fname <- file.path( opt$output.report.dir, paste( "data", opt$output.format, sep="." ) ) }
2efa1a284546 Uploaded
peter-waltman
parents:
diff changeset
52 if ( is.null( opt$num.k ) ) {
2efa1a284546 Uploaded
peter-waltman
parents:
diff changeset
53 opt$num.k <- -1
2efa1a284546 Uploaded
peter-waltman
parents:
diff changeset
54 } else {
2efa1a284546 Uploaded
peter-waltman
parents:
diff changeset
55 num.k <- as.integer( eval( parse( text=paste( "c(",gsub( "-", ":", gsub( ", |,", ",", opt$num.k ) ), ")" ) ) ) )
2efa1a284546 Uploaded
peter-waltman
parents:
diff changeset
56 num.k <- num.k[ ! is.na( num.k ) ]
2efa1a284546 Uploaded
peter-waltman
parents:
diff changeset
57 if ( length( opt$num.k ) == 0 ) stop( 'invalid input for k_range specified:', opt$num.k, "\n" )
2efa1a284546 Uploaded
peter-waltman
parents:
diff changeset
58
2efa1a284546 Uploaded
peter-waltman
parents:
diff changeset
59 num.k <- num.k[ ! num.k %in% 1 ] # strip out a k==1 since that doesn't make any sense
2efa1a284546 Uploaded
peter-waltman
parents:
diff changeset
60 opt$num.k <- num.k; rm( num.k )
2efa1a284546 Uploaded
peter-waltman
parents:
diff changeset
61 }
2efa1a284546 Uploaded
peter-waltman
parents:
diff changeset
62
2efa1a284546 Uploaded
peter-waltman
parents:
diff changeset
63 load( opt$dataset )
2efa1a284546 Uploaded
peter-waltman
parents:
diff changeset
64 ## if this is a clustering result w/cluster assignments ('raw' CCPLUS does not)
2efa1a284546 Uploaded
peter-waltman
parents:
diff changeset
65 if ( exists( 'cl' ) ) {
2efa1a284546 Uploaded
peter-waltman
parents:
diff changeset
66 k <- max( as.numeric( cl ) )
2efa1a284546 Uploaded
peter-waltman
parents:
diff changeset
67 cl <- matrix( cl, nc=1, dimnames=list( names(cl), k ) )
2efa1a284546 Uploaded
peter-waltman
parents:
diff changeset
68 if ( (length(opt$num.k)==1) && (opt$num.k == -1 ) ) opt$num.k <- k
2efa1a284546 Uploaded
peter-waltman
parents:
diff changeset
69
2efa1a284546 Uploaded
peter-waltman
parents:
diff changeset
70 ## if this is a one-off to produce a phenotype for the number of clusters that the user originally proposed
2efa1a284546 Uploaded
peter-waltman
parents:
diff changeset
71 if ( !opt$num.k[1] %in% c( -1, k ) ) {
2efa1a284546 Uploaded
peter-waltman
parents:
diff changeset
72
2efa1a284546 Uploaded
peter-waltman
parents:
diff changeset
73 if ( exists( 'partcl.res' ) || exists( 'select.result' ) ) {
2efa1a284546 Uploaded
peter-waltman
parents:
diff changeset
74
2efa1a284546 Uploaded
peter-waltman
parents:
diff changeset
75 if ( exists( 'partcl.res' ) ) {
2efa1a284546 Uploaded
peter-waltman
parents:
diff changeset
76 warning( 'The k_range value(s) specified are:',
2efa1a284546 Uploaded
peter-waltman
parents:
diff changeset
77 opt$num.k,
2efa1a284546 Uploaded
peter-waltman
parents:
diff changeset
78 "however k_range vals can not specify alternate k values for partition clusters. Using the K value that corresponds to this result instead\n" )
2efa1a284546 Uploaded
peter-waltman
parents:
diff changeset
79 } else {
2efa1a284546 Uploaded
peter-waltman
parents:
diff changeset
80 warning( 'The k_range value(s) specified are:',
2efa1a284546 Uploaded
peter-waltman
parents:
diff changeset
81 opt$num.k,
2efa1a284546 Uploaded
peter-waltman
parents:
diff changeset
82 "however k_range vals can not specify alternate k values for specific cluster results from CCPLUS (i.e. those from the Select K or Extract tools). To get alternate K values, re-run the dichotomizer on the 'raw' CCPLUS results. Using the K value that corresponds to this result instead\n" )
2efa1a284546 Uploaded
peter-waltman
parents:
diff changeset
83 }
2efa1a284546 Uploaded
peter-waltman
parents:
diff changeset
84
2efa1a284546 Uploaded
peter-waltman
parents:
diff changeset
85 opt$num.k <- k
2efa1a284546 Uploaded
peter-waltman
parents:
diff changeset
86 cl <- matrix( cl, nr=1, dimnames=list( k, names(cl) ) )
2efa1a284546 Uploaded
peter-waltman
parents:
diff changeset
87 } else {
2efa1a284546 Uploaded
peter-waltman
parents:
diff changeset
88 ## handle if this is a hclust result
2efa1a284546 Uploaded
peter-waltman
parents:
diff changeset
89
2efa1a284546 Uploaded
peter-waltman
parents:
diff changeset
90 opt$num.k <- opt$num.k[ opt$num.k < length( cl ) ]
2efa1a284546 Uploaded
peter-waltman
parents:
diff changeset
91 cl.samps <- rownames( cl )
2efa1a284546 Uploaded
peter-waltman
parents:
diff changeset
92 cl <- sapply( opt$num.k, function(i) cutree( treecl.res, i )[ cl.samps ] )
2efa1a284546 Uploaded
peter-waltman
parents:
diff changeset
93 colnames( cl ) <- opt$num.k
2efa1a284546 Uploaded
peter-waltman
parents:
diff changeset
94 }
2efa1a284546 Uploaded
peter-waltman
parents:
diff changeset
95 }
2efa1a284546 Uploaded
peter-waltman
parents:
diff changeset
96 } else if ( exists( 'results' ) ) {
2efa1a284546 Uploaded
peter-waltman
parents:
diff changeset
97 ## handle if this is a ccplus-raw result
2efa1a284546 Uploaded
peter-waltman
parents:
diff changeset
98 opt$num.k <- opt$num.k[ opt$num.k <= length( results ) ]
2efa1a284546 Uploaded
peter-waltman
parents:
diff changeset
99 cl <- sapply( results[ opt$num.k ], '[[', 'consensusClass' )
2efa1a284546 Uploaded
peter-waltman
parents:
diff changeset
100 colnames( cl ) <- opt$num.k
2efa1a284546 Uploaded
peter-waltman
parents:
diff changeset
101 }
2efa1a284546 Uploaded
peter-waltman
parents:
diff changeset
102
2efa1a284546 Uploaded
peter-waltman
parents:
diff changeset
103 pheno.mat <- lapply( 1:ncol(cl),
2efa1a284546 Uploaded
peter-waltman
parents:
diff changeset
104 function(i) {
2efa1a284546 Uploaded
peter-waltman
parents:
diff changeset
105 x <- cl[,i]
2efa1a284546 Uploaded
peter-waltman
parents:
diff changeset
106 cls <- ks <- sort( unique(x) )
2efa1a284546 Uploaded
peter-waltman
parents:
diff changeset
107 cls <- sapply( cls, function(y) as.numeric( x %in% y ) )
2efa1a284546 Uploaded
peter-waltman
parents:
diff changeset
108 colnames(cls) <- paste( "CLeq", ks, sep="" )
2efa1a284546 Uploaded
peter-waltman
parents:
diff changeset
109 rownames(cls) <- names(x)
2efa1a284546 Uploaded
peter-waltman
parents:
diff changeset
110 return(cls)
2efa1a284546 Uploaded
peter-waltman
parents:
diff changeset
111 }
2efa1a284546 Uploaded
peter-waltman
parents:
diff changeset
112 )
2efa1a284546 Uploaded
peter-waltman
parents:
diff changeset
113 names( pheno.mat ) <- opt$num.k
2efa1a284546 Uploaded
peter-waltman
parents:
diff changeset
114
2efa1a284546 Uploaded
peter-waltman
parents:
diff changeset
115 final.mat <- matrix( NA, nc=0, nrow=nrow(cl), dimnames=list( names(cl), NULL ) )
2efa1a284546 Uploaded
peter-waltman
parents:
diff changeset
116 for ( i in names( pheno.mat ) ) {
2efa1a284546 Uploaded
peter-waltman
parents:
diff changeset
117 colnames( pheno.mat[[i]] ) <- paste( "Keq", i, "_", colnames( pheno.mat[[i]] ), sep="" )
2efa1a284546 Uploaded
peter-waltman
parents:
diff changeset
118 final.mat <- cbind( final.mat, pheno.mat[[i]] )
2efa1a284546 Uploaded
peter-waltman
parents:
diff changeset
119 }
2efa1a284546 Uploaded
peter-waltman
parents:
diff changeset
120
2efa1a284546 Uploaded
peter-waltman
parents:
diff changeset
121 write.2.tab( final.mat, opt$output.fname )