1
|
1 #########################################################################
|
|
2 # SCRIPT INTENSITY CHECK #
|
|
3 # #
|
|
4 # Input: Data Matrix, VariableMetadata, SampleMetadata #
|
3
|
5 # Output: VariableMetadata, Graphics #
|
1
|
6 # #
|
|
7 # Dependencies: RcheckLibrary.R #
|
|
8 # #
|
|
9 #########################################################################
|
0
|
10
|
|
11
|
|
12 # Parameters (for dev)
|
|
13 if(FALSE){
|
|
14
|
|
15 rm(list = ls())
|
3
|
16 setwd("Y:\\Developpement\\Intensity check\\Pour tests\\Tests_global")
|
0
|
17
|
|
18 DM.name <- "DM_NA.tabular"
|
|
19 SM.name <- "SM_NA.tabular"
|
|
20 VM.name <- "vM_NA.tabular"
|
3
|
21 method <- "one_class"
|
|
22 chosen.stat <- "mean,sd,quartile,decile,NA"
|
6
|
23 class.col <- "uv"
|
3
|
24 test.fold <- "Yes"
|
|
25 class1 <- "Pools"
|
1
|
26 fold.frac <- "Top"
|
3
|
27 logarithm <- "log10"
|
0
|
28 VM.output <- "new_VM.txt"
|
|
29 graphs.output <- "Barplots_and_Boxplots.pdf"
|
|
30 }
|
|
31
|
1
|
32
|
|
33
|
|
34
|
3
|
35 intens_check <- function(DM.name, SM.name, VM.name, method, chosen.stat, class.col, test.fold, class1, fold.frac,
|
|
36 logarithm, VM.output, graphs.output){
|
1
|
37
|
3
|
38 # This function allows to check the intensities with various statistics, number of missing values and mean fold change
|
0
|
39 #
|
3
|
40 # Three methods proposed:
|
|
41 # - global: tests for each variable without distinction between samples
|
|
42 # - one class: one class versus all the remaining samples
|
|
43 # - each class: if the class columns contains at least three classes and you want to test each of them
|
0
|
44 #
|
|
45 # Parameters:
|
|
46 # DM.name, SM.name, VM.name: dataMatrix, sampleMetadata, variableMetadata files access
|
3
|
47 # method: "global", "one_class", "each_class"
|
|
48 # chosen.stat: character listing the chosen analysis (comma-separated)
|
|
49 # class.col: number of the sampleMetadata's column with classes (if method = one_class or each_class)
|
|
50 # test.fold: "yes" or "no" (if method = one_class or each_class)
|
|
51 # class1: name of the class (if method = one_class)
|
|
52 # fold.frac: "Top" -> class1/other or "Bottom" -> other/class1 (if method = one_class)
|
|
53 # logarithm: "log2", "log10" or "none" (if method = one_class or each_class)
|
1
|
54 # VM.output: output file's access (VM with new columns)
|
|
55 # graphs.output: pdf file's access with barplots for the proportion of NA and boxplots with the folds values
|
|
56
|
0
|
57
|
|
58
|
|
59
|
1
|
60 # Input ---------------------------------------------------------------------------------------------------
|
0
|
61
|
|
62 DM <- read.table(DM.name, header=TRUE, sep="\t", check.names=FALSE)
|
|
63 SM <- read.table(SM.name, header=TRUE, sep="\t", check.names=FALSE)
|
|
64 VM <- read.table(VM.name, header=TRUE, sep="\t", check.names=FALSE)
|
|
65
|
|
66
|
|
67
|
|
68 # Table match check with Rchecklibrary
|
|
69 table.check <- match3(DM, SM, VM)
|
|
70 check.err(table.check)
|
3
|
71
|
5
|
72 # Transposing the dataMatrix
|
0
|
73 rownames(DM) <- DM[,1]
|
|
74 var_names <- DM[,1]
|
|
75 DM <- DM[,-1]
|
|
76 DM <- data.frame(t(DM))
|
5
|
77
|
|
78 # Re-ordering the dataMatrix to match the sampleMetadata file order
|
|
79 DM <- merge(x=cbind(1:nrow(SM),SM), y=DM, by.x=2, by.y=0)
|
|
80 DM <- DM[order(DM[,2]),]
|
|
81 rownames(DM) <- DM[,1]
|
|
82 DM <- DM[,-c(1:(ncol(SM)+1))]
|
|
83
|
|
84
|
3
|
85 stat.list <- strsplit(chosen.stat,",")[[1]]
|
0
|
86
|
1
|
87
|
|
88 # check class.col, class1 and the number of classes ---------------------------------------------------------
|
0
|
89
|
3
|
90 #set 1 class for all samples in case of method = no_class
|
|
91 if(method=="no_class"){
|
|
92 c_class <- rep("global", length=nrow(DM))
|
|
93 classnames <- "global"
|
|
94 nb_class=1
|
|
95 test.fold <- "No"
|
0
|
96 }
|
3
|
97
|
|
98
|
|
99 if(method != "no_class"){
|
|
100
|
|
101 if(!(class.col %in% colnames(SM))){
|
6
|
102 stop("\n- - - - - - - - -\n", "The ",class.col, " column is not found in the specified sample metadata file.","\n- - - - - - - - -\n")
|
3
|
103 }
|
0
|
104
|
3
|
105 c_class <- SM[,class.col]
|
|
106 c_class <- as.factor(c_class)
|
|
107 nb_class <- nlevels(c_class)
|
|
108 classnames <- levels(c_class)
|
|
109
|
|
110 if((nb_class < 2)&&(test.fold=="Yes")){
|
6
|
111 err.1class <- c("\n The ",class.col, " column contains only one class, fold calculation could not be executed. \n")
|
3
|
112 cat(err.1class)
|
|
113 }
|
|
114
|
|
115 if((nb_class > (nrow(SM))/3)&&(method == "each_class")){
|
|
116 class.err <- c("\n There are too many classes, think about reducing the number of classes and excluding those
|
5
|
117 with few samples. \n")
|
3
|
118 cat(class.err)
|
|
119 }
|
|
120
|
|
121
|
|
122 if(method == "one_class"){
|
|
123 if(!(class1 %in% classnames)){
|
|
124 list.class1 <- c("\n Classes:",classnames,"\n")
|
|
125 cat(list.class1)
|
6
|
126 err.class1 <- c("The ",class1, " class does not appear in the ",class.col," column.")
|
3
|
127 stop("\n- - - - - - - - -\n", err.class1,"\n- - - - - - - - -\n")
|
|
128 }
|
0
|
129
|
3
|
130 #If method is "one_class", change others classes in "other"
|
|
131 for(i in 1:length(c_class)){
|
|
132 if(c_class[i]!=class1){
|
|
133 c_class <- as.character(c_class)
|
|
134 c_class[i] <- "Other"
|
|
135 c_class <- as.factor(c_class)
|
|
136 nb_class <- nlevels(c_class)
|
|
137 classnames <- c(class1,"Other")
|
|
138 }
|
|
139 }
|
|
140 }
|
|
141
|
1
|
142 }
|
3
|
143
|
|
144
|
|
145 # Statistics ------------------------------------------------------------------------------------------------
|
|
146
|
5
|
147 # Check whether the dataMatrix contains non-numeric values
|
|
148 if(!(is.numeric(as.matrix(DM)))){
|
|
149 # findchar definition
|
|
150 findchar <- function(myval){
|
|
151 if(is.na(myval)){
|
|
152 return("ok")
|
|
153 }else{
|
|
154 mytest <- as.character(myval)
|
|
155 if(is.na(as.numeric(mytest[1]))){
|
|
156 return("char")
|
|
157 }else{
|
|
158 return("ok")
|
|
159 }
|
|
160 }
|
|
161 }
|
|
162 # findchar application
|
|
163 chardiag <- suppressWarnings(apply(DM,2,vapply,findchar,"character"))
|
|
164 charlist <- which(chardiag == "char")
|
|
165 err.stock <- paste("\n- - - - - - - - -\nYour dataMatrix contains",
|
|
166 length(charlist),
|
|
167 "non-numeric value(s). To help you check your data, please find below a short overview:\n")
|
|
168 charover <- 1
|
|
169 while((length(err.stock)<10)&(length(err.stock)<(length(charlist)+1))){
|
|
170 charex <- paste("variable",colnames(DM)[ceiling(charlist[charover]/nrow(DM))],
|
|
171 "- sample",rownames(DM)[charlist[charover]-floor(charlist[charover]/nrow(DM))*nrow(DM)],
|
|
172 "~> value:",as.matrix(DM)[charlist[charover]],"\n")
|
|
173 err.stock <- c(err.stock,charex)
|
|
174 charover <- charover + 1
|
|
175 }
|
|
176 stop(c(err.stock,"The dataMatrix file is supposed to contain only numeric values.\n- - - - - - - - -\n"))
|
|
177 }
|
3
|
178
|
|
179 ### Initialization
|
|
180
|
|
181 DM <- cbind(c_class,DM)
|
|
182
|
|
183 stat.res <- t(DM[0,-1,drop=FALSE])
|
|
184 names <- NULL
|
|
185
|
|
186 mean.res <- NULL
|
|
187 mean.names <- NULL
|
0
|
188
|
3
|
189 sd.res <- NULL
|
|
190 sd.names <- NULL
|
|
191
|
|
192 med.res <- NULL
|
|
193 med.names <- NULL
|
|
194
|
|
195 quart.res <- NULL
|
|
196 quart.names <- NULL
|
|
197
|
|
198 dec.res <- NULL
|
|
199 dec.names <- NULL
|
|
200
|
|
201 NA.res <- NULL
|
|
202 NA.names <- NULL
|
|
203 pct_NA.res <- NULL
|
|
204 pct_NA.names <- NULL
|
|
205
|
|
206 fold.res <- NULL
|
|
207 fold.names <- NULL
|
|
208
|
|
209 if(("NA" %in% stat.list)||(test.fold=="Yes")){
|
|
210 graphs <- 1
|
|
211 }else{
|
|
212 graphs=0
|
|
213 }
|
|
214
|
|
215 data_bp <- data.frame() #table for NA barplot
|
|
216
|
0
|
217
|
1
|
218
|
3
|
219 ### Computation
|
0
|
220
|
1
|
221
|
3
|
222 for(j in 1:nb_class){
|
|
223
|
|
224 # Mean ---------
|
|
225
|
|
226 if("mean" %in% stat.list){
|
|
227 mean.res <- cbind(mean.res, colMeans(DM[which(DM$c_class==classnames[j]),-1],na.rm=TRUE))
|
|
228 mean.names <- cbind(mean.names, paste("Mean",classnames[j], sep="_"))
|
|
229 if(j == nb_class){
|
|
230 stat.res <- cbind(stat.res, mean.res)
|
|
231 names <- cbind(names, mean.names)
|
|
232 }
|
|
233 }
|
|
234
|
|
235 # Standard deviation -----
|
|
236
|
|
237 if("sd" %in% stat.list){
|
|
238 sd.res <- cbind(sd.res, apply(DM[which(DM$c_class==classnames[j]),-1],2,sd,na.rm=TRUE))
|
|
239 sd.names <- cbind(sd.names, paste("Sd",classnames[j], sep="_"))
|
|
240 if(j == nb_class){
|
|
241 stat.res <- cbind(stat.res, sd.res)
|
|
242 names <- cbind(names, sd.names)
|
|
243 }
|
|
244 }
|
|
245
|
|
246 # Median ---------
|
|
247
|
|
248 if(("median" %in% stat.list)&&(!("quartile" %in% stat.list))){
|
|
249 med.res <- cbind(med.res, apply(DM[which(DM$c_class==classnames[j]),-1],2,median,na.rm=TRUE))
|
|
250 med.names <- cbind(med.names, paste("Median",classnames[j], sep="_"))
|
|
251 if(j == nb_class){
|
|
252 stat.res <- cbind(stat.res, med.res)
|
|
253 names <- cbind(names, med.names)
|
|
254 }
|
|
255 }
|
|
256
|
|
257 # Quartiles ------
|
|
258
|
|
259 if("quartile" %in% stat.list){
|
|
260 quart.res <- cbind(quart.res,t(apply(DM[which(DM$c_class==classnames[j]),-1],2,quantile,na.rm=TRUE)))
|
|
261 quart.names <- cbind(quart.names, paste("Min",classnames[j], sep="_"),paste("Q1",classnames[j], sep="_"),
|
|
262 paste("Median",classnames[j],sep="_"),paste("Q3",classnames[j],sep="_"),
|
|
263 paste("Max",classnames[j],sep="_"))
|
|
264 if(j == nb_class){
|
|
265 stat.res <- cbind(stat.res, quart.res)
|
|
266 names <- cbind(names, quart.names)
|
|
267 }
|
|
268 }
|
|
269
|
|
270 # Deciles ------
|
|
271
|
|
272 if("decile" %in% stat.list){
|
|
273 dec.res <- cbind(dec.res,t(apply(DM[which(DM$c_class==classnames[j]),-1],2,quantile,na.rm=TRUE,seq(0,1,0.1))))
|
|
274 dec.names <- cbind(dec.names, t(matrix(paste((paste("D",seq(0,10,1),sep="")),classnames[j],sep="_"))))
|
|
275 if(j == nb_class){
|
|
276 stat.res <- cbind(stat.res, dec.res)
|
|
277 names <- cbind(names, dec.names)
|
|
278 }
|
|
279 }
|
|
280
|
|
281 # Missing values ------------
|
|
282
|
|
283 if("NA" %in% stat.list){
|
|
284
|
|
285 nb_NA <- apply(DM[which(DM$c_class==classnames[j]),-1],2,function(x) sum(is.na(x)))
|
|
286 pct_NA <- round(nb_NA/nrow(DM[which(DM$c_class==classnames[j]),-1])*100,digits=4)
|
|
287 NA.res <- cbind(NA.res,nb_NA)
|
|
288 pct_NA.res <- cbind(pct_NA.res,pct_NA)
|
|
289 NA.names <- cbind(NA.names, paste("NA",classnames[j], sep="_"))
|
|
290 pct_NA.names <- cbind(pct_NA.names,paste("Pct_NA", classnames[j], sep="_"))
|
|
291 if(j == nb_class){
|
|
292 stat.res <- cbind(stat.res, NA.res,pct_NA.res)
|
|
293 names <- cbind(names, NA.names,pct_NA.names)
|
|
294 }
|
0
|
295
|
3
|
296 #for barplots
|
|
297 Nb_NA_0_20 <- 0
|
|
298 Nb_NA_20_40 <- 0
|
|
299 Nb_NA_40_60 <- 0
|
|
300 Nb_NA_60_80 <- 0
|
|
301 Nb_NA_80_100 <- 0
|
|
302
|
|
303 for (i in 1:length(pct_NA)){
|
|
304
|
|
305 if ((0<=pct_NA[i])&(pct_NA[i]<20)){
|
|
306 Nb_NA_0_20=Nb_NA_0_20+1}
|
|
307
|
|
308 if ((20<=pct_NA[i])&(pct_NA[i]<40)){
|
|
309 Nb_NA_20_40=Nb_NA_20_40+1}
|
|
310
|
|
311 if ((40<=pct_NA[i])&(pct_NA[i]<60)){
|
|
312 Nb_NA_40_60=Nb_NA_40_60+1}
|
|
313
|
|
314 if ((60<=pct_NA[i])&(pct_NA[i]<80)){
|
|
315 Nb_NA_60_80=Nb_NA_60_80+1}
|
|
316
|
|
317 if ((80<=pct_NA[i])&(pct_NA[i]<=100)){
|
|
318 Nb_NA_80_100=Nb_NA_80_100+1}
|
|
319 }
|
|
320 data_bp[1,j] <- Nb_NA_0_20
|
|
321 data_bp[2,j] <- Nb_NA_20_40
|
|
322 data_bp[3,j] <- Nb_NA_40_60
|
|
323 data_bp[4,j] <- Nb_NA_60_80
|
|
324 data_bp[5,j] <- Nb_NA_80_100
|
|
325 rownames(data_bp) <- c("0%-20%", "20%-40%", "40%-60%", "60%-80%", "80%-100%")
|
|
326
|
|
327 if(j == nb_class){
|
|
328
|
|
329 # Alert message if there is no missing value in data matrix
|
|
330 sum_total <- sum(NA.res)
|
|
331 alerte <- NULL
|
|
332 if(sum_total==0){
|
|
333 alerte <- c(alerte, "Data Matrix contains no NA.\n")
|
|
334 }
|
|
335 if(length(alerte) != 0){
|
|
336 cat(alerte,"\n")
|
|
337 }
|
|
338
|
|
339
|
|
340 colnames(data_bp) <- classnames
|
|
341 data_bp <- as.matrix(data_bp)
|
|
342 }
|
|
343 }
|
|
344
|
|
345
|
|
346 # Mean fold change ------------
|
|
347
|
|
348 if(test.fold=="Yes"){
|
|
349 if(nb_class >= 2){
|
|
350 if(j!=nb_class){
|
|
351 ratio1 <- NULL
|
|
352 ratio2 <- NULL
|
|
353 if(method=="each_class"){
|
|
354 fold.frac <- "Top"
|
|
355 }
|
|
356 for(k in (j+1):nb_class) {
|
|
357 if(fold.frac=="Bottom"){
|
|
358 ratio1 <- classnames[k]
|
|
359 ratio2 <- classnames[j]
|
|
360 }else{
|
|
361 ratio1 <- classnames[j]
|
|
362 ratio2 <- classnames[k]
|
|
363 }
|
|
364 fold <- colMeans(DM[which(DM$c_class==ratio1),-1],na.rm=TRUE)/
|
|
365 colMeans(DM[which(DM$c_class==ratio2),-1],na.rm=TRUE)
|
|
366 if(logarithm=="log2"){
|
|
367 fold.res <- cbind(fold.res,log2(fold))
|
|
368 }else if(logarithm=="log10"){
|
|
369 fold.res <- cbind(fold.res,log10(fold))
|
|
370 }else{
|
|
371 fold.res <- cbind(fold.res, fold)
|
|
372 }
|
|
373 if(logarithm == "none"){
|
|
374 fold.names <- cbind(fold.names,paste("fold",ratio1,"VS", ratio2, sep="_"))
|
|
375 }else{
|
|
376 fold.names <- cbind(fold.names,paste(logarithm, "fold", ratio1, "VS", ratio2, sep="_"))
|
|
377 }
|
|
378 }
|
|
379
|
|
380 }else{
|
|
381 stat.res <- cbind(stat.res,fold.res)
|
|
382 names <- cbind(names, fold.names)
|
|
383 }
|
|
384 }
|
|
385 }
|
|
386
|
|
387 }
|
|
388
|
|
389 ############
|
|
390
|
|
391
|
|
392 # check columns names in variableMetadata
|
|
393
|
|
394 VM.names <- colnames(VM)
|
|
395 for (i in 1:length(VM.names)){
|
|
396 for (j in 1:length(names)){
|
|
397 if (VM.names[i]==names[j]){
|
|
398 names[j] <- paste(names[j], "2", sep="_")
|
0
|
399 }
|
|
400 }
|
|
401 }
|
3
|
402
|
|
403 colnames(stat.res) <- names
|
|
404
|
0
|
405
|
1
|
406
|
|
407
|
|
408
|
0
|
409
|
1
|
410 # Output ---------------------------------------------------------------------------------------------------
|
0
|
411
|
3
|
412 VM <-cbind(VM,stat.res)
|
0
|
413
|
|
414 write.table(VM, VM.output,sep="\t", quote=FALSE, row.names=FALSE)
|
|
415
|
3
|
416
|
|
417 ### graphics pdf
|
|
418
|
|
419 if(graphs == 1){
|
0
|
420
|
|
421 pdf(graphs.output)
|
3
|
422
|
0
|
423
|
1
|
424 #Barplots for NA
|
3
|
425 if("NA" %in% stat.list){
|
|
426 graph.colors <- c("green3","palegreen3","lightblue","orangered","red")
|
0
|
427 par(mar=c(5.1, 4.1, 4.1, 8.1), xpd=TRUE)
|
3
|
428
|
|
429 bp=barplot(data_bp, col=graph.colors, main="Proportion of NA", xlab="Classes", ylab="Variables")
|
|
430 legend("topright", fill=graph.colors,rownames(data_bp), inset=c(-0.3,0))
|
|
431
|
0
|
432 stock=0
|
|
433 for (i in 1:nrow(data_bp)){
|
1
|
434 text(bp, stock+data_bp[i,]/2, data_bp[i,], col="white", cex=0.7)
|
0
|
435 stock <- stock+data_bp[i,]
|
|
436 }
|
|
437
|
3
|
438 }
|
|
439
|
|
440 # Boxplots for fold test
|
1
|
441
|
3
|
442 if((test.fold=="Yes")&&(nb_class >= 2)){
|
|
443
|
|
444 clean_fold <- fold.res
|
1
|
445 for(i in 1:nrow(clean_fold)){
|
|
446 for(j in 1:ncol(clean_fold)){
|
|
447 if(is.infinite(clean_fold[i,j])){
|
|
448 clean_fold[i,j] <- NA
|
|
449 }
|
|
450 }
|
|
451 }
|
|
452 for (j in 1:ncol(clean_fold)){
|
|
453 title <- paste(fold.names[j])
|
3
|
454 boxplot(clean_fold[,j], main=title)
|
1
|
455 }
|
3
|
456 }
|
|
457
|
|
458 dev.off()
|
0
|
459
|
3
|
460 }else{
|
|
461 pdf(graphs.output)
|
|
462 plot.new()
|
|
463 legend("center","You did not select any option with graphical output.")
|
|
464 dev.off()
|
|
465 }
|
0
|
466
|
3
|
467 }
|
|
468
|
|
469
|
|
470
|
|
471
|
0
|
472
|
|
473 |