|
0
|
1 #!/usr/bin/Rscript
|
|
|
2
|
|
|
3 #usage, options and doc goes here
|
|
|
4 argspec <- c("normalize.r - takes any flat file and normalizes the rows or the columns using various normalizations (median_shift, mean_shift, t_statistic (z-score), exp_fit, normal_fit, weibull_0.5_fit, weibull_1_fit, weibull_1.5_fit, weibull_5_fit). Requires a single header line and a single cloumn of annotation.
|
|
|
5 Usage:
|
|
|
6 normalize.r input.tab norm_type norm_by > output.tab
|
|
|
7 Example:
|
|
|
8 Rscript normalize.r test_matrix.tab median_shift column > output.tab
|
|
|
9 Rscript normalize.r test_matrix.tab mean_shift row normals.tab > output.tab
|
|
|
10 Options:
|
|
|
11 input matrix (annotated by row and column names)
|
|
|
12 normalization type; available options:
|
|
|
13 median_shift - shifts all values by the median or the row/column if no normals are specified, otherwise shifts by the median of normals
|
|
|
14 mean_shift - shifts all values by the mean or the row/column if no normals are specified, otherwise shifts by the mean of normals
|
|
|
15 t_statistic - converts all values to z-scores; if normals are specified then converts to z-scores within normal and non-normal classes separately
|
|
|
16 exp_fit - (only by column) ranks data and transforms exponential CDF
|
|
|
17 normal_fit - (only by column) ranks data and transforms normal CDF
|
|
|
18 weibull_0.5_fit - (only by column) ranks data and transforms Weibull CDF with scale parameter = 1 and shape parameter = 0.5
|
|
|
19 weibull_1_fit - (only by column) ranks data and transforms Weibull CDF with scale parameter = 1 and shape parameter = 1
|
|
|
20 weibull_1.5_fit - (only by column) ranks data and transforms Weibull CDF with scale parameter = 1 and shape parameter = 1.5
|
|
|
21 weibull_5_fit - (only by column) ranks data and transforms Weibull CDF with scale parameter = 1 and shape parameter = 5
|
|
|
22 normalization by:
|
|
|
23 row
|
|
|
24 column
|
|
|
25 normals_file is an optional parameter which contains a list of column headers from the input matrix, which should be considered as normals
|
|
|
26 output file is specified through redirect character >")
|
|
|
27
|
|
|
28 read_matrix <- function(in_file){
|
|
|
29 header <- strsplit(readLines(con=in_file, n=1), "\t")[[1]]
|
|
|
30 cl.cols<- 1:length(header) > 1
|
|
|
31 data_matrix.df <- read.delim(in_file, header=TRUE, row.names=NULL, stringsAsFactors=FALSE, na.strings="NA", check.names=FALSE)
|
|
|
32 data_matrix <- as.matrix(data_matrix.df[,cl.cols])
|
|
|
33 rownames(data_matrix) <- data_matrix.df[,1]
|
|
|
34 return(data_matrix)
|
|
|
35
|
|
|
36 #read_mtrx <- as.matrix(read.table(in_file, header=TRUE, sep="", row.names=NULL, stringsAsFactors=FALSE, na.strings="NA")) #separate on white characters
|
|
|
37 #read_mtrx[,1]
|
|
|
38
|
|
|
39 #return(as.matrix(read.table(in_file, header=TRUE, sep="", row.names=1))) #separate on white characters
|
|
|
40 #mtrx <- read.delim(in_file, header=TRUE, sep="", row.names=NULL, stringsAsFactors=FALSE, na.strings="NA")
|
|
|
41 #print(mtrx[1,])
|
|
|
42 }
|
|
|
43
|
|
|
44 write_matrix <- function(data_matrix){
|
|
|
45 header <- append(c("Genes"), colnames(data_matrix))
|
|
|
46 write.table(t(header), stdout(), quote=FALSE, sep="\t", row.names=FALSE, col.names=FALSE)
|
|
|
47 write.table(data_matrix, stdout(), quote=FALSE, sep="\t", row.names=TRUE, col.names=FALSE)
|
|
|
48 }
|
|
|
49
|
|
|
50 read_normals <- function(in_file){
|
|
|
51 return(as.matrix(read.table(in_file, header=FALSE, sep="", as.is = TRUE))[, 1])
|
|
|
52 }
|
|
|
53
|
|
|
54 normalize <- function(data_matrix, norm_type, normals_list, tumors_list){
|
|
|
55 if(norm_type == 'MEDIAN_SHIFT'){
|
|
|
56 return(shift(data_matrix, 'MEDIAN', normals_list, tumors_list))
|
|
|
57 }
|
|
|
58 else if(norm_type == 'MEAN_SHIFT'){
|
|
|
59 return(shift(data_matrix, 'MEAN', normals_list, tumors_list))
|
|
|
60 }
|
|
|
61 else if(norm_type == 'T_STATISTIC'){
|
|
|
62 return(compute_z_score(data_matrix, normals_list, tumors_list))
|
|
|
63 }
|
|
|
64 else if(norm_type == 'EXPONENTIAL_FIT'){
|
|
|
65 return(fit_distribution(data_matrix, 'EXPONENTIAL'))
|
|
|
66 }
|
|
|
67 else if(norm_type == 'NORMAL_FIT'){
|
|
|
68 return(fit_distribution(data_matrix, 'NORMAL'))
|
|
|
69 }
|
|
|
70 else if(norm_type == 'WEIBULL_0.5_FIT'){
|
|
|
71 return(fit_distribution(data_matrix, 'WEIBULL_0.5'))
|
|
|
72 }
|
|
|
73 else if(norm_type == 'WEIBULL_1_FIT'){
|
|
|
74 return(fit_distribution(data_matrix, 'WEIBULL_1'))
|
|
|
75 }
|
|
|
76 else if(norm_type == 'WEIBULL_1.5_FIT'){
|
|
|
77 return(fit_distribution(data_matrix, 'WEIBULL_1.5'))
|
|
|
78 }
|
|
|
79 else if(norm_type == 'WEIBULL_5_FIT'){
|
|
|
80 return(fit_distribution(data_matrix, 'WEIBULL_5'))
|
|
|
81 }
|
|
|
82 }
|
|
|
83
|
|
|
84 shift <- function(data_matrix, shift_type, normals_list, tumors_list){
|
|
|
85 return(t(apply(data_matrix, 1, shift_normalize_row, norm_type=shift_type, normals_list=normals_list, tumors_list=tumors_list)))
|
|
|
86 }
|
|
|
87
|
|
|
88 shift_normalize_row <- function(data_row, norm_type, normals_list, tumors_list){
|
|
|
89 if(length(normals_list) == 0){ #no normals are specified
|
|
|
90 if(norm_type == 'MEDIAN'){
|
|
|
91 row_stat <- median(data_row)
|
|
|
92 }
|
|
|
93 else if(norm_type == 'MEAN'){
|
|
|
94 row_stat <- mean(data_row)
|
|
|
95 }
|
|
|
96 return(unlist(lapply(data_row, function(x){return(x - row_stat);})))
|
|
|
97 }
|
|
|
98 else{ #normals are specified
|
|
|
99 normal_values <- data_row[normals_list]
|
|
|
100 tumor_columns <- data_row[tumors_list]
|
|
|
101
|
|
|
102 if(norm_type == 'MEDIAN'){
|
|
|
103 row_stat <- median(normal_values)
|
|
|
104 }
|
|
|
105 else if(norm_type == 'MEAN'){
|
|
|
106 row_stat <- mean(normal_values)
|
|
|
107 }
|
|
|
108 return(unlist(lapply(tumor_columns, function(x){return(x - row_stat);})))
|
|
|
109 }
|
|
|
110 }
|
|
|
111
|
|
|
112 compute_z_score <- function(data_matrix, normals_list, tumors_list){
|
|
|
113 return(t(apply(data_matrix, 1, t_stat_normalize_row, normals_list=normals_list, tumors_list=tumors_list)))
|
|
|
114 }
|
|
|
115
|
|
|
116 t_stat_normalize_row <- function(data_row, normals_list, tumors_list){
|
|
|
117 if(length(normals_list) == 0){ #no normals are specified
|
|
|
118 row_mean <- mean(data_row)
|
|
|
119 row_sd <- sd(data_row)
|
|
|
120 return(unlist(lapply(data_row, function(x){return((x - row_mean)/row_sd);})))
|
|
|
121 }
|
|
|
122 else{ #normals are specified
|
|
|
123 normal_values <- data_row[normals_list]
|
|
|
124 normal_mean <- mean(normal_values)
|
|
|
125 normal_sd <- sd(normal_values)
|
|
|
126 normalized_normals <- unlist(lapply(normal_values, function(x){return((x - normal_mean)/normal_sd);}))
|
|
|
127
|
|
|
128 tumor_values <- data_row[tumors_list]
|
|
1
|
129 normalized_tumors <- unlist(lapply(tumor_values, function(x){return((x - normal_mean)/normal_sd);}))
|
|
0
|
130
|
|
|
131 return(append(normalized_normals, normalized_tumors))
|
|
|
132 }
|
|
|
133 }
|
|
|
134
|
|
|
135 rankNA <- function(col){ #originally written by Dan Carlin
|
|
|
136 col[!is.na(col)]<-(rank(col[!is.na(col)])/sum(!is.na(col)))-(1/sum(!is.na(col)))
|
|
|
137 return(col)
|
|
|
138 }
|
|
|
139
|
|
|
140 fit_distribution <- function(data_matrix, dist){
|
|
|
141 if(dist == 'EXPONENTIAL'){
|
|
|
142 ranked_data_matrix <- apply(data_matrix,2,rankNA) #idea by Dan Carlin
|
|
|
143 return(apply(ranked_data_matrix, c(1,2), qexp))
|
|
|
144 }
|
|
|
145 else if(dist == 'NORMAL'){
|
|
|
146 ranked_data_matrix <- apply(data_matrix,2,rankNA)
|
|
|
147 #return(apply(ranked_data_matrix, c(1,2), function(x){return(qnorm(mean=mean(x), sd=sd(x)));}))
|
|
1
|
148 return(apply(ranked_data_matrix, c(1,2), qnorm, mean=0, sd=1))
|
|
0
|
149 }
|
|
|
150 else if(dist == 'WEIBULL_0.5'){
|
|
|
151 ranked_data_matrix <- apply(data_matrix,2,rankNA)
|
|
|
152 return(apply(ranked_data_matrix, c(1,2), qweibull, scale=1, shape=0.5))
|
|
|
153 }
|
|
|
154 else if(dist == 'WEIBULL_1'){
|
|
|
155 ranked_data_matrix <- apply(data_matrix,2,rankNA)
|
|
|
156 return(apply(ranked_data_matrix, c(1,2), qweibull, scale=1, shape=1))
|
|
|
157 }
|
|
|
158 else if(dist == 'WEIBULL_1.5'){
|
|
|
159 ranked_data_matrix <- apply(data_matrix,2,rankNA)
|
|
|
160 return(apply(ranked_data_matrix, c(1,2), qweibull, scale=1, shape=1.5))
|
|
|
161 }
|
|
|
162 else if(dist == 'WEIBULL_5'){
|
|
|
163 ranked_data_matrix <- apply(data_matrix,2,rankNA)
|
|
|
164 return(apply(ranked_data_matrix, c(1,2), qweibull, scale=1, shape=5))
|
|
|
165 }
|
|
|
166 }
|
|
|
167
|
|
|
168 main <- function(argv) {
|
|
|
169 #determine if correct number of arguments are specified and if normals are specified
|
|
|
170 with_normals = FALSE
|
|
|
171
|
|
|
172 if(length(argv) == 1){
|
|
|
173 if(argv==c('--help')){
|
|
|
174 write(argspec, stderr());
|
|
|
175 q();
|
|
|
176 }
|
|
|
177 }
|
|
|
178
|
|
|
179 if(!(length(argv) == 3 || length(argv) == 4)){
|
|
|
180 write("ERROR: invalid number of arguments is specified", stderr());
|
|
|
181 q();
|
|
|
182 }
|
|
|
183
|
|
|
184 if(length(argv) == 4){
|
|
|
185 with_normals = TRUE
|
|
|
186 normals_file <- argv[4]
|
|
|
187 }
|
|
|
188
|
|
|
189 #store command line arguments in variables:
|
|
|
190 input_file <- argv[1]
|
|
|
191 norm_type <- toupper(argv[2])
|
|
|
192 norm_by <- toupper(argv[3])
|
|
|
193
|
|
|
194 #read the input file(s):
|
|
|
195 data_matrix <- read_matrix(input_file)
|
|
|
196
|
|
|
197 if(with_normals){
|
|
|
198 normals_list <- read_normals(normals_file)
|
|
|
199 normals_indices <- which(colnames(data_matrix) %in% normals_list)
|
|
|
200 tumor_indices <- which(!(colnames(data_matrix) %in% normals_list))
|
|
|
201 norm_by <- 'ROW'
|
|
|
202 }
|
|
|
203 else{
|
|
|
204 normals_indices <- c()
|
|
|
205 tumor_indices <- c()
|
|
|
206 }
|
|
|
207
|
|
|
208 #if normalize by columns then transpose the matrix:
|
|
|
209 if(norm_by == 'COLUMN'){
|
|
|
210 data_matrix <- t(data_matrix)
|
|
|
211 }
|
|
|
212
|
|
|
213 #normalize:
|
|
|
214 data_matrix <- normalize(data_matrix, norm_type, normals_indices, tumor_indices)
|
|
|
215
|
|
|
216 #if normalize by columns then transpose the matrix again since we normalized the transposed matrix by row:
|
|
|
217 if(norm_by == 'COLUMN'){
|
|
|
218 data_matrix <- t(data_matrix)
|
|
|
219 }
|
|
|
220
|
|
|
221 write_matrix(data_matrix)
|
|
|
222 #print(data_matrix)
|
|
|
223 }
|
|
|
224
|
|
|
225 main(commandArgs(TRUE))
|
