Mercurial > repos > labis-app > galaxy_proteomics

comparison t-test.R @ 0:ba070efb6f78 draft

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
planemo upload commit 13e72e84c523bda22bda792bbebf4720d28542d5-dirty
author labis-app
date Tue, 03 Jul 2018 17:34:13 -0400
parents
children
comparison
equal deleted inserted replaced
-1:000000000000 0:ba070efb6f78
1 #!/usr/bin/env Rscript
2
3 # t-test.R
4 # AUTHOR: Daniel Travieso
5 # E-mail: danielgtravieso@gmail.com
6 # LAST REVISED: April 2015
7 #
8 # Required packages to work: (getopt", "gtools")
9 # Laboratory of Mass Spectrometry at Brazilian Biosciences National Laboratory
10 # http://lnbio.cnpem.br/
11 # Copyright CC BY-NC-SA (c) 2014 Brazilian Center for Research in Energy and Materials
12 # All rights reserved.
13 require('gtools', quietly=TRUE);
14 require('getopt', quietly=TRUE);
15 #include and execute the read util script
16 library('read_util.R');
17 library('write_util.R');
18
19 #define de options input that the read_util$code will have
20 opt = matrix(c(
21 'inputfile_name', 'i', 1, 'character',
22 'type', 't', 1, 'character',
23 'outputfile_name', 'o', 1, 'character'
24 ),byrow=TRUE, ncol=4);
25
26 # parse de input
27 options = getopt(opt);
28
29 read_util <- read_function(options);
30
31 i<-1;
32 columns <- list();
33 aux <- c();
34 for (cat in read_util$diff_cat) {
35 col <- read_util$col_names[gsub(read_util$regex, "\\1", read_util$col_names) == cat]
36 aux <- c(aux, col);
37 columns[[i]] <- col;
38 i<-i+1;
39 }
40 # this is a filtered read_util$table to help with calculations
41 table_only_columns <- read_util$table[-1, aux]
42
43 # this loop computes the ttest result for each row
44 # and adds it to a vector
45 i <- 2;
46 ttestresult <- c("");
47 ttestsignificant <- c("");
48 if (length(read_util$diff_cat) < 2) {
49 print(sprintf("Can't calculate t-test. There is only one category for %s collumns", read_util$code));
50 q(1,save="no");
51 }
52
53 for (i in seq(2, nrow(table_only_columns)+1)) {
54 # the t-test arguments are the control values vector, the treatment values vector
55 # and some extra arguments. var.equal says it's a student t-test with stardard
56 # deviations assumed equal. mu=0 sets the hipothesis to be null.
57 ttestresult[i] <- t.test(table_only_columns[i-1, columns[[1]]],
58 table_only_columns[i-1, columns[[2]]], var.equal=TRUE, mu=0)$p.value;
59 if (is.na(ttestresult[i]))
60 ttestresult[i] = 1.0
61 }
62
63 # this defines if the p-value returned for each row is significant
64 ttestsignificant[ttestresult <= 0.05] <- "+"
65 ttestsignificant[ttestresult > 0.05] <- ""
66
67
68 # create two extra rows on the read_util$table, one for p-values and other
69 # for siginificance
70 #TODO: ou colocar perto da intensidade que se refere ou na 3ยช coluna
71 read_util$table[paste0("T.test.result.", read_util$code)] <- NA;
72 read_util$table[paste0("T.test.result.", read_util$code)] <- ttestresult;
73 read_util$table[paste0("T.test.significant.", read_util$code)] <- NA;
74 read_util$table[paste0("T.test.significant.", read_util$code)] <- ttestsignificant;
75
76
77
78
79 # write out the read_util$table
80 writeout(options$outputfile_name, read_util$table);