Mercurial > repos > marie-tremblay-metatoul > asca

comparison asca_w4m.R @ 0:93312041f1d5 draft default tip

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
planemo upload for repository https://github.com/workflow4metabolomics/ascaw4m commit 7ea9b0f8abc5a60c2c04fd2098788497f14766b6
author marie-tremblay-metatoul
date Fri, 21 Sep 2018 05:51:14 -0400
parents
children
comparison
equal deleted inserted replaced
-1:000000000000 0:93312041f1d5
1 asca_w4m <- function(datamatrix, samplemetadata, factors, variablemetadata, threshold, scaling="none", nPerm)
2 {
3 ## Transpose
4 # datamatrix <- t(datamatrix)
5
6 # Check sample ID's
7 rownames(datamatrix) <- make.names(rownames(datamatrix), unique = TRUE)
8 colnames(datamatrix) <- make.names(colnames(datamatrix), unique = TRUE)
9 rownames(samplemetadata) <- make.names(rownames(samplemetadata), unique = TRUE)
10 rownames(variablemetadata) <- make.names(rownames(variablemetadata), unique = TRUE)
11
12 if(!identical(rownames(datamatrix), rownames(samplemetadata)))
13 {
14 if(identical(sort(rownames(datamatrix)), sort(rownames(samplemetadata))))
15 {
16 cat("\n\nMessage: Re-ordering dataMatrix sample names to match sampleMetadata\n")
17 datamatrix <- datamatrix[rownames(samplemetadata), , drop = FALSE]
18 stopifnot(identical(sort(rownames(datamatrix)), sort(rownames(samplemetadata))))
19 }else {
20
21 cat("\n\nStop: The sample names of dataMatrix and sampleMetadata do not match:\n")
22 print(cbind.data.frame(indice = 1:nrow(datamatrix),
23 dataMatrix=rownames(datamatrix),
24 sampleMetadata=rownames(samplemetadata))[rownames(datamatrix) != rownames(samplemetadata), , drop = FALSE])
25 }
26 }
27
28 # Check feature ID's
29 if(!identical(colnames(datamatrix), rownames(variablemetadata)))
30 {
31 if(identical(sort(colnames(datamatrix)), sort(rownames(variablemetadata))))
32 {
33 cat("\n\nMessage: Re-ordering dataMatrix variable names to match variableMetadata\n")
34 datamatrix <- datamatrix[, rownames(variablemetadata), drop = FALSE]
35 stopifnot(identical(sort(colnames(datamatrix)), sort(rownames(variablemetadata))))
36 }else {
37 cat("\n\nStop: The variable names of dataMatrix and variableMetadata do not match:\n")
38 print(cbind.data.frame(indice = 1:ncol(datamatrix),
39 dataMatrix=colnames(datamatrix),
40 variableMetadata=rownames(variablemetadata))[colnames(datamatrix) != rownames(variablemetadata), , drop = FALSE])
41 }
42 }
43
44 # Design
45 design <- data.matrix(samplemetadata[, colnames(samplemetadata) %in% factors])
46
47 # Scaling if scaling!=none
48 datamatrix <- prep(datamatrix, scaling)
49
50 # Computation of the A-SCA model
51 data.asca <- ASCA.Calculate_w4m(datamatrix, design, scaling=scaling)
52
53 # Permutation test
54 data.asca.permutation <- ASCA.DoPermutationTest(data.asca[[1]], perm=nPerm)
55 p <- c(data.asca.permutation, 0)
56
57
58 # % of explained variance
59 ssq <- (data.asca[[2]]$summary.ssq)
60 ssq <- cbind(round(rbind(ssq[2], ssq[3],ssq[4],ssq[5])*100, 2), p)
61 rownames(ssq) <- c(factors[1], factors[2], "Interaction", "Residuals")
62 colnames(ssq) <- c("% of explained variance", "Permutation p-value")
63
64 # Add Scores and loadings at the end of meatadata files
65 noms <- colnames(samplemetadata)
66 samplemetadata <- cbind(samplemetadata, (data.asca[[1]]$'1'$means.matrix + data.asca[[1]]$remainder) %*% data.asca[[1]]$'1'$svd$v[, 1:2],
67 (data.asca[[1]]$'2'$means.matrix + data.asca[[1]]$remainder) %*% data.asca[[1]]$'2'$svd$v[, 1:2],
68 (data.asca[[1]]$'12'$means.matrix + data.asca[[1]]$remainder) %*% data.asca[[1]]$'12'$svd$v[, 1:2])
69 colnames(samplemetadata) <- c(noms, paste(factors[1],"XSCOR-p1", sep="_"), paste(factors[1],"XSCOR-p2", sep="_"),
70 paste(factors[2],"XSCOR-p1", sep="_"), paste(factors[2],"XSCOR-p2", sep="_"),
71 "Interact_XSCOR-p1", "Interact_XSCOR-p2")
72
73 noms <- colnames(variablemetadata)
74 variablemetadata <- cbind(variablemetadata, data.asca[[1]]$'1'$svd$v[, 1:2], data.asca[[1]]$'2'$svd$v[, 1:2], data.asca[[1]]$'12'$svd$v[, 1:2])
75 colnames(variablemetadata) <- c(noms, paste(factors[1],"XLOAD-p1", sep="_"), paste(factors[1],"XLOAD-p2", sep="_"),
76 paste(factors[2],"XLOAD-p1", sep="_"), paste(factors[2],"XLOAD-p2", sep="_"),
77 "Interact_XLOAD-p1", "Interact_XLOAD-p2")
78
79 l <- list(data.asca[[1]], data.asca.permutation, ssq, samplemetadata, variablemetadata)
80 names(l) <- c("ASCA","p-values", "ssq", "samplemetadata", "variablemetadata")
81 return(l)
82 }