Mercurial > repos > mmonsoor > probmetab

annotate export.class.table-color-graph.R @ 1:c4bd0c40eb3b draft

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
planemo upload commit 1c73dd5ed7b5b8595963e0e1b6f8427978d8d5f5-dirty
author lecorguille
date Mon, 04 Jul 2016 09:51:51 -0400
parents e13ec2c3fabe
children
Ignore whitespace changes - Everywhere: Within whitespace: At end of lines:
rev   line source
0
e13ec2c3fabe planemo upload commit 25fd6a739741295e3f434e0be0286dee61e06825
mmonsoor
parents:
diff changeset
1 #' export.class.table
e13ec2c3fabe planemo upload commit 25fd6a739741295e3f434e0be0286dee61e06825
mmonsoor
parents:
diff changeset
2 #'
e13ec2c3fabe planemo upload commit 25fd6a739741295e3f434e0be0286dee61e06825
mmonsoor
parents:
diff changeset
3 #' Builds a matrix with the probability for all mass to candidate compounds
e13ec2c3fabe planemo upload commit 25fd6a739741295e3f434e0be0286dee61e06825
mmonsoor
parents:
diff changeset
4 #' assignments, by averaging the number of assignments obtained by the gibbs sampler algorithm
e13ec2c3fabe planemo upload commit 25fd6a739741295e3f434e0be0286dee61e06825
mmonsoor
parents:
diff changeset
5 #' or ordering the compound candidates with the likelihood matrix.
e13ec2c3fabe planemo upload commit 25fd6a739741295e3f434e0be0286dee61e06825
mmonsoor
parents:
diff changeset
6 #'
e13ec2c3fabe planemo upload commit 25fd6a739741295e3f434e0be0286dee61e06825
mmonsoor
parents:
diff changeset
7 #' @param gibbsL a list of attributions and probabilities from gibbs.samp function.
e13ec2c3fabe planemo upload commit 25fd6a739741295e3f434e0be0286dee61e06825
mmonsoor
parents:
diff changeset
8 #' @param reactionM data.frame with compound annotation information.
e13ec2c3fabe planemo upload commit 25fd6a739741295e3f434e0be0286dee61e06825
mmonsoor
parents:
diff changeset
9 #' @param molIon non redundant ion annotation.
e13ec2c3fabe planemo upload commit 25fd6a739741295e3f434e0be0286dee61e06825
mmonsoor
parents:
diff changeset
10 #' @param probM optionally to gibbsL, a matrix of likelihoods.
e13ec2c3fabe planemo upload commit 25fd6a739741295e3f434e0be0286dee61e06825
mmonsoor
parents:
diff changeset
11 #' @param html logical, indicating whether a html file should be generated. This parameter uses the raw data to plot EICs and may be time consuming.
e13ec2c3fabe planemo upload commit 25fd6a739741295e3f434e0be0286dee61e06825
mmonsoor
parents:
diff changeset
12 #' @param filename html file name, the default is "test".
e13ec2c3fabe planemo upload commit 25fd6a739741295e3f434e0be0286dee61e06825
mmonsoor
parents:
diff changeset
13 #' @param burnIn how many samples of the gibbs sampler should be discarded.
e13ec2c3fabe planemo upload commit 25fd6a739741295e3f434e0be0286dee61e06825
mmonsoor
parents:
diff changeset
14 #' @param linkPattern which pattern should be linked to compound id, for now we have
e13ec2c3fabe planemo upload commit 25fd6a739741295e3f434e0be0286dee61e06825
mmonsoor
parents:
diff changeset
15 #' implemented "kegg", "pubchem" and "chebi" patterns.
e13ec2c3fabe planemo upload commit 25fd6a739741295e3f434e0be0286dee61e06825
mmonsoor
parents:
diff changeset
16 #' @param m.test statistical test to compare mean differences. This option
e13ec2c3fabe planemo upload commit 25fd6a739741295e3f434e0be0286dee61e06825
mmonsoor
parents:
diff changeset
17 #' is only available to single acquisition mode analysis, with options
e13ec2c3fabe planemo upload commit 25fd6a739741295e3f434e0be0286dee61e06825
mmonsoor
parents:
diff changeset
18 #' "t.test" and "anova".
e13ec2c3fabe planemo upload commit 25fd6a739741295e3f434e0be0286dee61e06825
mmonsoor
parents:
diff changeset
19 #' @param class1 if the m.test is "t.test" first class to compare in the test,
e13ec2c3fabe planemo upload commit 25fd6a739741295e3f434e0be0286dee61e06825
mmonsoor
parents:
diff changeset
20 #' according with xcmsSet phenoData.
e13ec2c3fabe planemo upload commit 25fd6a739741295e3f434e0be0286dee61e06825
mmonsoor
parents:
diff changeset
21 #' @param class2 if the m.test is "t.test" second class to compare in the test,
e13ec2c3fabe planemo upload commit 25fd6a739741295e3f434e0be0286dee61e06825
mmonsoor
parents:
diff changeset
22 #' according with xcmsSet phenoData.
e13ec2c3fabe planemo upload commit 25fd6a739741295e3f434e0be0286dee61e06825
mmonsoor
parents:
diff changeset
23 #' @param norm logical, if TRUE performs median normalization from (Anal. Chem. 2011, 83, 5864-5872).
e13ec2c3fabe planemo upload commit 25fd6a739741295e3f434e0be0286dee61e06825
mmonsoor
parents:
diff changeset
24 #' @param DB data.frame table used to search compounds, with the field name to be incorporated in the html table.
e13ec2c3fabe planemo upload commit 25fd6a739741295e3f434e0be0286dee61e06825
mmonsoor
parents:
diff changeset
25 #' @param prob how to calculate the probability to attribute a mass to a compound.
e13ec2c3fabe planemo upload commit 25fd6a739741295e3f434e0be0286dee61e06825
mmonsoor
parents:
diff changeset
26 #' Default is "count", which divide the number of times each identity was
e13ec2c3fabe planemo upload commit 25fd6a739741295e3f434e0be0286dee61e06825
mmonsoor
parents:
diff changeset
27 #' was attributed by the number of samples. Optionally the user could
e13ec2c3fabe planemo upload commit 25fd6a739741295e3f434e0be0286dee61e06825
mmonsoor
parents:
diff changeset
28 #' choose to use the mean of the probabilities of the identity, "mean".
e13ec2c3fabe planemo upload commit 25fd6a739741295e3f434e0be0286dee61e06825
mmonsoor
parents:
diff changeset
29 #' @return A list with a matrix "classTable" with attributions and probabilities and
e13ec2c3fabe planemo upload commit 25fd6a739741295e3f434e0be0286dee61e06825
mmonsoor
parents:
diff changeset
30 #' indexes of selected masses from xcms peak table.
e13ec2c3fabe planemo upload commit 25fd6a739741295e3f434e0be0286dee61e06825
mmonsoor
parents:
diff changeset
31 #'
e13ec2c3fabe planemo upload commit 25fd6a739741295e3f434e0be0286dee61e06825
mmonsoor
parents:
diff changeset
32 #' @export
e13ec2c3fabe planemo upload commit 25fd6a739741295e3f434e0be0286dee61e06825
mmonsoor
parents:
diff changeset
33
e13ec2c3fabe planemo upload commit 25fd6a739741295e3f434e0be0286dee61e06825
mmonsoor
parents:
diff changeset
34 export.class.table <- function(gibbsL=NULL, reactionM, molIon=NULL, probM=NULL, html=FALSE, filename="test", burnIn=3000, linkPattern="kegg", m.test="none", class1=NULL, class2=NULL, norm=FALSE, DB, prob="count", colorplot=FALSE, addLink=NULL) {
e13ec2c3fabe planemo upload commit 25fd6a739741295e3f434e0be0286dee61e06825
mmonsoor
parents:
diff changeset
35
e13ec2c3fabe planemo upload commit 25fd6a739741295e3f434e0be0286dee61e06825
mmonsoor
parents:
diff changeset
36 plotEIC <- function (xcmsObject, figidx, pngidx, colorplot, mode=NULL) {
e13ec2c3fabe planemo upload commit 25fd6a739741295e3f434e0be0286dee61e06825
mmonsoor
parents:
diff changeset
37 dir.create(paste(filename,"_fig",sep=""))
e13ec2c3fabe planemo upload commit 25fd6a739741295e3f434e0be0286dee61e06825
mmonsoor
parents:
diff changeset
38 gt<-groups(xcmsObject)
e13ec2c3fabe planemo upload commit 25fd6a739741295e3f434e0be0286dee61e06825
mmonsoor
parents:
diff changeset
39 if(colorplot==TRUE){
e13ec2c3fabe planemo upload commit 25fd6a739741295e3f434e0be0286dee61e06825
mmonsoor
parents:
diff changeset
40 gt2 <- gt[figidx,]
e13ec2c3fabe planemo upload commit 25fd6a739741295e3f434e0be0286dee61e06825
mmonsoor
parents:
diff changeset
41 rgt <- gt2[,c("rtmin","rtmax")]
e13ec2c3fabe planemo upload commit 25fd6a739741295e3f434e0be0286dee61e06825
mmonsoor
parents:
diff changeset
42 rgt[,1] <- rgt[,1]-100
e13ec2c3fabe planemo upload commit 25fd6a739741295e3f434e0be0286dee61e06825
mmonsoor
parents:
diff changeset
43 rgt[,2] <- rgt[,2]+100
e13ec2c3fabe planemo upload commit 25fd6a739741295e3f434e0be0286dee61e06825
mmonsoor
parents:
diff changeset
44 #require(doMC)
e13ec2c3fabe planemo upload commit 25fd6a739741295e3f434e0be0286dee61e06825
mmonsoor
parents:
diff changeset
45 #registerDoMC()
e13ec2c3fabe planemo upload commit 25fd6a739741295e3f434e0be0286dee61e06825
mmonsoor
parents:
diff changeset
46 #system.time(
e13ec2c3fabe planemo upload commit 25fd6a739741295e3f434e0be0286dee61e06825
mmonsoor
parents:
diff changeset
47 #foreach(i=1:nrow(gt2)) %dopar% {
e13ec2c3fabe planemo upload commit 25fd6a739741295e3f434e0be0286dee61e06825
mmonsoor
parents:
diff changeset
48 for(i in 1:nrow(gt2)){
e13ec2c3fabe planemo upload commit 25fd6a739741295e3f434e0be0286dee61e06825
mmonsoor
parents:
diff changeset
49 groupidx1 <- which(gt[,"rtmed"] > rgt[i,1] & gt[,"rtmed"] < rgt[i,2] & gt[,"mzmed"]> gt2[i,"mzmed"]-0.1 & gt[,"mzmed"]< gt2[i,"mzmed"]+0.1)
e13ec2c3fabe planemo upload commit 25fd6a739741295e3f434e0be0286dee61e06825
mmonsoor
parents:
diff changeset
50 eiccor <- getEIC(xcmsObject, groupidx = groupidx1)
e13ec2c3fabe planemo upload commit 25fd6a739741295e3f434e0be0286dee61e06825
mmonsoor
parents:
diff changeset
51 png(paste(filename, "_fig/", sprintf("%003d", i), ".png", sep=""))
e13ec2c3fabe planemo upload commit 25fd6a739741295e3f434e0be0286dee61e06825
mmonsoor
parents:
diff changeset
52 plot(eiccor, xcmsObject, groupidx = 1)
e13ec2c3fabe planemo upload commit 25fd6a739741295e3f434e0be0286dee61e06825
mmonsoor
parents:
diff changeset
53 dev.off()
e13ec2c3fabe planemo upload commit 25fd6a739741295e3f434e0be0286dee61e06825
mmonsoor
parents:
diff changeset
54 }
e13ec2c3fabe planemo upload commit 25fd6a739741295e3f434e0be0286dee61e06825
mmonsoor
parents:
diff changeset
55 } else {
e13ec2c3fabe planemo upload commit 25fd6a739741295e3f434e0be0286dee61e06825
mmonsoor
parents:
diff changeset
56 gt <- gt[figidx,]
e13ec2c3fabe planemo upload commit 25fd6a739741295e3f434e0be0286dee61e06825
mmonsoor
parents:
diff changeset
57 rgt <- gt[,c("rtmin","rtmax")]
e13ec2c3fabe planemo upload commit 25fd6a739741295e3f434e0be0286dee61e06825
mmonsoor
parents:
diff changeset
58 rgt[,1] <- rgt[,1]-100
e13ec2c3fabe planemo upload commit 25fd6a739741295e3f434e0be0286dee61e06825
mmonsoor
parents:
diff changeset
59 rgt[,2] <- rgt[,2]+100
e13ec2c3fabe planemo upload commit 25fd6a739741295e3f434e0be0286dee61e06825
mmonsoor
parents:
diff changeset
60
e13ec2c3fabe planemo upload commit 25fd6a739741295e3f434e0be0286dee61e06825
mmonsoor
parents:
diff changeset
61 eics <- getEIC(xcmsObject, mzrange=gt, rtrange =rgt, groupidx = 1:nrow(gt))
e13ec2c3fabe planemo upload commit 25fd6a739741295e3f434e0be0286dee61e06825
mmonsoor
parents:
diff changeset
62 png(file.path(paste(filename, "_fig/%003d.png", sep="")), height=768, width=1024)
e13ec2c3fabe planemo upload commit 25fd6a739741295e3f434e0be0286dee61e06825
mmonsoor
parents:
diff changeset
63 #png(file.path(paste(filename, "_fig/", pngidx, sep="")), h=768, w=1024)
e13ec2c3fabe planemo upload commit 25fd6a739741295e3f434e0be0286dee61e06825
mmonsoor
parents:
diff changeset
64 plot(eics, xcmsObject)
e13ec2c3fabe planemo upload commit 25fd6a739741295e3f434e0be0286dee61e06825
mmonsoor
parents:
diff changeset
65 dev.off()
e13ec2c3fabe planemo upload commit 25fd6a739741295e3f434e0be0286dee61e06825
mmonsoor
parents:
diff changeset
66 }
e13ec2c3fabe planemo upload commit 25fd6a739741295e3f434e0be0286dee61e06825
mmonsoor
parents:
diff changeset
67 if(!is.null(mode)) {
e13ec2c3fabe planemo upload commit 25fd6a739741295e3f434e0be0286dee61e06825
mmonsoor
parents:
diff changeset
68 pngs <- dir(paste(filename, "_fig/", sep=""))
e13ec2c3fabe planemo upload commit 25fd6a739741295e3f434e0be0286dee61e06825
mmonsoor
parents:
diff changeset
69 if(length(grep("pos|neg" , pngs))) pngs <- pngs[-grep("pos|neg" , pngs)]
e13ec2c3fabe planemo upload commit 25fd6a739741295e3f434e0be0286dee61e06825
mmonsoor
parents:
diff changeset
70 opng <- as.numeric(sub(".png","", pngs))
e13ec2c3fabe planemo upload commit 25fd6a739741295e3f434e0be0286dee61e06825
mmonsoor
parents:
diff changeset
71 pngs <- pngs[order(opng)]
e13ec2c3fabe planemo upload commit 25fd6a739741295e3f434e0be0286dee61e06825
mmonsoor
parents:
diff changeset
72 name1 <- paste(filename, "_fig/",pngs, sep="")
e13ec2c3fabe planemo upload commit 25fd6a739741295e3f434e0be0286dee61e06825
mmonsoor
parents:
diff changeset
73 name2 <- paste(filename, "_fig/",pngidx, mode, ".png", sep="")
e13ec2c3fabe planemo upload commit 25fd6a739741295e3f434e0be0286dee61e06825
mmonsoor
parents:
diff changeset
74 for(i in 1:length(name1)) file.rename(name1[i], name2[i])
e13ec2c3fabe planemo upload commit 25fd6a739741295e3f434e0be0286dee61e06825
mmonsoor
parents:
diff changeset
75 }
e13ec2c3fabe planemo upload commit 25fd6a739741295e3f434e0be0286dee61e06825
mmonsoor
parents:
diff changeset
76
e13ec2c3fabe planemo upload commit 25fd6a739741295e3f434e0be0286dee61e06825
mmonsoor
parents:
diff changeset
77 }
e13ec2c3fabe planemo upload commit 25fd6a739741295e3f434e0be0286dee61e06825
mmonsoor
parents:
diff changeset
78 allion <- molIon$molIon[molIon$molIon[,"isotope"]==0,]
e13ec2c3fabe planemo upload commit 25fd6a739741295e3f434e0be0286dee61e06825
mmonsoor
parents:
diff changeset
79 ReactMatrix <- reactionM[reactionM[,5]!="unknown",]
e13ec2c3fabe planemo upload commit 25fd6a739741295e3f434e0be0286dee61e06825
mmonsoor
parents:
diff changeset
80 x <- apply(unique(ReactMatrix[,c(2, 3)]), 2, as.numeric) # Have to look for all pairs
e13ec2c3fabe planemo upload commit 25fd6a739741295e3f434e0be0286dee61e06825
mmonsoor
parents:
diff changeset
81 y <- as.numeric(ReactMatrix[,4])
e13ec2c3fabe planemo upload commit 25fd6a739741295e3f434e0be0286dee61e06825
mmonsoor
parents:
diff changeset
82 prob_mean_ma <- matrix(0, nrow = length(y), ncol = nrow(x))
e13ec2c3fabe planemo upload commit 25fd6a739741295e3f434e0be0286dee61e06825
mmonsoor
parents:
diff changeset
83 # z_average <- matrix(0, nrow = length(y), ncol = length(x))
e13ec2c3fabe planemo upload commit 25fd6a739741295e3f434e0be0286dee61e06825
mmonsoor
parents:
diff changeset
84
e13ec2c3fabe planemo upload commit 25fd6a739741295e3f434e0be0286dee61e06825
mmonsoor
parents:
diff changeset
85 if (!is.null(gibbsL)){
e13ec2c3fabe planemo upload commit 25fd6a739741295e3f434e0be0286dee61e06825
mmonsoor
parents:
diff changeset
86 prob_table <- gibbsL$prob_table[,-c(1:burnIn)]
e13ec2c3fabe planemo upload commit 25fd6a739741295e3f434e0be0286dee61e06825
mmonsoor
parents:
diff changeset
87 class_table <- gibbsL$class_table[,-c(1:burnIn)]
e13ec2c3fabe planemo upload commit 25fd6a739741295e3f434e0be0286dee61e06825
mmonsoor
parents:
diff changeset
88 #indList <- tapply(1:nrow(ReactMatrix), as.numeric(ReactMatrix[,1]), function(x) x)
e13ec2c3fabe planemo upload commit 25fd6a739741295e3f434e0be0286dee61e06825
mmonsoor
parents:
diff changeset
89 coords <- tapply(1:nrow(ReactMatrix), ReactMatrix[,"molIonID"], function(x) x)
e13ec2c3fabe planemo upload commit 25fd6a739741295e3f434e0be0286dee61e06825
mmonsoor
parents:
diff changeset
90 coords2 <- unlist(lapply(coords, function(x) rep(x[1], length(x))))
e13ec2c3fabe planemo upload commit 25fd6a739741295e3f434e0be0286dee61e06825
mmonsoor
parents:
diff changeset
91 indList <- coords[order(unique(coords2))]
e13ec2c3fabe planemo upload commit 25fd6a739741295e3f434e0be0286dee61e06825
mmonsoor
parents:
diff changeset
92 fillMatrix <- function(j,i) {
e13ec2c3fabe planemo upload commit 25fd6a739741295e3f434e0be0286dee61e06825
mmonsoor
parents:
diff changeset
93 idp <- which(class_table[i,] == j)
e13ec2c3fabe planemo upload commit 25fd6a739741295e3f434e0be0286dee61e06825
mmonsoor
parents:
diff changeset
94 if(prob=="count") prob_mean_ma[j,i] <<- length(idp)/ncol(class_table)
e13ec2c3fabe planemo upload commit 25fd6a739741295e3f434e0be0286dee61e06825
mmonsoor
parents:
diff changeset
95 if(prob=="mean") prob_mean_ma[j,i] <<- mean(prob_table[i,idp])
e13ec2c3fabe planemo upload commit 25fd6a739741295e3f434e0be0286dee61e06825
mmonsoor
parents:
diff changeset
96 }
e13ec2c3fabe planemo upload commit 25fd6a739741295e3f434e0be0286dee61e06825
mmonsoor
parents:
diff changeset
97
e13ec2c3fabe planemo upload commit 25fd6a739741295e3f434e0be0286dee61e06825
mmonsoor
parents:
diff changeset
98
e13ec2c3fabe planemo upload commit 25fd6a739741295e3f434e0be0286dee61e06825
mmonsoor
parents:
diff changeset
99 for ( i in 1:nrow(x) ) {
e13ec2c3fabe planemo upload commit 25fd6a739741295e3f434e0be0286dee61e06825
mmonsoor
parents:
diff changeset
100
e13ec2c3fabe planemo upload commit 25fd6a739741295e3f434e0be0286dee61e06825
mmonsoor
parents:
diff changeset
101 sapply(indList[[i]], fillMatrix, i)
e13ec2c3fabe planemo upload commit 25fd6a739741295e3f434e0be0286dee61e06825
mmonsoor
parents:
diff changeset
102 }
e13ec2c3fabe planemo upload commit 25fd6a739741295e3f434e0be0286dee61e06825
mmonsoor
parents:
diff changeset
103 if(sum(prob_mean_ma=="NaN")) prob_mean_ma[prob_mean_ma=="NaN"] <- 0
e13ec2c3fabe planemo upload commit 25fd6a739741295e3f434e0be0286dee61e06825
mmonsoor
parents:
diff changeset
104 # for ( i in 1:nrow(x) ) {
e13ec2c3fabe planemo upload commit 25fd6a739741295e3f434e0be0286dee61e06825
mmonsoor
parents:
diff changeset
105 # for ( j in 1:length(y) ) {
e13ec2c3fabe planemo upload commit 25fd6a739741295e3f434e0be0286dee61e06825
mmonsoor
parents:
diff changeset
106 # idp <- which(class_table[i,] == j)
e13ec2c3fabe planemo upload commit 25fd6a739741295e3f434e0be0286dee61e06825
mmonsoor
parents:
diff changeset
107 # prob_mean_ma[j,i] <- mean(prob_table[i,idp])
e13ec2c3fabe planemo upload commit 25fd6a739741295e3f434e0be0286dee61e06825
mmonsoor
parents:
diff changeset
108 # # this is an alternative way to calculate the probabilities, should try latter, and compare results
e13ec2c3fabe planemo upload commit 25fd6a739741295e3f434e0be0286dee61e06825
mmonsoor
parents:
diff changeset
109 # #prob_mean_ma[j,i] <- length(idp)/ncol(class_table)
e13ec2c3fabe planemo upload commit 25fd6a739741295e3f434e0be0286dee61e06825
mmonsoor
parents:
diff changeset
110 # if ( prob_mean_ma[j,i] == "NaN" ) prob_mean_ma[j,i] <- 0
e13ec2c3fabe planemo upload commit 25fd6a739741295e3f434e0be0286dee61e06825
mmonsoor
parents:
diff changeset
111 # }
e13ec2c3fabe planemo upload commit 25fd6a739741295e3f434e0be0286dee61e06825
mmonsoor
parents:
diff changeset
112 # # do I still need this matrix?
e13ec2c3fabe planemo upload commit 25fd6a739741295e3f434e0be0286dee61e06825
mmonsoor
parents:
diff changeset
113 # k=which(prob_mean_ma[,i]==max(prob_mean_ma[,i]))
e13ec2c3fabe planemo upload commit 25fd6a739741295e3f434e0be0286dee61e06825
mmonsoor
parents:
diff changeset
114 # z_average[k[1],i]=1
e13ec2c3fabe planemo upload commit 25fd6a739741295e3f434e0be0286dee61e06825
mmonsoor
parents:
diff changeset
115 # }
e13ec2c3fabe planemo upload commit 25fd6a739741295e3f434e0be0286dee61e06825
mmonsoor
parents:
diff changeset
116 }
e13ec2c3fabe planemo upload commit 25fd6a739741295e3f434e0be0286dee61e06825
mmonsoor
parents:
diff changeset
117 else {
e13ec2c3fabe planemo upload commit 25fd6a739741295e3f434e0be0286dee61e06825
mmonsoor
parents:
diff changeset
118 prob_mean_ma <- probM
e13ec2c3fabe planemo upload commit 25fd6a739741295e3f434e0be0286dee61e06825
mmonsoor
parents:
diff changeset
119 }
e13ec2c3fabe planemo upload commit 25fd6a739741295e3f434e0be0286dee61e06825
mmonsoor
parents:
diff changeset
120 # think about natural probabilities
e13ec2c3fabe planemo upload commit 25fd6a739741295e3f434e0be0286dee61e06825
mmonsoor
parents:
diff changeset
121 # prob_mean_ma[prob_mean_ma[,1]!=0,1]/sum(prob_mean_ma[prob_mean_ma[,1]!=0,1])
e13ec2c3fabe planemo upload commit 25fd6a739741295e3f434e0be0286dee61e06825
mmonsoor
parents:
diff changeset
122 prob_mean_ma <- apply(prob_mean_ma, 2, function(x){ x[x!=0] <- x[x!=0]/sum(x[x!=0]); return(x)} )
e13ec2c3fabe planemo upload commit 25fd6a739741295e3f434e0be0286dee61e06825
mmonsoor
parents:
diff changeset
123
e13ec2c3fabe planemo upload commit 25fd6a739741295e3f434e0be0286dee61e06825
mmonsoor
parents:
diff changeset
124 # create a dir to figures
e13ec2c3fabe planemo upload commit 25fd6a739741295e3f434e0be0286dee61e06825
mmonsoor
parents:
diff changeset
125 lpattern <- function(type){
e13ec2c3fabe planemo upload commit 25fd6a739741295e3f434e0be0286dee61e06825
mmonsoor
parents:
diff changeset
126 switch(type,
e13ec2c3fabe planemo upload commit 25fd6a739741295e3f434e0be0286dee61e06825
mmonsoor
parents:
diff changeset
127 kegg = "http://www.genome.jp/dbget-bin/www_bget?",
e13ec2c3fabe planemo upload commit 25fd6a739741295e3f434e0be0286dee61e06825
mmonsoor
parents:
diff changeset
128 chebi = "http://www.ebi.ac.uk/chebi/searchId.do;EFB7DFF9E88306BBCD6AB78B32664A85?chebiId=",
e13ec2c3fabe planemo upload commit 25fd6a739741295e3f434e0be0286dee61e06825
mmonsoor
parents:
diff changeset
129 pubchem = "http://www.ncbi.nlm.nih.gov/pccompound/?term="
e13ec2c3fabe planemo upload commit 25fd6a739741295e3f434e0be0286dee61e06825
mmonsoor
parents:
diff changeset
130 )
e13ec2c3fabe planemo upload commit 25fd6a739741295e3f434e0be0286dee61e06825
mmonsoor
parents:
diff changeset
131 }
e13ec2c3fabe planemo upload commit 25fd6a739741295e3f434e0be0286dee61e06825
mmonsoor
parents:
diff changeset
132 linkURL <- lpattern(linkPattern)
e13ec2c3fabe planemo upload commit 25fd6a739741295e3f434e0be0286dee61e06825
mmonsoor
parents:
diff changeset
133 fig <- paste("file://", getwd(), paste("/",filename,"_fig/",sep=""), sep="")
e13ec2c3fabe planemo upload commit 25fd6a739741295e3f434e0be0286dee61e06825
mmonsoor
parents:
diff changeset
134 if(!is.null(molIon$cameraobj)) {
e13ec2c3fabe planemo upload commit 25fd6a739741295e3f434e0be0286dee61e06825
mmonsoor
parents:
diff changeset
135 figidx <- c("")
e13ec2c3fabe planemo upload commit 25fd6a739741295e3f434e0be0286dee61e06825
mmonsoor
parents:
diff changeset
136 coords <- gsub("(^\\d)","X\\1",rownames(molIon$cameraobj@xcmsSet@phenoData))
e13ec2c3fabe planemo upload commit 25fd6a739741295e3f434e0be0286dee61e06825
mmonsoor
parents:
diff changeset
137 # experimental! Which set of characters????
e13ec2c3fabe planemo upload commit 25fd6a739741295e3f434e0be0286dee61e06825
mmonsoor
parents:
diff changeset
138 coords <- gsub("-|\\,|~","\\.",coords)
e13ec2c3fabe planemo upload commit 25fd6a739741295e3f434e0be0286dee61e06825
mmonsoor
parents:
diff changeset
139 coords <- gsub("\\s+","\\.",coords)
e13ec2c3fabe planemo upload commit 25fd6a739741295e3f434e0be0286dee61e06825
mmonsoor
parents:
diff changeset
140 peaklist <- getPeaklist(molIon$cameraobj)
e13ec2c3fabe planemo upload commit 25fd6a739741295e3f434e0be0286dee61e06825
mmonsoor
parents:
diff changeset
141 rpeaklist <- peaklist[,c("mz","rt","isotopes","adduct","pcgroup")]
e13ec2c3fabe planemo upload commit 25fd6a739741295e3f434e0be0286dee61e06825
mmonsoor
parents:
diff changeset
142 }
e13ec2c3fabe planemo upload commit 25fd6a739741295e3f434e0be0286dee61e06825
mmonsoor
parents:
diff changeset
143 else {
e13ec2c3fabe planemo upload commit 25fd6a739741295e3f434e0be0286dee61e06825
mmonsoor
parents:
diff changeset
144 figidx <- c("","")
e13ec2c3fabe planemo upload commit 25fd6a739741295e3f434e0be0286dee61e06825
mmonsoor
parents:
diff changeset
145 coordsP <- gsub("(^\\d)","X\\1",rownames(molIon$pos@xcmsSet@phenoData))
e13ec2c3fabe planemo upload commit 25fd6a739741295e3f434e0be0286dee61e06825
mmonsoor
parents:
diff changeset
146 # experimental! Which set of characters????
e13ec2c3fabe planemo upload commit 25fd6a739741295e3f434e0be0286dee61e06825
mmonsoor
parents:
diff changeset
147 coordsP <- gsub("-|\\,|~","\\.",coordsP)
e13ec2c3fabe planemo upload commit 25fd6a739741295e3f434e0be0286dee61e06825
mmonsoor
parents:
diff changeset
148 coordsP <- gsub("\\s+","\\.",coordsP)
e13ec2c3fabe planemo upload commit 25fd6a739741295e3f434e0be0286dee61e06825
mmonsoor
parents:
diff changeset
149 coordsN <- gsub("(^\\d)","X\\1",rownames(molIon$neg@xcmsSet@phenoData))
e13ec2c3fabe planemo upload commit 25fd6a739741295e3f434e0be0286dee61e06825
mmonsoor
parents:
diff changeset
150 # experimental! Which set of characters????
e13ec2c3fabe planemo upload commit 25fd6a739741295e3f434e0be0286dee61e06825
mmonsoor
parents:
diff changeset
151 coordsN <- gsub("-|\\,|~","\\.",coordsN)
e13ec2c3fabe planemo upload commit 25fd6a739741295e3f434e0be0286dee61e06825
mmonsoor
parents:
diff changeset
152 coordsN <- gsub("\\s+","\\.",coordsN)
e13ec2c3fabe planemo upload commit 25fd6a739741295e3f434e0be0286dee61e06825
mmonsoor
parents:
diff changeset
153 coords <- coordsP
e13ec2c3fabe planemo upload commit 25fd6a739741295e3f434e0be0286dee61e06825
mmonsoor
parents:
diff changeset
154 if(length(coordsP)!=length(coordsN)) cat("\n Warning: The number of samples are different\n")
e13ec2c3fabe planemo upload commit 25fd6a739741295e3f434e0be0286dee61e06825
mmonsoor
parents:
diff changeset
155
e13ec2c3fabe planemo upload commit 25fd6a739741295e3f434e0be0286dee61e06825
mmonsoor
parents:
diff changeset
156 peaklistP <- getPeaklist(molIon$pos)
e13ec2c3fabe planemo upload commit 25fd6a739741295e3f434e0be0286dee61e06825
mmonsoor
parents:
diff changeset
157 rpeaklistP <- peaklistP[,c("mz","rt","isotopes","adduct","pcgroup")]
e13ec2c3fabe planemo upload commit 25fd6a739741295e3f434e0be0286dee61e06825
mmonsoor
parents:
diff changeset
158 peaklistN <- getPeaklist(molIon$neg)
e13ec2c3fabe planemo upload commit 25fd6a739741295e3f434e0be0286dee61e06825
mmonsoor
parents:
diff changeset
159 rpeaklistN <- peaklistN[,c("mz","rt","isotopes","adduct","pcgroup")]
e13ec2c3fabe planemo upload commit 25fd6a739741295e3f434e0be0286dee61e06825
mmonsoor
parents:
diff changeset
160 }
e13ec2c3fabe planemo upload commit 25fd6a739741295e3f434e0be0286dee61e06825
mmonsoor
parents:
diff changeset
161
e13ec2c3fabe planemo upload commit 25fd6a739741295e3f434e0be0286dee61e06825
mmonsoor
parents:
diff changeset
162 # if(sum(is.na(peaklist))) {
e13ec2c3fabe planemo upload commit 25fd6a739741295e3f434e0be0286dee61e06825
mmonsoor
parents:
diff changeset
163 # cat("\nWarning: NAs Found in peaklist\n\nSubstituting for \"ones\"\n")
e13ec2c3fabe planemo upload commit 25fd6a739741295e3f434e0be0286dee61e06825
mmonsoor
parents:
diff changeset
164 # na.ids <- which(is.na(peaklist),arr.ind=TRUE)
e13ec2c3fabe planemo upload commit 25fd6a739741295e3f434e0be0286dee61e06825
mmonsoor
parents:
diff changeset
165 # for(l in 1:nrow(na.ids)){
e13ec2c3fabe planemo upload commit 25fd6a739741295e3f434e0be0286dee61e06825
mmonsoor
parents:
diff changeset
166 # peaklist[na.ids[l,][1], na.ids[l,][2]] <- 1
e13ec2c3fabe planemo upload commit 25fd6a739741295e3f434e0be0286dee61e06825
mmonsoor
parents:
diff changeset
167 # }
e13ec2c3fabe planemo upload commit 25fd6a739741295e3f434e0be0286dee61e06825
mmonsoor
parents:
diff changeset
168 # }
e13ec2c3fabe planemo upload commit 25fd6a739741295e3f434e0be0286dee61e06825
mmonsoor
parents:
diff changeset
169 #
e13ec2c3fabe planemo upload commit 25fd6a739741295e3f434e0be0286dee61e06825
mmonsoor
parents:
diff changeset
170
e13ec2c3fabe planemo upload commit 25fd6a739741295e3f434e0be0286dee61e06825
mmonsoor
parents:
diff changeset
171 ans <- matrix("", nrow=1, ncol=7+length(coords))
e13ec2c3fabe planemo upload commit 25fd6a739741295e3f434e0be0286dee61e06825
mmonsoor
parents:
diff changeset
172 unq <- unique(ReactMatrix[,2:3])
e13ec2c3fabe planemo upload commit 25fd6a739741295e3f434e0be0286dee61e06825
mmonsoor
parents:
diff changeset
173 for (i in 1:nrow(unq)) {
e13ec2c3fabe planemo upload commit 25fd6a739741295e3f434e0be0286dee61e06825
mmonsoor
parents:
diff changeset
174 coord <- which(ReactMatrix[,2]==unq[i,1] & ReactMatrix[,3]==unq[i,2])
e13ec2c3fabe planemo upload commit 25fd6a739741295e3f434e0be0286dee61e06825
mmonsoor
parents:
diff changeset
175 coord2 <- which(allion[,2]==unq[i,1] & allion[,1]==unq[i,2])
e13ec2c3fabe planemo upload commit 25fd6a739741295e3f434e0be0286dee61e06825
mmonsoor
parents:
diff changeset
176 # idx2 <- unique(which(allion[,1] %in% reactionM[reactionM[,5]=="unknown",2]))
e13ec2c3fabe planemo upload commit 25fd6a739741295e3f434e0be0286dee61e06825
mmonsoor
parents:
diff changeset
177 # work with the higher intensities for a given ion annotation, not necessarily the right one
e13ec2c3fabe planemo upload commit 25fd6a739741295e3f434e0be0286dee61e06825
mmonsoor
parents:
diff changeset
178
e13ec2c3fabe planemo upload commit 25fd6a739741295e3f434e0be0286dee61e06825
mmonsoor
parents:
diff changeset
179 if(!is.null(molIon$cameraobj)) {
e13ec2c3fabe planemo upload commit 25fd6a739741295e3f434e0be0286dee61e06825
mmonsoor
parents:
diff changeset
180 idx <- as.vector(unlist(sapply(allion[coord2,"trace"],
e13ec2c3fabe planemo upload commit 25fd6a739741295e3f434e0be0286dee61e06825
mmonsoor
parents:
diff changeset
181 function(x) {
e13ec2c3fabe planemo upload commit 25fd6a739741295e3f434e0be0286dee61e06825
mmonsoor
parents:
diff changeset
182 x <- as.matrix(x)
e13ec2c3fabe planemo upload commit 25fd6a739741295e3f434e0be0286dee61e06825
mmonsoor
parents:
diff changeset
183 raw <- strsplit(x,";")[[1]]
e13ec2c3fabe planemo upload commit 25fd6a739741295e3f434e0be0286dee61e06825
mmonsoor
parents:
diff changeset
184 mraw <- apply(peaklist[raw, coords], 1, mean)
e13ec2c3fabe planemo upload commit 25fd6a739741295e3f434e0be0286dee61e06825
mmonsoor
parents:
diff changeset
185 raw[which.max(mraw)]
e13ec2c3fabe planemo upload commit 25fd6a739741295e3f434e0be0286dee61e06825
mmonsoor
parents:
diff changeset
186 }
e13ec2c3fabe planemo upload commit 25fd6a739741295e3f434e0be0286dee61e06825
mmonsoor
parents:
diff changeset
187
e13ec2c3fabe planemo upload commit 25fd6a739741295e3f434e0be0286dee61e06825
mmonsoor
parents:
diff changeset
188 )
e13ec2c3fabe planemo upload commit 25fd6a739741295e3f434e0be0286dee61e06825
mmonsoor
parents:
diff changeset
189 )
e13ec2c3fabe planemo upload commit 25fd6a739741295e3f434e0be0286dee61e06825
mmonsoor
parents:
diff changeset
190 )
e13ec2c3fabe planemo upload commit 25fd6a739741295e3f434e0be0286dee61e06825
mmonsoor
parents:
diff changeset
191
e13ec2c3fabe planemo upload commit 25fd6a739741295e3f434e0be0286dee61e06825
mmonsoor
parents:
diff changeset
192 idx <- unique(idx)
e13ec2c3fabe planemo upload commit 25fd6a739741295e3f434e0be0286dee61e06825
mmonsoor
parents:
diff changeset
193 figidx <- append(figidx,idx)
e13ec2c3fabe planemo upload commit 25fd6a739741295e3f434e0be0286dee61e06825
mmonsoor
parents:
diff changeset
194 }
e13ec2c3fabe planemo upload commit 25fd6a739741295e3f434e0be0286dee61e06825
mmonsoor
parents:
diff changeset
195 else {
e13ec2c3fabe planemo upload commit 25fd6a739741295e3f434e0be0286dee61e06825
mmonsoor
parents:
diff changeset
196 idx <- c()
e13ec2c3fabe planemo upload commit 25fd6a739741295e3f434e0be0286dee61e06825
mmonsoor
parents:
diff changeset
197
e13ec2c3fabe planemo upload commit 25fd6a739741295e3f434e0be0286dee61e06825
mmonsoor
parents:
diff changeset
198 for(l in 1:nrow( allion[coord2,c("trace","comb")])) {
e13ec2c3fabe planemo upload commit 25fd6a739741295e3f434e0be0286dee61e06825
mmonsoor
parents:
diff changeset
199 x <- as.matrix(allion[coord2,c("trace","comb")][l,])
e13ec2c3fabe planemo upload commit 25fd6a739741295e3f434e0be0286dee61e06825
mmonsoor
parents:
diff changeset
200 raw <- strsplit(x[1],";")[[1]]
e13ec2c3fabe planemo upload commit 25fd6a739741295e3f434e0be0286dee61e06825
mmonsoor
parents:
diff changeset
201 if(x[2]!="neg"){
e13ec2c3fabe planemo upload commit 25fd6a739741295e3f434e0be0286dee61e06825
mmonsoor
parents:
diff changeset
202 mraw <- apply(peaklistP[raw, coordsP], 1, mean, na.rm=TRUE)
e13ec2c3fabe planemo upload commit 25fd6a739741295e3f434e0be0286dee61e06825
mmonsoor
parents:
diff changeset
203 }
e13ec2c3fabe planemo upload commit 25fd6a739741295e3f434e0be0286dee61e06825
mmonsoor
parents:
diff changeset
204 else {
e13ec2c3fabe planemo upload commit 25fd6a739741295e3f434e0be0286dee61e06825
mmonsoor
parents:
diff changeset
205
e13ec2c3fabe planemo upload commit 25fd6a739741295e3f434e0be0286dee61e06825
mmonsoor
parents:
diff changeset
206 mraw <- apply(peaklistN[raw, coordsN], 1, mean, na.rm=TRUE)
e13ec2c3fabe planemo upload commit 25fd6a739741295e3f434e0be0286dee61e06825
mmonsoor
parents:
diff changeset
207 }
e13ec2c3fabe planemo upload commit 25fd6a739741295e3f434e0be0286dee61e06825
mmonsoor
parents:
diff changeset
208 idx <- c(idx, raw[which.max(mraw)])
e13ec2c3fabe planemo upload commit 25fd6a739741295e3f434e0be0286dee61e06825
mmonsoor
parents:
diff changeset
209 }
e13ec2c3fabe planemo upload commit 25fd6a739741295e3f434e0be0286dee61e06825
mmonsoor
parents:
diff changeset
210
e13ec2c3fabe planemo upload commit 25fd6a739741295e3f434e0be0286dee61e06825
mmonsoor
parents:
diff changeset
211
e13ec2c3fabe planemo upload commit 25fd6a739741295e3f434e0be0286dee61e06825
mmonsoor
parents:
diff changeset
212 idx <- unique(idx)
e13ec2c3fabe planemo upload commit 25fd6a739741295e3f434e0be0286dee61e06825
mmonsoor
parents:
diff changeset
213 figidx <- rbind(figidx,c(idx,allion[coord2,"comb"][1]))
e13ec2c3fabe planemo upload commit 25fd6a739741295e3f434e0be0286dee61e06825
mmonsoor
parents:
diff changeset
214 }
e13ec2c3fabe planemo upload commit 25fd6a739741295e3f434e0be0286dee61e06825
mmonsoor
parents:
diff changeset
215 #figidx <- append(figidx,strsplit(allion[coord2,5], ";")[[1]][1])
e13ec2c3fabe planemo upload commit 25fd6a739741295e3f434e0be0286dee61e06825
mmonsoor
parents:
diff changeset
216 ans1 <- matrix("", nrow=length(coord), ncol=7+length(coords))
e13ec2c3fabe planemo upload commit 25fd6a739741295e3f434e0be0286dee61e06825
mmonsoor
parents:
diff changeset
217 ans1[,2]<-as.matrix(ReactMatrix[coord,5])
e13ec2c3fabe planemo upload commit 25fd6a739741295e3f434e0be0286dee61e06825
mmonsoor
parents:
diff changeset
218 prob <- as.matrix(prob_mean_ma[coord, i]) # need to change and compare a pair of mass/rt
e13ec2c3fabe planemo upload commit 25fd6a739741295e3f434e0be0286dee61e06825
mmonsoor
parents:
diff changeset
219 # number figs
e13ec2c3fabe planemo upload commit 25fd6a739741295e3f434e0be0286dee61e06825
mmonsoor
parents:
diff changeset
220 if ( i >= 100 ) { ans1[1,6]=i }
e13ec2c3fabe planemo upload commit 25fd6a739741295e3f434e0be0286dee61e06825
mmonsoor
parents:
diff changeset
221 else { if ( i >= 10 ) { ans1[1,6]=paste(0,i, sep="") } else { ans1[1,6]=paste("00",i, sep="") } }
e13ec2c3fabe planemo upload commit 25fd6a739741295e3f434e0be0286dee61e06825
mmonsoor
parents:
diff changeset
222
e13ec2c3fabe planemo upload commit 25fd6a739741295e3f434e0be0286dee61e06825
mmonsoor
parents:
diff changeset
223 if (sum(prob)>0) {
e13ec2c3fabe planemo upload commit 25fd6a739741295e3f434e0be0286dee61e06825
mmonsoor
parents:
diff changeset
224 #prob <- prob/sum(prob)
e13ec2c3fabe planemo upload commit 25fd6a739741295e3f434e0be0286dee61e06825
mmonsoor
parents:
diff changeset
225 o <- order(prob, decreasing=TRUE)
e13ec2c3fabe planemo upload commit 25fd6a739741295e3f434e0be0286dee61e06825
mmonsoor
parents:
diff changeset
226 ans1[,-6] <- ans1[o,-6]
e13ec2c3fabe planemo upload commit 25fd6a739741295e3f434e0be0286dee61e06825
mmonsoor
parents:
diff changeset
227 ans1 <- matrix(ans1, nrow=length(o))
e13ec2c3fabe planemo upload commit 25fd6a739741295e3f434e0be0286dee61e06825
mmonsoor
parents:
diff changeset
228 ans1[1,1] <- ReactMatrix[coord[1],3]
e13ec2c3fabe planemo upload commit 25fd6a739741295e3f434e0be0286dee61e06825
mmonsoor
parents:
diff changeset
229 #ans1[,3] <- round(prob/min(prob[prob!=0]), 3)[o]
e13ec2c3fabe planemo upload commit 25fd6a739741295e3f434e0be0286dee61e06825
mmonsoor
parents:
diff changeset
230 ans1[,3] <- round(prob, 3)[o]
e13ec2c3fabe planemo upload commit 25fd6a739741295e3f434e0be0286dee61e06825
mmonsoor
parents:
diff changeset
231 if (length(prob[prob!=0])>1) {
e13ec2c3fabe planemo upload commit 25fd6a739741295e3f434e0be0286dee61e06825
mmonsoor
parents:
diff changeset
232 entropy <- -sum(prob[prob!=0]*log(prob[prob!=0], length(prob[prob!=0])))
e13ec2c3fabe planemo upload commit 25fd6a739741295e3f434e0be0286dee61e06825
mmonsoor
parents:
diff changeset
233 }
e13ec2c3fabe planemo upload commit 25fd6a739741295e3f434e0be0286dee61e06825
mmonsoor
parents:
diff changeset
234 else { entropy <- 0
e13ec2c3fabe planemo upload commit 25fd6a739741295e3f434e0be0286dee61e06825
mmonsoor
parents:
diff changeset
235 }
e13ec2c3fabe planemo upload commit 25fd6a739741295e3f434e0be0286dee61e06825
mmonsoor
parents:
diff changeset
236 ans1[1,4] <- round(entropy, 3)
e13ec2c3fabe planemo upload commit 25fd6a739741295e3f434e0be0286dee61e06825
mmonsoor
parents:
diff changeset
237 }
e13ec2c3fabe planemo upload commit 25fd6a739741295e3f434e0be0286dee61e06825
mmonsoor
parents:
diff changeset
238 else {
e13ec2c3fabe planemo upload commit 25fd6a739741295e3f434e0be0286dee61e06825
mmonsoor
parents:
diff changeset
239 ans1[1,1] <- ReactMatrix[coord[1],3]
e13ec2c3fabe planemo upload commit 25fd6a739741295e3f434e0be0286dee61e06825
mmonsoor
parents:
diff changeset
240 ans1[1,3] <- "undef"
e13ec2c3fabe planemo upload commit 25fd6a739741295e3f434e0be0286dee61e06825
mmonsoor
parents:
diff changeset
241 }
e13ec2c3fabe planemo upload commit 25fd6a739741295e3f434e0be0286dee61e06825
mmonsoor
parents:
diff changeset
242
e13ec2c3fabe planemo upload commit 25fd6a739741295e3f434e0be0286dee61e06825
mmonsoor
parents:
diff changeset
243 if(!is.null(molIon$cameraobj)) {
e13ec2c3fabe planemo upload commit 25fd6a739741295e3f434e0be0286dee61e06825
mmonsoor
parents:
diff changeset
244 ans1[1,7] <- apply(rpeaklist[idx,], 1, function(x) paste(x[c(1,2,3,4)], collapse="#"))
e13ec2c3fabe planemo upload commit 25fd6a739741295e3f434e0be0286dee61e06825
mmonsoor
parents:
diff changeset
245 ans1[1,8:ncol(ans1)] <- as.matrix(peaklist[idx, coords])
e13ec2c3fabe planemo upload commit 25fd6a739741295e3f434e0be0286dee61e06825
mmonsoor
parents:
diff changeset
246 }
e13ec2c3fabe planemo upload commit 25fd6a739741295e3f434e0be0286dee61e06825
mmonsoor
parents:
diff changeset
247 else {
e13ec2c3fabe planemo upload commit 25fd6a739741295e3f434e0be0286dee61e06825
mmonsoor
parents:
diff changeset
248 if(allion[coord2,"comb"]=="pos"|allion[coord2,"comb"]=="both") {
e13ec2c3fabe planemo upload commit 25fd6a739741295e3f434e0be0286dee61e06825
mmonsoor
parents:
diff changeset
249 ans1[1,7] <- apply(rpeaklistP[idx,], 1, function(x) paste(x[c(1,2,3,4)], collapse="#"))
e13ec2c3fabe planemo upload commit 25fd6a739741295e3f434e0be0286dee61e06825
mmonsoor
parents:
diff changeset
250 ans1[1,8:ncol(ans1)] <- as.matrix(peaklistP[idx, coordsP])
e13ec2c3fabe planemo upload commit 25fd6a739741295e3f434e0be0286dee61e06825
mmonsoor
parents:
diff changeset
251 }
e13ec2c3fabe planemo upload commit 25fd6a739741295e3f434e0be0286dee61e06825
mmonsoor
parents:
diff changeset
252 else {
e13ec2c3fabe planemo upload commit 25fd6a739741295e3f434e0be0286dee61e06825
mmonsoor
parents:
diff changeset
253 ans1[1,7] <- apply(rpeaklistN[idx,], 1, function(x) paste(x[c(1,2,3,4)], collapse="#"))
e13ec2c3fabe planemo upload commit 25fd6a739741295e3f434e0be0286dee61e06825
mmonsoor
parents:
diff changeset
254 ans1[1,8:ncol(ans1)] <- as.matrix(peaklistN[idx, coordsN])
e13ec2c3fabe planemo upload commit 25fd6a739741295e3f434e0be0286dee61e06825
mmonsoor
parents:
diff changeset
255 }
e13ec2c3fabe planemo upload commit 25fd6a739741295e3f434e0be0286dee61e06825
mmonsoor
parents:
diff changeset
256 }
e13ec2c3fabe planemo upload commit 25fd6a739741295e3f434e0be0286dee61e06825
mmonsoor
parents:
diff changeset
257 ans <- rbind(ans, as.matrix(ans1))
e13ec2c3fabe planemo upload commit 25fd6a739741295e3f434e0be0286dee61e06825
mmonsoor
parents:
diff changeset
258 }
e13ec2c3fabe planemo upload commit 25fd6a739741295e3f434e0be0286dee61e06825
mmonsoor
parents:
diff changeset
259 ans <- ans[-1,]
e13ec2c3fabe planemo upload commit 25fd6a739741295e3f434e0be0286dee61e06825
mmonsoor
parents:
diff changeset
260 # this option should change according with the bank
e13ec2c3fabe planemo upload commit 25fd6a739741295e3f434e0be0286dee61e06825
mmonsoor
parents:
diff changeset
261 if(html) {
e13ec2c3fabe planemo upload commit 25fd6a739741295e3f434e0be0286dee61e06825
mmonsoor
parents:
diff changeset
262 nid <- unlist(sapply(ans[,2], function(x) which(DB$id==x)))
e13ec2c3fabe planemo upload commit 25fd6a739741295e3f434e0be0286dee61e06825
mmonsoor
parents:
diff changeset
263 #ans[,2] <- as.character(DB$name[nid])
e13ec2c3fabe planemo upload commit 25fd6a739741295e3f434e0be0286dee61e06825
mmonsoor
parents:
diff changeset
264 }
e13ec2c3fabe planemo upload commit 25fd6a739741295e3f434e0be0286dee61e06825
mmonsoor
parents:
diff changeset
265 unk <- reactionM[reactionM[,5]=="unknown",]
e13ec2c3fabe planemo upload commit 25fd6a739741295e3f434e0be0286dee61e06825
mmonsoor
parents:
diff changeset
266 ans1 <- matrix("", nrow=nrow(unk), ncol=7+length(coords))
e13ec2c3fabe planemo upload commit 25fd6a739741295e3f434e0be0286dee61e06825
mmonsoor
parents:
diff changeset
267 ans1[,1] <- unk[,3]
e13ec2c3fabe planemo upload commit 25fd6a739741295e3f434e0be0286dee61e06825
mmonsoor
parents:
diff changeset
268 ans1[,2] <- unk[,5]
e13ec2c3fabe planemo upload commit 25fd6a739741295e3f434e0be0286dee61e06825
mmonsoor
parents:
diff changeset
269 for(j in 1:nrow(ans1)) {
e13ec2c3fabe planemo upload commit 25fd6a739741295e3f434e0be0286dee61e06825
mmonsoor
parents:
diff changeset
270 i <- j + max(as.numeric(ans[,6]),na.rm=TRUE)
e13ec2c3fabe planemo upload commit 25fd6a739741295e3f434e0be0286dee61e06825
mmonsoor
parents:
diff changeset
271 if ( i >= 100 ) { ans1[j,6]=i }
e13ec2c3fabe planemo upload commit 25fd6a739741295e3f434e0be0286dee61e06825
mmonsoor
parents:
diff changeset
272 else { if ( i >= 10 ) { ans1[j,6]=paste(0,i, sep="") } else { ans1[j,6]=paste("00",i, sep="") } }
e13ec2c3fabe planemo upload commit 25fd6a739741295e3f434e0be0286dee61e06825
mmonsoor
parents:
diff changeset
273 }
e13ec2c3fabe planemo upload commit 25fd6a739741295e3f434e0be0286dee61e06825
mmonsoor
parents:
diff changeset
274 # this step try to recover ids of ion annotation for masses without database annotation
e13ec2c3fabe planemo upload commit 25fd6a739741295e3f434e0be0286dee61e06825
mmonsoor
parents:
diff changeset
275 idx2 <- c(); #for(m in 1:nrow(allion)) if(sum(allion[m,2]==as.numeric(unk[,2])) & sum(allion[m,1]==as.numeric(unk[,3]))) idx2 <- append(idx2, m)
e13ec2c3fabe planemo upload commit 25fd6a739741295e3f434e0be0286dee61e06825
mmonsoor
parents:
diff changeset
276 # temp changes made 03/03/2014 have to check carefuly
e13ec2c3fabe planemo upload commit 25fd6a739741295e3f434e0be0286dee61e06825
mmonsoor
parents:
diff changeset
277 lidx <- lapply(1:nrow(allion), function(m) which(allion[m,2]==unk[,2] & allion[m,1]==unk[,3]))
e13ec2c3fabe planemo upload commit 25fd6a739741295e3f434e0be0286dee61e06825
mmonsoor
parents:
diff changeset
278 idx2 <- which(lapply(lidx, length)>0)
e13ec2c3fabe planemo upload commit 25fd6a739741295e3f434e0be0286dee61e06825
mmonsoor
parents:
diff changeset
279
e13ec2c3fabe planemo upload commit 25fd6a739741295e3f434e0be0286dee61e06825
mmonsoor
parents:
diff changeset
280 if(!is.null(molIon$cameraobj)) {
e13ec2c3fabe planemo upload commit 25fd6a739741295e3f434e0be0286dee61e06825
mmonsoor
parents:
diff changeset
281 idx <- as.vector(unlist(sapply(allion[idx2,"trace"],
e13ec2c3fabe planemo upload commit 25fd6a739741295e3f434e0be0286dee61e06825
mmonsoor
parents:
diff changeset
282 function(x) {
e13ec2c3fabe planemo upload commit 25fd6a739741295e3f434e0be0286dee61e06825
mmonsoor
parents:
diff changeset
283 x <- as.matrix(x)
e13ec2c3fabe planemo upload commit 25fd6a739741295e3f434e0be0286dee61e06825
mmonsoor
parents:
diff changeset
284 raw <- strsplit(x,";")[[1]]
e13ec2c3fabe planemo upload commit 25fd6a739741295e3f434e0be0286dee61e06825
mmonsoor
parents:
diff changeset
285 mraw <- apply(peaklist[raw, coords], 1, mean)
e13ec2c3fabe planemo upload commit 25fd6a739741295e3f434e0be0286dee61e06825
mmonsoor
parents:
diff changeset
286 raw[which.max(mraw)]
e13ec2c3fabe planemo upload commit 25fd6a739741295e3f434e0be0286dee61e06825
mmonsoor
parents:
diff changeset
287 }
e13ec2c3fabe planemo upload commit 25fd6a739741295e3f434e0be0286dee61e06825
mmonsoor
parents:
diff changeset
288
e13ec2c3fabe planemo upload commit 25fd6a739741295e3f434e0be0286dee61e06825
mmonsoor
parents:
diff changeset
289 )
e13ec2c3fabe planemo upload commit 25fd6a739741295e3f434e0be0286dee61e06825
mmonsoor
parents:
diff changeset
290 )
e13ec2c3fabe planemo upload commit 25fd6a739741295e3f434e0be0286dee61e06825
mmonsoor
parents:
diff changeset
291 )
e13ec2c3fabe planemo upload commit 25fd6a739741295e3f434e0be0286dee61e06825
mmonsoor
parents:
diff changeset
292 }
e13ec2c3fabe planemo upload commit 25fd6a739741295e3f434e0be0286dee61e06825
mmonsoor
parents:
diff changeset
293
e13ec2c3fabe planemo upload commit 25fd6a739741295e3f434e0be0286dee61e06825
mmonsoor
parents:
diff changeset
294 else {
e13ec2c3fabe planemo upload commit 25fd6a739741295e3f434e0be0286dee61e06825
mmonsoor
parents:
diff changeset
295 # don't know what happened here with apply
e13ec2c3fabe planemo upload commit 25fd6a739741295e3f434e0be0286dee61e06825
mmonsoor
parents:
diff changeset
296 idx <- c()
e13ec2c3fabe planemo upload commit 25fd6a739741295e3f434e0be0286dee61e06825
mmonsoor
parents:
diff changeset
297
e13ec2c3fabe planemo upload commit 25fd6a739741295e3f434e0be0286dee61e06825
mmonsoor
parents:
diff changeset
298 for(i in 1:nrow( allion[idx2,c("trace","comb")])) {
e13ec2c3fabe planemo upload commit 25fd6a739741295e3f434e0be0286dee61e06825
mmonsoor
parents:
diff changeset
299 x <- as.matrix(allion[idx2,c("trace","comb")][i,])
e13ec2c3fabe planemo upload commit 25fd6a739741295e3f434e0be0286dee61e06825
mmonsoor
parents:
diff changeset
300 raw <- strsplit(x[1],";")[[1]]
e13ec2c3fabe planemo upload commit 25fd6a739741295e3f434e0be0286dee61e06825
mmonsoor
parents:
diff changeset
301 if(x[2]!="neg"){
e13ec2c3fabe planemo upload commit 25fd6a739741295e3f434e0be0286dee61e06825
mmonsoor
parents:
diff changeset
302 mraw <- apply(peaklistP[raw, coordsP], 1, mean, na.rm=TRUE)
e13ec2c3fabe planemo upload commit 25fd6a739741295e3f434e0be0286dee61e06825
mmonsoor
parents:
diff changeset
303 }
e13ec2c3fabe planemo upload commit 25fd6a739741295e3f434e0be0286dee61e06825
mmonsoor
parents:
diff changeset
304 else {
e13ec2c3fabe planemo upload commit 25fd6a739741295e3f434e0be0286dee61e06825
mmonsoor
parents:
diff changeset
305
e13ec2c3fabe planemo upload commit 25fd6a739741295e3f434e0be0286dee61e06825
mmonsoor
parents:
diff changeset
306 mraw <- apply(peaklistN[raw, coordsN], 1, mean, na.rm=TRUE)
e13ec2c3fabe planemo upload commit 25fd6a739741295e3f434e0be0286dee61e06825
mmonsoor
parents:
diff changeset
307 }
e13ec2c3fabe planemo upload commit 25fd6a739741295e3f434e0be0286dee61e06825
mmonsoor
parents:
diff changeset
308 idx <- c(idx, raw[which.max(mraw)])
e13ec2c3fabe planemo upload commit 25fd6a739741295e3f434e0be0286dee61e06825
mmonsoor
parents:
diff changeset
309 }
e13ec2c3fabe planemo upload commit 25fd6a739741295e3f434e0be0286dee61e06825
mmonsoor
parents:
diff changeset
310
e13ec2c3fabe planemo upload commit 25fd6a739741295e3f434e0be0286dee61e06825
mmonsoor
parents:
diff changeset
311
e13ec2c3fabe planemo upload commit 25fd6a739741295e3f434e0be0286dee61e06825
mmonsoor
parents:
diff changeset
312
e13ec2c3fabe planemo upload commit 25fd6a739741295e3f434e0be0286dee61e06825
mmonsoor
parents:
diff changeset
313 tmpidx <- cbind(idx,allion[idx2,"comb"])
e13ec2c3fabe planemo upload commit 25fd6a739741295e3f434e0be0286dee61e06825
mmonsoor
parents:
diff changeset
314 }
e13ec2c3fabe planemo upload commit 25fd6a739741295e3f434e0be0286dee61e06825
mmonsoor
parents:
diff changeset
315 if(!is.null(molIon$cameraobj)) {
e13ec2c3fabe planemo upload commit 25fd6a739741295e3f434e0be0286dee61e06825
mmonsoor
parents:
diff changeset
316 ans1[,7] <- apply(rpeaklist[idx,], 1, function(x) paste(x[c(1,2,3,4)], collapse="#"))
e13ec2c3fabe planemo upload commit 25fd6a739741295e3f434e0be0286dee61e06825
mmonsoor
parents:
diff changeset
317 ans1[,8:ncol(ans1)] <- as.matrix(peaklist[idx, coords])
e13ec2c3fabe planemo upload commit 25fd6a739741295e3f434e0be0286dee61e06825
mmonsoor
parents:
diff changeset
318 }
e13ec2c3fabe planemo upload commit 25fd6a739741295e3f434e0be0286dee61e06825
mmonsoor
parents:
diff changeset
319 else {
e13ec2c3fabe planemo upload commit 25fd6a739741295e3f434e0be0286dee61e06825
mmonsoor
parents:
diff changeset
320 idxP <- tmpidx[tmpidx[,2]!="neg",1]
e13ec2c3fabe planemo upload commit 25fd6a739741295e3f434e0be0286dee61e06825
mmonsoor
parents:
diff changeset
321 ans1[1:length(idxP),7] <- apply(rpeaklistP[idxP,], 1, function(x) paste(x[c(1,2,3,4)], collapse="#"))
e13ec2c3fabe planemo upload commit 25fd6a739741295e3f434e0be0286dee61e06825
mmonsoor
parents:
diff changeset
322 ans1[1:length(idxP),8:ncol(ans1)] <- as.matrix(peaklistP[idxP, coordsP])
e13ec2c3fabe planemo upload commit 25fd6a739741295e3f434e0be0286dee61e06825
mmonsoor
parents:
diff changeset
323 idxN <- tmpidx[tmpidx[,2]=="neg",1]
e13ec2c3fabe planemo upload commit 25fd6a739741295e3f434e0be0286dee61e06825
mmonsoor
parents:
diff changeset
324 ans1[(length(idxP)+1):nrow(ans1),7] <- apply(rpeaklistN[idxN,], 1, function(x) paste(x[c(1,2,3,4)], collapse="#"))
e13ec2c3fabe planemo upload commit 25fd6a739741295e3f434e0be0286dee61e06825
mmonsoor
parents:
diff changeset
325 ans1[(length(idxP)+1):nrow(ans1),8:ncol(ans1)] <- as.matrix(peaklistN[idxN, coordsN])
e13ec2c3fabe planemo upload commit 25fd6a739741295e3f434e0be0286dee61e06825
mmonsoor
parents:
diff changeset
326 }
e13ec2c3fabe planemo upload commit 25fd6a739741295e3f434e0be0286dee61e06825
mmonsoor
parents:
diff changeset
327 ans <- rbind(ans, as.matrix(ans1))
e13ec2c3fabe planemo upload commit 25fd6a739741295e3f434e0be0286dee61e06825
mmonsoor
parents:
diff changeset
328
e13ec2c3fabe planemo upload commit 25fd6a739741295e3f434e0be0286dee61e06825
mmonsoor
parents:
diff changeset
329 if(!is.null(molIon$cameraobj)) {
e13ec2c3fabe planemo upload commit 25fd6a739741295e3f434e0be0286dee61e06825
mmonsoor
parents:
diff changeset
330 figidx <- c(figidx,idx)
e13ec2c3fabe planemo upload commit 25fd6a739741295e3f434e0be0286dee61e06825
mmonsoor
parents:
diff changeset
331 figidx <- as.numeric(figidx[-1])
e13ec2c3fabe planemo upload commit 25fd6a739741295e3f434e0be0286dee61e06825
mmonsoor
parents:
diff changeset
332 }
e13ec2c3fabe planemo upload commit 25fd6a739741295e3f434e0be0286dee61e06825
mmonsoor
parents:
diff changeset
333 else {
e13ec2c3fabe planemo upload commit 25fd6a739741295e3f434e0be0286dee61e06825
mmonsoor
parents:
diff changeset
334 figidx <- rbind(figidx,tmpidx)
e13ec2c3fabe planemo upload commit 25fd6a739741295e3f434e0be0286dee61e06825
mmonsoor
parents:
diff changeset
335 allidx <- figidx[-1,]
e13ec2c3fabe planemo upload commit 25fd6a739741295e3f434e0be0286dee61e06825
mmonsoor
parents:
diff changeset
336 allidx <- cbind(allidx, ans[ans[,6]!="",6])
e13ec2c3fabe planemo upload commit 25fd6a739741295e3f434e0be0286dee61e06825
mmonsoor
parents:
diff changeset
337 figidx <- as.numeric(figidx[-1,1])
e13ec2c3fabe planemo upload commit 25fd6a739741295e3f434e0be0286dee61e06825
mmonsoor
parents:
diff changeset
338 }
e13ec2c3fabe planemo upload commit 25fd6a739741295e3f434e0be0286dee61e06825
mmonsoor
parents:
diff changeset
339
e13ec2c3fabe planemo upload commit 25fd6a739741295e3f434e0be0286dee61e06825
mmonsoor
parents:
diff changeset
340
e13ec2c3fabe planemo upload commit 25fd6a739741295e3f434e0be0286dee61e06825
mmonsoor
parents:
diff changeset
341 if(m.test=="none") {
e13ec2c3fabe planemo upload commit 25fd6a739741295e3f434e0be0286dee61e06825
mmonsoor
parents:
diff changeset
342 testname <- "none"
e13ec2c3fabe planemo upload commit 25fd6a739741295e3f434e0be0286dee61e06825
mmonsoor
parents:
diff changeset
343 #testname <- "Formula"
e13ec2c3fabe planemo upload commit 25fd6a739741295e3f434e0be0286dee61e06825
mmonsoor
parents:
diff changeset
344 #ans[ans[,2]!="unknown",][,5] <- as.character(DB$formula[nid])
e13ec2c3fabe planemo upload commit 25fd6a739741295e3f434e0be0286dee61e06825
mmonsoor
parents:
diff changeset
345 }
e13ec2c3fabe planemo upload commit 25fd6a739741295e3f434e0be0286dee61e06825
mmonsoor
parents:
diff changeset
346 if(m.test=="t.test") {
e13ec2c3fabe planemo upload commit 25fd6a739741295e3f434e0be0286dee61e06825
mmonsoor
parents:
diff changeset
347 normalize.medFC <- function(mat) {
e13ec2c3fabe planemo upload commit 25fd6a739741295e3f434e0be0286dee61e06825
mmonsoor
parents:
diff changeset
348 # Perform median fold change normalisation
e13ec2c3fabe planemo upload commit 25fd6a739741295e3f434e0be0286dee61e06825
mmonsoor
parents:
diff changeset
349 # X - data set [Variables & Samples]
e13ec2c3fabe planemo upload commit 25fd6a739741295e3f434e0be0286dee61e06825
mmonsoor
parents:
diff changeset
350 medSam <- apply(mat, 1, median)
e13ec2c3fabe planemo upload commit 25fd6a739741295e3f434e0be0286dee61e06825
mmonsoor
parents:
diff changeset
351 medSam[which(medSam==0)] <- 0.0001
e13ec2c3fabe planemo upload commit 25fd6a739741295e3f434e0be0286dee61e06825
mmonsoor
parents:
diff changeset
352 mat <- apply(mat, 2, function(mat, medSam){
e13ec2c3fabe planemo upload commit 25fd6a739741295e3f434e0be0286dee61e06825
mmonsoor
parents:
diff changeset
353 medFDiSmpl <- mat/medSam
e13ec2c3fabe planemo upload commit 25fd6a739741295e3f434e0be0286dee61e06825
mmonsoor
parents:
diff changeset
354 vec<-mat/median(medFDiSmpl)
e13ec2c3fabe planemo upload commit 25fd6a739741295e3f434e0be0286dee61e06825
mmonsoor
parents:
diff changeset
355 return(vec)
e13ec2c3fabe planemo upload commit 25fd6a739741295e3f434e0be0286dee61e06825
mmonsoor
parents:
diff changeset
356 }, medSam)
e13ec2c3fabe planemo upload commit 25fd6a739741295e3f434e0be0286dee61e06825
mmonsoor
parents:
diff changeset
357 return (mat)
e13ec2c3fabe planemo upload commit 25fd6a739741295e3f434e0be0286dee61e06825
mmonsoor
parents:
diff changeset
358 }
e13ec2c3fabe planemo upload commit 25fd6a739741295e3f434e0be0286dee61e06825
mmonsoor
parents:
diff changeset
359 # this piece of code was copied from xcms
e13ec2c3fabe planemo upload commit 25fd6a739741295e3f434e0be0286dee61e06825
mmonsoor
parents:
diff changeset
360 pval <- function(X, classlabel, teststat) {
e13ec2c3fabe planemo upload commit 25fd6a739741295e3f434e0be0286dee61e06825
mmonsoor
parents:
diff changeset
361
e13ec2c3fabe planemo upload commit 25fd6a739741295e3f434e0be0286dee61e06825
mmonsoor
parents:
diff changeset
362 n1 <- rowSums(!is.na(X[,classlabel == 0]))
e13ec2c3fabe planemo upload commit 25fd6a739741295e3f434e0be0286dee61e06825
mmonsoor
parents:
diff changeset
363 n2 <- rowSums(!is.na(X[,classlabel == 1]))
e13ec2c3fabe planemo upload commit 25fd6a739741295e3f434e0be0286dee61e06825
mmonsoor
parents:
diff changeset
364 A <- apply(X[,classlabel == 0], 1, sd, na.rm=TRUE)^2/n1 ## sd(t(X[,classlabel == 0]), na.rm = TRUE)^2/n1
e13ec2c3fabe planemo upload commit 25fd6a739741295e3f434e0be0286dee61e06825
mmonsoor
parents:
diff changeset
365 B <- apply(X[,classlabel == 1], 1, sd, na.rm=TRUE)^2/n2 ## sd(t(X[,classlabel == 1]), na.rm = TRUE)^2/n2
e13ec2c3fabe planemo upload commit 25fd6a739741295e3f434e0be0286dee61e06825
mmonsoor
parents:
diff changeset
366 df <- (A+B)^2/(A^2/(n1-1)+B^2/(n2-1))
e13ec2c3fabe planemo upload commit 25fd6a739741295e3f434e0be0286dee61e06825
mmonsoor
parents:
diff changeset
367
e13ec2c3fabe planemo upload commit 25fd6a739741295e3f434e0be0286dee61e06825
mmonsoor
parents:
diff changeset
368 pvalue <- 2 * (1 - pt(abs(teststat), df))
e13ec2c3fabe planemo upload commit 25fd6a739741295e3f434e0be0286dee61e06825
mmonsoor
parents:
diff changeset
369 invisible(pvalue)
e13ec2c3fabe planemo upload commit 25fd6a739741295e3f434e0be0286dee61e06825
mmonsoor
parents:
diff changeset
370 }
e13ec2c3fabe planemo upload commit 25fd6a739741295e3f434e0be0286dee61e06825
mmonsoor
parents:
diff changeset
371
e13ec2c3fabe planemo upload commit 25fd6a739741295e3f434e0be0286dee61e06825
mmonsoor
parents:
diff changeset
372 c1 <- grep(class1, molIon$cameraobj@xcmsSet@phenoData[,1])
e13ec2c3fabe planemo upload commit 25fd6a739741295e3f434e0be0286dee61e06825
mmonsoor
parents:
diff changeset
373 c2 <- grep(class2, molIon$cameraobj@xcmsSet@phenoData[,1])
e13ec2c3fabe planemo upload commit 25fd6a739741295e3f434e0be0286dee61e06825
mmonsoor
parents:
diff changeset
374 testclab <- c(rep(0, length(c1)), rep(1, length(c2)))
e13ec2c3fabe planemo upload commit 25fd6a739741295e3f434e0be0286dee61e06825
mmonsoor
parents:
diff changeset
375 testval <- groupval(molIon$cameraobj@xcmsSet, "medret", "into")
e13ec2c3fabe planemo upload commit 25fd6a739741295e3f434e0be0286dee61e06825
mmonsoor
parents:
diff changeset
376 if(norm) testval <- normalize.medFC(testval)
e13ec2c3fabe planemo upload commit 25fd6a739741295e3f434e0be0286dee61e06825
mmonsoor
parents:
diff changeset
377 tstat <- mt.teststat(testval, testclab)
e13ec2c3fabe planemo upload commit 25fd6a739741295e3f434e0be0286dee61e06825
mmonsoor
parents:
diff changeset
378 pvalue <- pval(testval, testclab, tstat)
e13ec2c3fabe planemo upload commit 25fd6a739741295e3f434e0be0286dee61e06825
mmonsoor
parents:
diff changeset
379
e13ec2c3fabe planemo upload commit 25fd6a739741295e3f434e0be0286dee61e06825
mmonsoor
parents:
diff changeset
380 #
e13ec2c3fabe planemo upload commit 25fd6a739741295e3f434e0be0286dee61e06825
mmonsoor
parents:
diff changeset
381 # rport <- diffreport(molIon$cameraobj@xcmsSet, class1=class1, class2= class2, sortpval=FALSE)
e13ec2c3fabe planemo upload commit 25fd6a739741295e3f434e0be0286dee61e06825
mmonsoor
parents:
diff changeset
382 # ans[ans[,6]!="",5] <- rport[figidx, "pvalue"]
e13ec2c3fabe planemo upload commit 25fd6a739741295e3f434e0be0286dee61e06825
mmonsoor
parents:
diff changeset
383 ans[ans[,6]!="",5] <- pvalue[figidx]
e13ec2c3fabe planemo upload commit 25fd6a739741295e3f434e0be0286dee61e06825
mmonsoor
parents:
diff changeset
384 testname <- "t.test p-value"
e13ec2c3fabe planemo upload commit 25fd6a739741295e3f434e0be0286dee61e06825
mmonsoor
parents:
diff changeset
385 }
e13ec2c3fabe planemo upload commit 25fd6a739741295e3f434e0be0286dee61e06825
mmonsoor
parents:
diff changeset
386 if(m.test=="anova"){
e13ec2c3fabe planemo upload commit 25fd6a739741295e3f434e0be0286dee61e06825
mmonsoor
parents:
diff changeset
387 class <- molIon$cameraobj@xcmsSet@phenoData
e13ec2c3fabe planemo upload commit 25fd6a739741295e3f434e0be0286dee61e06825
mmonsoor
parents:
diff changeset
388 getPvalue <- function(dataidx) {
e13ec2c3fabe planemo upload commit 25fd6a739741295e3f434e0be0286dee61e06825
mmonsoor
parents:
diff changeset
389 aov.data <- data.frame(resp=as.numeric(peaklist[dataidx,coords]), class=class)
e13ec2c3fabe planemo upload commit 25fd6a739741295e3f434e0be0286dee61e06825
mmonsoor
parents:
diff changeset
390 anova(aov(resp~class, aov.data))$Pr[1]
e13ec2c3fabe planemo upload commit 25fd6a739741295e3f434e0be0286dee61e06825
mmonsoor
parents:
diff changeset
391 }
e13ec2c3fabe planemo upload commit 25fd6a739741295e3f434e0be0286dee61e06825
mmonsoor
parents:
diff changeset
392 testname <- "anova p-value"
e13ec2c3fabe planemo upload commit 25fd6a739741295e3f434e0be0286dee61e06825
mmonsoor
parents:
diff changeset
393 ans[ans[,6]!="",5] <- sapply(figidx, getPvalue)
e13ec2c3fabe planemo upload commit 25fd6a739741295e3f434e0be0286dee61e06825
mmonsoor
parents:
diff changeset
394 }
e13ec2c3fabe planemo upload commit 25fd6a739741295e3f434e0be0286dee61e06825
mmonsoor
parents:
diff changeset
395
e13ec2c3fabe planemo upload commit 25fd6a739741295e3f434e0be0286dee61e06825
mmonsoor
parents:
diff changeset
396 header <- matrix(c("Proposed Mass","Most probable Compound","Probability","Entropy", testname,"EIC-plot", "Ion annotation",coords), nrow=1 , ncol=7+length(coords) )
e13ec2c3fabe planemo upload commit 25fd6a739741295e3f434e0be0286dee61e06825
mmonsoor
parents:
diff changeset
397 ans <- rbind(header, ans)
e13ec2c3fabe planemo upload commit 25fd6a739741295e3f434e0be0286dee61e06825
mmonsoor
parents:
diff changeset
398
e13ec2c3fabe planemo upload commit 25fd6a739741295e3f434e0be0286dee61e06825
mmonsoor
parents:
diff changeset
399
e13ec2c3fabe planemo upload commit 25fd6a739741295e3f434e0be0286dee61e06825
mmonsoor
parents:
diff changeset
400 # additional field
e13ec2c3fabe planemo upload commit 25fd6a739741295e3f434e0be0286dee61e06825
mmonsoor
parents:
diff changeset
401 # ans <- cbind(ans[,1:2], ans[,2], ans[,3:ncol(ans)])
e13ec2c3fabe planemo upload commit 25fd6a739741295e3f434e0be0286dee61e06825
mmonsoor
parents:
diff changeset
402 #ans[ans[,3]!="unknown",][-1,3] <- as.character(DB$sbml.id[nid])
e13ec2c3fabe planemo upload commit 25fd6a739741295e3f434e0be0286dee61e06825
mmonsoor
parents:
diff changeset
403
e13ec2c3fabe planemo upload commit 25fd6a739741295e3f434e0be0286dee61e06825
mmonsoor
parents:
diff changeset
404 if(html) {
e13ec2c3fabe planemo upload commit 25fd6a739741295e3f434e0be0286dee61e06825
mmonsoor
parents:
diff changeset
405 #require(hwriter)
e13ec2c3fabe planemo upload commit 25fd6a739741295e3f434e0be0286dee61e06825
mmonsoor
parents:
diff changeset
406 ansb <- ans
e13ec2c3fabe planemo upload commit 25fd6a739741295e3f434e0be0286dee61e06825
mmonsoor
parents:
diff changeset
407 ans[ans[,2]!="unknown",][-1,2] <- as.character(DB$name[nid])
e13ec2c3fabe planemo upload commit 25fd6a739741295e3f434e0be0286dee61e06825
mmonsoor
parents:
diff changeset
408 if(linkPattern=="pubchem") ansb <- ans
e13ec2c3fabe planemo upload commit 25fd6a739741295e3f434e0be0286dee61e06825
mmonsoor
parents:
diff changeset
409
e13ec2c3fabe planemo upload commit 25fd6a739741295e3f434e0be0286dee61e06825
mmonsoor
parents:
diff changeset
410 hyper=matrix(paste(linkURL, ansb[-1,2], sep=""),ncol=1 )
e13ec2c3fabe planemo upload commit 25fd6a739741295e3f434e0be0286dee61e06825
mmonsoor
parents:
diff changeset
411 if(!is.null(molIon$cameraobj)) {
e13ec2c3fabe planemo upload commit 25fd6a739741295e3f434e0be0286dee61e06825
mmonsoor
parents:
diff changeset
412 hyper1=matrix(paste(fig, ans[-1,6],".png", sep=""),ncol=1 )
e13ec2c3fabe planemo upload commit 25fd6a739741295e3f434e0be0286dee61e06825
mmonsoor
parents:
diff changeset
413 }
e13ec2c3fabe planemo upload commit 25fd6a739741295e3f434e0be0286dee61e06825
mmonsoor
parents:
diff changeset
414 else {
e13ec2c3fabe planemo upload commit 25fd6a739741295e3f434e0be0286dee61e06825
mmonsoor
parents:
diff changeset
415 hyper1 <- ans[-1,6]
e13ec2c3fabe planemo upload commit 25fd6a739741295e3f434e0be0286dee61e06825
mmonsoor
parents:
diff changeset
416 hyper1[ans[-1,6]!=""][allidx[,2]!="neg"] <- paste(hyper1[ans[-1,6]!=""][allidx[,2]!="neg"], "pos", sep="")
e13ec2c3fabe planemo upload commit 25fd6a739741295e3f434e0be0286dee61e06825
mmonsoor
parents:
diff changeset
417 hyper1[ans[-1,6]!=""][allidx[,2]=="neg"] <- paste(hyper1[ans[-1,6]!=""][allidx[,2]=="neg"], "neg", sep="")
e13ec2c3fabe planemo upload commit 25fd6a739741295e3f434e0be0286dee61e06825
mmonsoor
parents:
diff changeset
418 hyper1=matrix(paste(fig, hyper1,".png", sep=""),ncol=1 )
e13ec2c3fabe planemo upload commit 25fd6a739741295e3f434e0be0286dee61e06825
mmonsoor
parents:
diff changeset
419 }
e13ec2c3fabe planemo upload commit 25fd6a739741295e3f434e0be0286dee61e06825
mmonsoor
parents:
diff changeset
420 p=openPage(paste(filename,".html",sep=""))
e13ec2c3fabe planemo upload commit 25fd6a739741295e3f434e0be0286dee61e06825
mmonsoor
parents:
diff changeset
421 ans2 <- ans[,1:7]
e13ec2c3fabe planemo upload commit 25fd6a739741295e3f434e0be0286dee61e06825
mmonsoor
parents:
diff changeset
422 link <- cbind(matrix(NA,nrow(ans2),1),rbind(NA,hyper),matrix(NA,nrow(ans2),3),rbind(NA,hyper1),matrix(NA,nrow(ans2),1))
e13ec2c3fabe planemo upload commit 25fd6a739741295e3f434e0be0286dee61e06825
mmonsoor
parents:
diff changeset
423 # This block is responsible to add as many columns as the user
e13ec2c3fabe planemo upload commit 25fd6a739741295e3f434e0be0286dee61e06825
mmonsoor
parents:
diff changeset
424 # wants
e13ec2c3fabe planemo upload commit 25fd6a739741295e3f434e0be0286dee61e06825
mmonsoor
parents:
diff changeset
425 if(!is.null(addLink)){
e13ec2c3fabe planemo upload commit 25fd6a739741295e3f434e0be0286dee61e06825
mmonsoor
parents:
diff changeset
426 for(l in 1:length(addLink$link)) {
e13ec2c3fabe planemo upload commit 25fd6a739741295e3f434e0be0286dee61e06825
mmonsoor
parents:
diff changeset
427 link <- cbind(link, rbind(NA, addLink[[l]]))
e13ec2c3fabe planemo upload commit 25fd6a739741295e3f434e0be0286dee61e06825
mmonsoor
parents:
diff changeset
428 }
e13ec2c3fabe planemo upload commit 25fd6a739741295e3f434e0be0286dee61e06825
mmonsoor
parents:
diff changeset
429 for(a in 1:length(addLink$ans)) {
e13ec2c3fabe planemo upload commit 25fd6a739741295e3f434e0be0286dee61e06825
mmonsoor
parents:
diff changeset
430 ans2 <- cbind(ans2,addLink$ans[[a]])
e13ec2c3fabe planemo upload commit 25fd6a739741295e3f434e0be0286dee61e06825
mmonsoor
parents:
diff changeset
431 #colnames(ans2)[7+a] <- addLink$ans[[a]][1]
e13ec2c3fabe planemo upload commit 25fd6a739741295e3f434e0be0286dee61e06825
mmonsoor
parents:
diff changeset
432 }
e13ec2c3fabe planemo upload commit 25fd6a739741295e3f434e0be0286dee61e06825
mmonsoor
parents:
diff changeset
433 }
e13ec2c3fabe planemo upload commit 25fd6a739741295e3f434e0be0286dee61e06825
mmonsoor
parents:
diff changeset
434 hwrite(ans2, p,row.bgcolor='#ffdc98', link=link )
e13ec2c3fabe planemo upload commit 25fd6a739741295e3f434e0be0286dee61e06825
mmonsoor
parents:
diff changeset
435 closePage(p)
e13ec2c3fabe planemo upload commit 25fd6a739741295e3f434e0be0286dee61e06825
mmonsoor
parents:
diff changeset
436 if(!is.null(molIon$cameraobj)) {
e13ec2c3fabe planemo upload commit 25fd6a739741295e3f434e0be0286dee61e06825
mmonsoor
parents:
diff changeset
437 plotEIC(molIon$cameraobj@xcmsSet, figidx, ans[ans[,6]!="",6][-1], colorplot=colorplot)
e13ec2c3fabe planemo upload commit 25fd6a739741295e3f434e0be0286dee61e06825
mmonsoor
parents:
diff changeset
438 }
e13ec2c3fabe planemo upload commit 25fd6a739741295e3f434e0be0286dee61e06825
mmonsoor
parents:
diff changeset
439 else {
e13ec2c3fabe planemo upload commit 25fd6a739741295e3f434e0be0286dee61e06825
mmonsoor
parents:
diff changeset
440 dataidxP <- as.numeric(allidx[allidx[,2]!="neg",1])
e13ec2c3fabe planemo upload commit 25fd6a739741295e3f434e0be0286dee61e06825
mmonsoor
parents:
diff changeset
441 pngidxP <- allidx[allidx[,2]!="neg",3]
e13ec2c3fabe planemo upload commit 25fd6a739741295e3f434e0be0286dee61e06825
mmonsoor
parents:
diff changeset
442 plotEIC(molIon$pos@xcmsSet, dataidxP, pngidxP, "pos", colorplot=colorplot)
e13ec2c3fabe planemo upload commit 25fd6a739741295e3f434e0be0286dee61e06825
mmonsoor
parents:
diff changeset
443 dataidxN <- as.numeric(allidx[allidx[,2]=="neg",1])
e13ec2c3fabe planemo upload commit 25fd6a739741295e3f434e0be0286dee61e06825
mmonsoor
parents:
diff changeset
444 pngidxN <- allidx[allidx[,2]=="neg",3]
e13ec2c3fabe planemo upload commit 25fd6a739741295e3f434e0be0286dee61e06825
mmonsoor
parents:
diff changeset
445 plotEIC(molIon$neg@xcmsSet, dataidxN, pngidxN, "neg", colorplot=colorplot)
e13ec2c3fabe planemo upload commit 25fd6a739741295e3f434e0be0286dee61e06825
mmonsoor
parents:
diff changeset
446 }
e13ec2c3fabe planemo upload commit 25fd6a739741295e3f434e0be0286dee61e06825
mmonsoor
parents:
diff changeset
447
e13ec2c3fabe planemo upload commit 25fd6a739741295e3f434e0be0286dee61e06825
mmonsoor
parents:
diff changeset
448 }
e13ec2c3fabe planemo upload commit 25fd6a739741295e3f434e0be0286dee61e06825
mmonsoor
parents:
diff changeset
449 else {
e13ec2c3fabe planemo upload commit 25fd6a739741295e3f434e0be0286dee61e06825
mmonsoor
parents:
diff changeset
450 ansb <- ans
e13ec2c3fabe planemo upload commit 25fd6a739741295e3f434e0be0286dee61e06825
mmonsoor
parents:
diff changeset
451 }
e13ec2c3fabe planemo upload commit 25fd6a739741295e3f434e0be0286dee61e06825
mmonsoor
parents:
diff changeset
452 colnames(ansb) <- ansb[1,]
e13ec2c3fabe planemo upload commit 25fd6a739741295e3f434e0be0286dee61e06825
mmonsoor
parents:
diff changeset
453 ansb <- ansb[-1,]
e13ec2c3fabe planemo upload commit 25fd6a739741295e3f434e0be0286dee61e06825
mmonsoor
parents:
diff changeset
454 return(list(classTable=ansb, figidx=figidx))
e13ec2c3fabe planemo upload commit 25fd6a739741295e3f434e0be0286dee61e06825
mmonsoor
parents:
diff changeset
455 }