|
5
|
1 #!/usr/bin/env Rscript
|
|
|
2 ## Script by Peter Waltman
|
|
|
3 ##
|
|
|
4 ## License under Creative Commons Attribution 3.0 Unported (CC BY 3.0)
|
|
|
5 ##
|
|
|
6 argspec <- c("TBD
|
|
|
7
|
|
|
8
|
|
|
9 Usage:
|
|
|
10 normalize.matrix.R -d <data.file>
|
|
|
11 Optional:
|
|
|
12 -o <output.name>
|
|
|
13 TBD
|
|
|
14 \n\n")
|
|
|
15
|
|
|
16 args <- commandArgs(TRUE)
|
|
|
17 if ( length( args ) == 1 && args =="--help") {
|
|
|
18 write(argspec, stderr())
|
|
|
19 q();
|
|
|
20 }
|
|
|
21
|
|
|
22 ## some helper fn's
|
|
|
23 write.2.tab <- function( mat,
|
|
|
24 fname ) {
|
|
|
25 mat <- rbind( colnames( mat ), mat )
|
|
|
26 mat <- cbind( c( "ID", rownames( mat )[-1] ),
|
|
|
27 mat )
|
|
|
28 write.table( mat, fname, sep="\t", row.names=FALSE, col.names=FALSE, quote=FALSE )
|
|
|
29 }
|
|
|
30
|
|
|
31 lib.load.quiet <- function( package ) {
|
|
|
32 package <- as.character(substitute(package))
|
|
|
33 suppressPackageStartupMessages( do.call( "library", list( package=package ) ) )
|
|
|
34 }
|
|
|
35 lib.load.quiet(getopt)
|
|
|
36
|
|
|
37
|
|
|
38 spec <- matrix( c( "data.fname", "d", 1, "character",
|
|
|
39 "data.fname2", "D", 1, "character",
|
|
|
40 "ctr.method", "c", 2, "character",
|
|
|
41 "output.fname", "o", 2, "character",
|
|
|
42 "return.normals", "r", 2, "character",
|
|
|
43 "out.norm.fname", "O", 2, "character"
|
|
|
44 ),
|
|
|
45 nc=4,
|
|
|
46 byrow=TRUE
|
|
|
47 )
|
|
|
48
|
|
|
49 opt <- getopt( spec=spec )
|
|
|
50
|
|
|
51 #set some reasonable defaults for the options that are needed,
|
|
|
52 #but were not specified.
|
|
|
53 if ( is.null( opt$output.fname ) ) { opt$output.fname <- 'merge_merge.tumors.tab' }
|
|
|
54 if ( is.null( opt$ctr.method ) ) { opt$ctr.method <- 'median' }
|
|
|
55 if ( is.null( opt$return.normals ) ) {
|
|
|
56 opt$return.normals <- FALSE
|
|
|
57 } else {
|
|
|
58 opt$return.normals <- ( tolower( opt$return.normals ) %in% "yes" )
|
|
|
59 }
|
|
|
60 if ( is.null( opt$out.norm.fname ) ) { opt$out.norm.fname <- 'merge_merge.normals.tab' }
|
|
|
61
|
|
|
62
|
|
|
63 data <- as.matrix( read.delim( opt$data.fname, header=T, row.names=1 , check.names=FALSE ) )
|
|
|
64 normer.data <- as.matrix( read.delim( opt$data.fname2, header=T, row.names=1 , check.names=FALSE ) )
|
|
|
65 if ( ! all( rownames(data) %in% rownames(normer.data) ) ) {
|
|
|
66 gene.ovp <- rownames(data)
|
|
|
67 gene.ovp <- gene.ovp[ gene.ovp %in% rownames(normer.data) ]
|
|
|
68 if ( length( gene.ovp ) == 0 ) stop( 'No gene overlap between specified matrices\n' )
|
|
|
69
|
|
|
70 data <- data[ gene.ovp, , drop=F ]
|
|
|
71 normer.data <- normer.data[ gene.ovp, , drop=F ]
|
|
|
72 }
|
|
|
73
|
|
|
74 my.ctr.vals <- apply( normer.data, 1, get( opt$ctr.method ), na.rm=T )
|
|
|
75
|
|
|
76 data <- apply( data, 2, function(x) x-my.ctr.vals )
|
|
|
77
|
|
|
78 write.2.tab( data, opt$output.fname )
|
|
|
79 if ( opt$return.normals ) {
|
|
|
80 normer.data <- normer.data - my.ctr.vals
|
|
|
81 write.2.tab( normer.data, opt$out.norm.fname )
|
|
|
82 }
|
|
|
83
|
|
|
84
|
