Mercurial > repos > mmonsoor > probmetab

comparison lib.r @ 3:abcfa1648b66 draft

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
planemo upload commit c897279aa8cae0a4ad889bb05b143f32d2b6d712
author lecorguille
date Fri, 07 Apr 2017 07:14:12 -0400
parents e13ec2c3fabe
children 52b222a626b0
comparison
equal deleted inserted replaced
2:3a9914b37f2f 3:abcfa1648b66
3 # Contributors: Yann Guitton and Jean-francois Martin 3 # Contributors: Yann Guitton and Jean-francois Martin
4 4
5 5
6 ##Main probmetab function launch by the Galaxy ProbMetab wrapper 6 ##Main probmetab function launch by the Galaxy ProbMetab wrapper
7 probmetab = function(xa, xaP, xaN, variableMetadata, variableMetadataP, variableMetadataN, listArguments){ 7 probmetab = function(xa, xaP, xaN, variableMetadata, variableMetadataP, variableMetadataN, listArguments){
8 ##ONE MODE ACQUISITION## 8 ##ONE MODE ACQUISITION##
9 if(listArguments[["mode_acquisition"]]=="one") { 9 if(listArguments[["mode_acquisition"]]=="one") {
10 comb=NULL 10 comb=NULL
11 11
12 #Get the polarity from xa object 12 #Get the polarity from xa object
13 polarity=xa@polarity 13 polarity=xa@polarity
14 #SNR option 14 #SNR option
15 if ("xsetnofill" %in% names(listArguments)) { 15 if ("xsetnofill" %in% names(listArguments)) {
16 load(listArguments[["xsetnofill"]]) 16 load(listArguments[["xsetnofill"]])
17 xsetnofill=xset 17 xsetnofill=xset
18 } 18 }
19 else{ 19 else{
20 xsetnofill=NULL 20 xsetnofill=NULL
21 } 21 }
22 #Exclude samples 22 #Exclude samples
23 if ("toexclude" %in% names(listArguments)) { 23 if ("toexclude" %in% names(listArguments)) {
24 toexclude=listArguments[["toexclude"]] 24 toexclude=listArguments[["toexclude"]]
25 } 25 }
26 else { 26 else {
27 toexclude=NULL 27 toexclude=NULL
28 } 28 }
29 ionAnnot=get.annot(xa, polarity=polarity, allowMiss=listArguments[["allowMiss"]],xset=xsetnofill,toexclude=toexclude) 29 ionAnnot=get.annot(xa, polarity=polarity, allowMiss=listArguments[["allowMiss"]],xset=xsetnofill,toexclude=toexclude)
30 comb=NULL 30 comb=NULL
31 } 31 }
32 32
33 ##TWO MODES ACQUISITION## 33 ##TWO MODES ACQUISITION##
34 #Mode annotatediffreport 34 #Mode annotatediffreport
35 else if(listArguments[["inputs_mode"]]=="two"){ 35 else if(listArguments[["inputs_mode"]]=="two"){
36 ##Prepare the objects that will be used for the get.annot function 36 ##Prepare the objects that will be used for the get.annot function
37 comb=1 37 comb=1
38 38
39 39
40 xsetPnofill=NULL 40 xsetPnofill=NULL
41 xsetNnofill=NULL 41 xsetNnofill=NULL
42 # TODO: a reactiver 42 # TODO: a reactiver
43 #if ("xsetPnofill" %in% names(listArguments)) { 43 #if ("xsetPnofill" %in% names(listArguments)) {
44 # load(listArguments[["xsetPnofill"]]) 44 # load(listArguments[["xsetPnofill"]])
45 # xsetPnofill=xset 45 # xsetPnofill=xset
46 #} 46 #}
47 #if ("xsetNnofill" %in% names(listArguments)) { 47 #if ("xsetNnofill" %in% names(listArguments)) {
48 # load(listArguments[["xsetNnofill"]]) 48 # load(listArguments[["xsetNnofill"]])
49 # xsetNnofill=xset 49 # xsetNnofill=xset
50 #} 50 #}
51 # include CAMERA non-annotated compounds, and snr retrieval 51 # include CAMERA non-annotated compounds, and snr retrieval
52 # comb 2+ - used on Table 1 52 # comb 2+ - used on Table 1
53 ionAnnotP2plus = get.annot(axP, allowMiss=listArguments[["allowMiss"]], xset=xsetPnofill,toexclude=listArguments[["toexclude"]]) 53 ionAnnotP2plus = get.annot(xaP, allowMiss=listArguments[["allowMiss"]], xset=xsetPnofill,toexclude=listArguments[["toexclude"]])
54 ionAnnotN2plus = get.annot(axN, polarity="negative", allowMiss=listArguments[["allowMiss"]], xset=xsetNnofill,toexclude=listArguments[["toexclude"]]) 54 ionAnnotN2plus = get.annot(xaN, polarity="negative", allowMiss=listArguments[["allowMiss"]], xset=xsetNnofill,toexclude=listArguments[["toexclude"]])
55 ionAnnot = combineMolIon(ionAnnotP2plus, ionAnnotN2plus) 55 ionAnnot = combineMolIon(ionAnnotP2plus, ionAnnotN2plus)
56 print(sum(ionAnnot$molIon[,3]==1)) 56 print(sum(ionAnnot$molIon[,3]==1))
57 print(sum(ionAnnot$molIon[,3]==0)) 57 print(sum(ionAnnot$molIon[,3]==0))
58 write.table(ionAnnot[1], sep="\t", quote=FALSE, row.names=FALSE, file="CombineMolIon.tsv") 58 write.table(ionAnnot[1], sep="\t", quote=FALSE, row.names=FALSE, file="CombineMolIon.tsv")
59 #Merge variableMetadata Negative and positive acquisitions mode 59 #Merge variableMetadata Negative and positive acquisitions mode
60 60
61 61
62 #Mode combinexsannos TODO bug avec tableau issus de combinexsannos 62 #Mode combinexsannos TODO bug avec tableau issus de combinexsannos
63 #else { 63 #else {
64 #load(listArguments[["image_combinexsannos"]]) 64 #load(listArguments[["image_combinexsannos"]])
65 #image_combinexsannos=cAnnot 65 #image_combinexsannos=cAnnot
66 ##Prepare the objects that will be used for the combineMolIon function 66 ##Prepare the objects that will be used for the combineMolIon function
67 #load(listArguments[["image_pos"]]) 67 #load(listArguments[["image_pos"]])
68 #image_pos=xa 68 #image_pos=xa
69 #ionAnnot=combineMolIon(peaklist=cAnnot, cameraobj=image_pos, polarity="pos") 69 #ionAnnot=combineMolIon(peaklist=cAnnot, cameraobj=image_pos, polarity="pos")
70 #} 70 #}
71 71
72 } 72 }
73 73
74 ##DATABASE MATCHING## 74 ##DATABASE MATCHING##
75 if (listArguments[["kegg_db"]]=="KEGG"){ 75 if (listArguments[["kegg_db"]]=="KEGG"){
76 DB=build.database.kegg(orgID = NULL) 76 DB=build.database.kegg(orgID = NULL)
77 } 77 }
78 else{ 78 else{
79 table_list <<- NULL 79 table_list <<- NULL
80 ids=strsplit(listArguments[["kegg_db"]],",") 80 ids=strsplit(listArguments[["kegg_db"]],",")
81 ids=ids[[1]] 81 ids=ids[[1]]
82 if(length(ids)>1){ 82 if(length(ids)>1){
83 for(i in 1:length(ids)){ 83 for(i in 1:length(ids)){
84 table_list[[i]] <- build.database.kegg(ids[i]) 84 table_list[[i]] <- build.database.kegg(ids[i])
85 } 85 }
86 db_table=do.call("rbind",table_list) 86 db_table=do.call("rbind",table_list)
87 DB=unique(db_table) 87 DB=unique(db_table)
88 } 88 }
89 else{ 89 else{
90 DB=build.database.kegg(listArguments[["kegg_db"]]) 90 DB=build.database.kegg(listArguments[["kegg_db"]])
91 } 91 }
92 } 92 }
93 #Matching des mass exactes mesurees avec les masses des compounds KEGG (pas M+H ou M-H) 93 #Matching des mass exactes mesurees avec les masses des compounds KEGG (pas M+H ou M-H)
94 reactionM = create.reactionM(DB, ionAnnot, ppm.tol=listArguments[["ppm_tol"]]) 94 reactionM = create.reactionM(DB, ionAnnot, ppm.tol=listArguments[["ppm_tol"]])
95 ##PROBABILITY RANKING## 95 ##PROBABILITY RANKING##
96 # number of masses with candidates inside the fixed mass window 96 # number of masses with candidates inside the fixed mass window
97 # and masses with more than one candidate 97 # and masses with more than one candidate
98 length(unique(reactionM[reactionM[,"id"]!="unknown",1])) 98 length(unique(reactionM[reactionM[,"id"]!="unknown",1]))
99 sum(table(reactionM[reactionM[,"id"]!="unknown",1])>1) 99 sum(table(reactionM[reactionM[,"id"]!="unknown",1])>1)
100 #if (listArguments[["useIso"]]){ 100 #if (listArguments[["useIso"]]){
101 #BUG TODO 101 #BUG TODO
102 # Calculate the ratio between observed and theoretical isotopic patterns. 102 # Calculate the ratio between observed and theoretical isotopic patterns.
103 # If you don't have an assessment of carbon offset to carbon number prediction 103 # If you don't have an assessment of carbon offset to carbon number prediction
104 # skip this step and use the reactionM as input to weigthM function. 104 # skip this step and use the reactionM as input to weigthM function.
105 #isoPatt < incorporate.isotopes(comb2plus, reactionM, , samp=12:23, DB=DB) 105 #isoPatt < incorporate.isotopes(comb2plus, reactionM, , samp=12:23, DB=DB)
106 # calculate the likelihood of each mass to compound assignment using mass accuracy,and isotopic pattern, when present 106 # calculate the likelihood of each mass to compound assignment using mass accuracy,and isotopic pattern, when present
107 #wl < weightM(isoPatt,intervals=seq(0,1000,by=500), offset=c(3.115712, 3.434146, 2.350798)) 107 #wl < weightM(isoPatt,intervals=seq(0,1000,by=500), offset=c(3.115712, 3.434146, 2.350798))
108 108
109 #isoPatt=incorporate.isotopes(ionAnnot, reactionM,comb=comb,var=listArguments[["var"]],DB=DB) 109 #isoPatt=incorporate.isotopes(ionAnnot, reactionM,comb=comb,var=listArguments[["var"]],DB=DB)
110 110
111 #wl = weightM(reactionM, useIso=true) 111 #wl = weightM(reactionM, useIso=true)
112 #} 112 #}
113 #else { 113 #else {
114 #wl = weightM(reactionM, useIso=FALSE) 114 #wl = weightM(reactionM, useIso=FALSE)
115 #} 115 #}
116 wl =weightM(reactionM, useIso=FALSE) 116 wl =weightM(reactionM, useIso=FALSE)
117 w = design.connection(reactionM) 117 w = design.connection(reactionM)
118 # Probability calculations 118 # Probability calculations
119 x = 1:ncol(wl$wm) 119 x = 1:ncol(wl$wm)
120 y = 1:nrow(wl$wm) 120 y = 1:nrow(wl$wm)
121 conn = gibbs.samp(x, y, 5000, w, wl$wm) 121 conn = gibbs.samp(x, y, 5000, w, wl$wm)
122 ansConn = export.class.table(conn, reactionM, ionAnnot, DB=DB,html=listArguments[["html"]],filename="AnalysisExample",prob=listArguments[["prob"]]) 122 ansConn = export.class.table(conn, reactionM, ionAnnot, DB=DB,html=listArguments[["html"]],filename="AnalysisExample",prob=listArguments[["prob"]])
123 if(listArguments[["html"]]){ 123 if(listArguments[["html"]]){
124 #Zip the EICS plot 124 #Zip the EICS plot
125 system(paste('zip -r "Analysis_Report.zip" "AnalysisExample_fig"')) 125 system(paste('zip -r "Analysis_Report.zip" "AnalysisExample_fig"'))
126 } 126 }
127 127
128 # calculate the correlations and partial correlations and cross reference then with reactions 128 # calculate the correlations and partial correlations and cross reference then with reactions
129 mw=which(w==1,arr.ind=TRUE) 129 mw=which(w==1,arr.ind=TRUE)
130 #reac2cor function : Use the intensity of putative molecules in repeated samples to calculate correlations and partial 130 #reac2cor function : Use the intensity of putative molecules in repeated samples to calculate correlations and partial
131 #correlation in a user defined threshold of false discovery rate for significance testing. After the 131 #correlation in a user defined threshold of false discovery rate for significance testing. After the
132 #correlation test the function also overlay significant correlations with all putative reactions between 132 #correlation test the function also overlay significant correlations with all putative reactions between
133 #two masses. 133 #two masses.
134 #It generates a list of estimated correlations and reactions. 134 #It generates a list of estimated correlations and reactions.
135 corList=reac2cor(mw,ansConn$classTable,listArguments[["opt"]],listArguments[["corths"]],listArguments[["corprob"]],listArguments[["pcorprob"]]) 135 corList=reac2cor(mw,ansConn$classTable,listArguments[["opt"]],listArguments[["corths"]],listArguments[["corprob"]],listArguments[["pcorprob"]])
136 ans=list("ansConn"=ansConn,"corList"=corList) 136 ans=list("ansConn"=ansConn,"corList"=corList)
137 #Generate the siff table for CytoScape 137 #Generate the siff table for CytoScape
138 cytoscape_output(corList,ansConn) 138 cytoscape_output(corList,ansConn)
139 139
140 140
141 #Execute the merge_probmetab function to merge the variableMetadata table and annotations from ProbMetab results 141 #Execute the merge_probmetab function to merge the variableMetadata table and annotations from ProbMetab results
142 142
143 if(listArguments[["mode_acquisition"]]=="one") { 143 if(listArguments[["mode_acquisition"]]=="one") {
144 #Retrocompatibility with previous annotateDiffreport variableMetadata dataframe (must replace mzmed column by mz, and rtmed by rt) 144 #Retrocompatibility with previous annotateDiffreport variableMetadata dataframe (must replace mzmed column by mz, and rtmed by rt)
145 names(variableMetadata)[names(variableMetadata)=="mzmed"] <- "mz" 145 names(variableMetadata)[names(variableMetadata)=="mzmed"] <- "mz"
146 names(variableMetadata)[names(variableMetadata)=="rtmed"] <- "rt" 146 names(variableMetadata)[names(variableMetadata)=="rtmed"] <- "rt"
147 variableM=merge_probmetab(variableMetadata, ansConn) 147 variableM=merge_probmetab(variableMetadata, ansConn)
148 write.table(variableM, sep="\t", quote=FALSE, row.names=FALSE, file="variableMetadata.tsv") 148 write.table(variableM, sep="\t", quote=FALSE, row.names=FALSE, file="variableMetadata.tsv")
149 } else if (listArguments[["mode_acquisition"]]=="two") { 149 } else if (listArguments[["mode_acquisition"]]=="two") {
150 #Retrocompatibility with previous annotateDiffreport variableMetadata dataframe (must replace mzmed column by mz, and rtmed by rt) 150 #Retrocompatibility with previous annotateDiffreport variableMetadata dataframe (must replace mzmed column by mz, and rtmed by rt)
151 names(variableMetadataP)[names(variableMetadata)=="mzmed"] <- "mz" 151 names(variableMetadataP)[names(variableMetadata)=="mzmed"] <- "mz"
152 names(variableMetadataP)[names(variableMetadata)=="rtmed"] <- "rt" 152 names(variableMetadataP)[names(variableMetadata)=="rtmed"] <- "rt"
153 names(variableMetadataN)[names(variableMetadata)=="mzmed"] <- "mz" 153 names(variableMetadataN)[names(variableMetadata)=="mzmed"] <- "mz"
154 names(variableMetadataN)[names(variableMetadata)=="rtmed"] <- "rt" 154 names(variableMetadataN)[names(variableMetadata)=="rtmed"] <- "rt"
155 variableMP=merge_probmetab(variableMetadataP, ansConn) 155 variableMP=merge_probmetab(variableMetadataP, ansConn)
156 write.table(variableMP, sep="\t", quote=FALSE, row.names=FALSE, file="variableMetadata_Positive.tsv") 156 write.table(variableMP, sep="\t", quote=FALSE, row.names=FALSE, file="variableMetadata_Positive.tsv")
157 variableMN=merge_probmetab(variableMetadataN, ansConn) 157 variableMN=merge_probmetab(variableMetadataN, ansConn)
158 write.table(variableMN, sep="\t", quote=FALSE, row.names=FALSE, file="variableMetadata_Negative.tsv") 158 write.table(variableMN, sep="\t", quote=FALSE, row.names=FALSE, file="variableMetadata_Negative.tsv")
159 } 159 }
160 160
161 return(ans) 161 return(ans)
162 162
163 } 163 }
164 164
165 ##Function that generates a siff table for CytoScape 165 ##Function that generates a siff table for CytoScape
166 cytoscape_output=function(corList,ansConn){ 166 cytoscape_output=function(corList,ansConn){
167 signif_cor=as.data.frame(corList$signif.cor) 167 signif_cor=as.data.frame(corList$signif.cor)
168 classTable=as.data.frame(ansConn$classTable) 168 classTable=as.data.frame(ansConn$classTable)
169 #Siff table 169 #Siff table
170 siff_table=cbind(signif_cor["node1"],signif_cor["cor"],signif_cor["node2"]) 170 siff_table=cbind(signif_cor["node1"],signif_cor["cor"],signif_cor["node2"])
171 #attribute table output for Cytoscape 171 #attribute table output for Cytoscape
172 172
173 ## START Code part from the export2cytoscape function of ProbMetab written by Ricardo R. Silva 173 ## START Code part from the export2cytoscape function of ProbMetab written by Ricardo R. Silva
174 for (i in 1:nrow(classTable)) if (classTable[i, 1] == ""){ 174 for (i in 1:nrow(classTable)) if (classTable[i, 1] == ""){
175 classTable[i, c(1, 4, 6, 7)] <- classTable[i - 1, c(1, 4, 6, 7)] 175 classTable[i, c(1, 4, 6, 7)] <- classTable[i - 1, c(1, 4, 6, 7)]
176 } 176 }
177 msel <- as.matrix(classTable[, 1:7]) 177 msel <- as.matrix(classTable[, 1:7])
178 msel <- cbind(msel[, 6], msel[,-6]) 178 msel <- cbind(msel[, 6], msel[,-6])
179 colnames(msel)[1] <- "Id" 179 colnames(msel)[1] <- "Id"
180 msel[, 1] <- sub("^\\s+", "", msel[, 1]) 180 msel[, 1] <- sub("^\\s+", "", msel[, 1])
181 colnames(msel)[1] <- "Id" 181 colnames(msel)[1] <- "Id"
182 ids <- unique(msel[, 1]) 182 ids <- unique(msel[, 1])
183 attrMatrix <- matrix("", nrow = length(ids), ncol = ncol(msel)-1) 183 attrMatrix <- matrix("", nrow = length(ids), ncol = ncol(msel)-1)
184 for (i in 1:length(ids)) { 184 for (i in 1:length(ids)) {
185 attrMatrix[i, 1] <- unique(msel[msel[, 1] == ids[i], 185 attrMatrix[i, 1] <- unique(msel[msel[, 1] == ids[i],
186 2]) 186 2])
187 attrMatrix[i, 2] <- paste("[", paste(msel[msel[, 187 attrMatrix[i, 2] <- paste("[", paste(msel[msel[,
188 1] == ids[i], 3], collapse = ", "), "]", sep = "") 188 1] == ids[i], 3], collapse = ", "), "]", sep = "")
189 attrMatrix[i, 3] <- paste("[", paste(msel[msel[, 189 attrMatrix[i, 3] <- paste("[", paste(msel[msel[,
190 1] == ids[i], 4], collapse = ", "), "]", sep = "") 190 1] == ids[i], 4], collapse = ", "), "]", sep = "")
191 attrMatrix[i, 4] <- unique(msel[msel[, 1] == ids[i], 191 attrMatrix[i, 4] <- unique(msel[msel[, 1] == ids[i],
192 5]) 192 5])
193 attrMatrix[i, 5] <- paste("[", paste(msel[msel[, 193 attrMatrix[i, 5] <- paste("[", paste(msel[msel[,
194 1] == ids[i], 6], collapse = ", "), "]", sep = "") 194 1] == ids[i], 6], collapse = ", "), "]", sep = "")
195 attrMatrix[i, 6] <- unique(msel[msel[, 1] == ids[i], 195 attrMatrix[i, 6] <- unique(msel[msel[, 1] == ids[i],
196 7]) 196 7])
197 } 197 }
198 ids <- as.numeric(unique(msel[, 1])) 198 ids <- as.numeric(unique(msel[, 1]))
199 attrMatrix <- cbind(ids, attrMatrix) 199 attrMatrix <- cbind(ids, attrMatrix)
200 colnames(attrMatrix) <- colnames(msel) 200 colnames(attrMatrix) <- colnames(msel)
201 ## END Code part from the export2cytoscape function of ProbMetab writieen by Ricardo R. Silva 201 ## END Code part from the export2cytoscape function of ProbMetab writieen by Ricardo R. Silva
202 write.table(attrMatrix, sep="\t", quote=FALSE, row.names=FALSE, file="Analysis_Report.tsv") 202 write.table(attrMatrix, sep="\t", quote=FALSE, row.names=FALSE, file="Analysis_Report.tsv")
203 write.table(siff_table, sep="\t", quote=FALSE, row.names=FALSE, file="sif.tsv") 203 write.table(siff_table, sep="\t", quote=FALSE, row.names=FALSE, file="sif.tsv")
204 204
205 return(attrMatrix) 205 return(attrMatrix)
206 } 206 }
207 207
208 ##Functions written by Jean-Francois Martin 208 ##Functions written by Jean-Francois Martin
209 209
210 deter_ioni <- function (aninfo, pm) 210 deter_ioni <- function (aninfo, pm)
211 { 211 {
212 # determine ionisation in ProbMetab result file, used in function merge_probmetab 212 # determine ionisation in ProbMetab result file, used in function merge_probmetab
213 # input : for 1 ion, aninfo = string with m/z rt and CAMERA annotation from ProbMetab result file 213 # input : for 1 ion, aninfo = string with m/z rt and CAMERA annotation from ProbMetab result file
214 # if the difference between m/z and the probmetab proposed mass is ~1 we use the sign (positive or negative) of this diference 214 # if the difference between m/z and the probmetab proposed mass is ~1 we use the sign (positive or negative) of this diference
215 # to define the type of ionisation 215 # to define the type of ionisation
216 # If adduct or fragment was detected, therefore diff >>1 and so, we search for substring "+" ou "2+" ou "3+" ou "-"... 216 # If adduct or fragment was detected, therefore diff >>1 and so, we search for substring "+" ou "2+" ou "3+" ou "-"...
217 # to define the type of ionisation 217 # to define the type of ionisation
218 # aninfo : vecteur of character resulting of the parsing(sep="#") of the probmetab annotation 218 # aninfo : vecteur of character resulting of the parsing(sep="#") of the probmetab annotation
219 if (round(abs(as.numeric(aninfo[1]) - pm),0) ==1) { 219 if (round(abs(as.numeric(aninfo[1]) - pm),0) ==1) {
220 if (as.numeric(aninfo[1]) - pm <0) {esi <- "n"} else {esi <- "p"} 220 if (as.numeric(aninfo[1]) - pm <0) {esi <- "n"} else {esi <- "p"}
221 } else 221 } else
222 if (!is.na(aninfo[4])) { 222 if (!is.na(aninfo[4])) {
223 anstr <- aninfo[4] 223 anstr <- aninfo[4]
224 # cat(anstr) 224 # cat(anstr)
225 if ((grepl("]+",anstr,fixed=T)==T) || (grepl("]2+",anstr,fixed=T)==T) || (grepl("]3+",anstr,fixed=T)==T)) { esi <- "p"} 225 if ((grepl("]+",anstr,fixed=T)==T) || (grepl("]2+",anstr,fixed=T)==T) || (grepl("]3+",anstr,fixed=T)==T)) { esi <- "p"}
226 else 226 else
227 if ((grepl("]-",anstr,fixed=T)==T) || (grepl("]2-",anstr,fixed=T)==T) || (grepl("]3-",anstr,fixed=T)==T)) { esi <- "n"} 227 if ((grepl("]-",anstr,fixed=T)==T) || (grepl("]2-",anstr,fixed=T)==T) || (grepl("]3-",anstr,fixed=T)==T)) { esi <- "n"}
228 # cat(" ioni ",esi,"\n") 228 # cat(" ioni ",esi,"\n")
229 } else 229 } else
230 { esi <- "u"} 230 { esi <- "u"}
231 231
232 return(esi) 232 return(esi)
233 } 233 }
234 234
235 235
236 merge_probmetab <- function(metaVar,ansConn) { 236 merge_probmetab <- function(metaVar,ansConn) {
237 ## Parse ProbMetab information result file and merge in variable_metaData initial file 237 ## Parse ProbMetab information result file and merge in variable_metaData initial file
238 ## inputs : 238 ## inputs :
239 ## metaVar : data.frame of metadataVariable input of probmetab function 239 ## metaVar : data.frame of metadataVariable input of probmetab function
240 ## ansConn : data.frame of ProbMetab result 240 ## ansConn : data.frame of ProbMetab result
241 ## output : dataframe with Probmetab results merge with variableMetadata 241 ## output : dataframe with Probmetab results merge with variableMetadata
242 ## Constante 242 ## Constante
243 ## iannot : indice de la colonne annotation dans le resultat de probMetab 243 ## iannot : indice de la colonne annotation dans le resultat de probMetab
244 iannot <- 4 244 iannot <- 4
245 245
246 ## definition of an unique identification of ions mz with 3 decimals and rt(sec) with 1 decimal to avoid 246 ## definition of an unique identification of ions mz with 3 decimals and rt(sec) with 1 decimal to avoid
247 ## duplicate ions name in the diffreport result file 247 ## duplicate ions name in the diffreport result file
248 ions <- paste ("M",round(metaVar$mz,3),"T",round(metaVar$rt,1),sep="") 248 ions <- paste ("M",round(metaVar$mz,3),"T",round(metaVar$rt,1),sep="")
249 metaVar <- data.frame(ions,metaVar) 249 metaVar <- data.frame(ions,metaVar)
250 250
251 ###### Result data.frame from ProbMetab result list 251 ###### Result data.frame from ProbMetab result list
252 an_ini <- ansConn$classTable 252 an_ini <- ansConn$classTable
253 253
254 ## Suppression of rows without mz and rt or unknown and columns of intensities 254 ## Suppression of rows without mz and rt or unknown and columns of intensities
255 ## COLUMNS SUBSCRIPTS HAVE TO BE CHECKED WITh DIFFERENT RESULTS FILES 255 ## COLUMNS SUBSCRIPTS HAVE TO BE CHECKED WITh DIFFERENT RESULTS FILES
256 an <- an_ini[(an_ini[,2]!="unknown"),c(1,2,3,7)] 256 an <- an_ini[(an_ini[,2]!="unknown"),c(1,2,3,7)]
257 ## initialisation of vectors receiving the result of the parse of the column annotation (subscrip iannot) 257 ## initialisation of vectors receiving the result of the parse of the column annotation (subscrip iannot)
258 mz <- rep(0,dim(an)[1]) 258 mz <- rep(0,dim(an)[1])
259 rt <- rep(0,dim(an)[1]) 259 rt <- rep(0,dim(an)[1])
260 propmz <- rep(0,dim(an)[1]) 260 propmz <- rep(0,dim(an)[1])
261 ioni <- rep("u",dim(an)[1]) 261 ioni <- rep("u",dim(an)[1])
262 262
263 ## parse the column annotation and define ionisation mode 263 ## parse the column annotation and define ionisation mode
264 for (i in 1:dim(an)[1]) { 264 for (i in 1:dim(an)[1]) {
265 if (an[i,1] != "") { 265 if (an[i,1] != "") {
266 info_mzrt <- unlist(strsplit(an[i,iannot],"#")) 266 info_mzrt <- unlist(strsplit(an[i,iannot],"#"))
267 propmz[i] <- as.numeric(an[i,1]) 267 propmz[i] <- as.numeric(an[i,1])
274 mz[i] <- as.numeric(mz[i-1]) 274 mz[i] <- as.numeric(mz[i-1])
275 rt[i] <- as.numeric(rt[i-1]) 275 rt[i] <- as.numeric(rt[i-1])
276 ioni[i] <- ioni[i-1] 276 ioni[i] <- ioni[i-1]
277 } 277 }
278 } 278 }
279 279
280 ## definition of an unique identification of ions : mz with 3 decimals and rt(sec) with 1 decimal 280 ## definition of an unique identification of ions : mz with 3 decimals and rt(sec) with 1 decimal
281 ## The same as for the metadataVariable data.frame to match with. 281 ## The same as for the metadataVariable data.frame to match with.
282 ions <- paste ("M",round(mz,3),"T",round(rt,1),sep="") 282 ions <- paste ("M",round(mz,3),"T",round(rt,1),sep="")
283 an <- data.frame(ions,ioni,propmz,mz,rt,an) 283 an <- data.frame(ions,ioni,propmz,mz,rt,an)
284 284
285 ## transposition of the different probmetab annotations which are in different rows in the initial result data.frame 285 ## transposition of the different probmetab annotations which are in different rows in the initial result data.frame
286 ## on only 1 row separated with a ";" 286 ## on only 1 row separated with a ";"
287 li <- as.matrix(table(an$propmz)) 287 li <- as.matrix(table(an$propmz))
288 li <- data.frame(dimnames(li)[1],li) 288 li <- data.frame(dimnames(li)[1],li)
289 dimnames(li)[[2]][1] <- "propmz" 289 dimnames(li)[[2]][1] <- "propmz"
296 c <- c + 1 296 c <- c + 1
297 suban <- an[an$propmz==li[c,1],] 297 suban <- an[an$propmz==li[c,1],]
298 ions[c] <- as.character(suban[1,1]) 298 ions[c] <- as.character(suban[1,1])
299 propmz[c] <- as.numeric(suban[1,3]) 299 propmz[c] <- as.numeric(suban[1,3])
300 mpc[c] <- paste(suban[,7],collapse=";") 300 mpc[c] <- paste(suban[,7],collapse=";")
301 proba[c] <- paste(as.character(suban[,8]),collapse=";") 301 proba[c] <- paste(as.character(suban[,8]),collapse=";")
302 } 302 }
303 303
304 ## Creation of the data.frame with 1 row per ions 304 ## Creation of the data.frame with 1 row per ions
305 anc <- data.frame(ions,propmz,mpc,proba) 305 anc <- data.frame(ions,propmz,mpc,proba)
306 anc <- anc[order(anc[,1]),] 306 anc <- anc[order(anc[,1]),]
307 307
308 metaVarFinal <- merge(metaVar, anc, by.x=1, by.y=1, all.x=T, all.y=T) 308 metaVarFinal <- merge(metaVar, anc, by.x=1, by.y=1, all.x=T, all.y=T)
309 metaVarFinal <- metaVarFinal[,-1] 309 metaVarFinal <- metaVarFinal[,-1]
310 #write.table(metaVarFinal,file="res.txt", sep="\t", row.names=F, quote=F) 310 #write.table(metaVarFinal,file="res.txt", sep="\t", row.names=F, quote=F)
311 311
312 return (metaVarFinal) 312 return (metaVarFinal)
313 } 313 }
314 314
315 # RETROCOMPATIBILITE avec ancienne version de annotate 315 # RETROCOMPATIBILITE avec ancienne version de annotate
316 getVariableMetadata = function(xa) { 316 getVariableMetadata = function(xa) {
317 # --- variableMetadata --- 317 # --- variableMetadata ---
318 peakList=getPeaklist(xa) 318 peakList=getPeaklist(xa)
319 peakList=cbind(groupnames(xa@xcmsSet),peakList); colnames(peakList)[1] = c("name"); 319 peakList=cbind(groupnames(xa@xcmsSet),peakList); colnames(peakList)[1] = c("name");
320 variableMetadata=peakList[,!(colnames(peakList) %in% c(sampnames(xa@xcmsSet)))] 320 variableMetadata=peakList[,!(colnames(peakList) %in% c(sampnames(xa@xcmsSet)))]
321 variableMetadata$name= paste("M",round(variableMetadata$mz),"T",round(variableMetadata$rt),sep="") 321 variableMetadata$name= groupnames(xa@xcmsSet)
322 return (variableMetadata) 322 return (variableMetadata)
323 } 323 }
324
325
326 # This function get the raw file path from the arguments
327 getRawfilePathFromArguments <- function(singlefile, zipfile, listArguments) {
328 if (!is.null(listArguments[["zipfile"]])) zipfile = listArguments[["zipfile"]]
329 if (!is.null(listArguments[["zipfilePositive"]])) zipfile = listArguments[["zipfilePositive"]]
330 if (!is.null(listArguments[["zipfileNegative"]])) zipfile = listArguments[["zipfileNegative"]]
331
332 if (!is.null(listArguments[["singlefile_galaxyPath"]])) {
333 singlefile_galaxyPaths = listArguments[["singlefile_galaxyPath"]];
334 singlefile_sampleNames = listArguments[["singlefile_sampleName"]]
335 }
336 if (!is.null(listArguments[["singlefile_galaxyPathPositive"]])) {
337 singlefile_galaxyPaths = listArguments[["singlefile_galaxyPathPositive"]];
338 singlefile_sampleNames = listArguments[["singlefile_sampleNamePositive"]]
339 }
340 if (!is.null(listArguments[["singlefile_galaxyPathNegative"]])) {
341 singlefile_galaxyPaths = listArguments[["singlefile_galaxyPathNegative"]];
342 singlefile_sampleNames = listArguments[["singlefile_sampleNameNegative"]]
343 }
344 if (exists("singlefile_galaxyPaths")){
345 singlefile_galaxyPaths = unlist(strsplit(singlefile_galaxyPaths,","))
346 singlefile_sampleNames = unlist(strsplit(singlefile_sampleNames,","))
347
348 singlefile=NULL
349 for (singlefile_galaxyPath_i in seq(1:length(singlefile_galaxyPaths))) {
350 singlefile_galaxyPath=singlefile_galaxyPaths[singlefile_galaxyPath_i]
351 singlefile_sampleName=singlefile_sampleNames[singlefile_galaxyPath_i]
352 singlefile[[singlefile_sampleName]] = singlefile_galaxyPath
353 }
354 }
355 return(list(zipfile=zipfile, singlefile=singlefile))
356 }
357
358
359 # This function retrieve the raw file in the working directory
360 # - if zipfile: unzip the file with its directory tree
361 # - if singlefiles: set symlink with the good filename
362 retrieveRawfileInTheWorkingDirectory <- function(singlefile, zipfile) {
363
364 if(!is.null(singlefile) && (length("singlefile")>0)) {
365 for (singlefile_sampleName in names(singlefile)) {
366 singlefile_galaxyPath = singlefile[[singlefile_sampleName]]
367 if(!file.exists(singlefile_galaxyPath)){
368 error_message=paste("Cannot access the sample:",singlefile_sampleName,"located:",singlefile_galaxyPath,". Please, contact your administrator ... if you have one!")
369 print(error_message); stop(error_message)
370 }
371
372 file.symlink(singlefile_galaxyPath,singlefile_sampleName)
373 }
374 directory = "."
375
376 }
377 if(!is.null(zipfile) && (zipfile!="")) {
378 if(!file.exists(zipfile)){
379 error_message=paste("Cannot access the Zip file:",zipfile,". Please, contact your administrator ... if you have one!")
380 print(error_message)
381 stop(error_message)
382 }
383
384 #list all file in the zip file
385 #zip_files=unzip(zipfile,list=T)[,"Name"]
386
387 #unzip
388 suppressWarnings(unzip(zipfile, unzip="unzip"))
389
390 #get the directory name
391 filesInZip=unzip(zipfile, list=T);
392 directories=unique(unlist(lapply(strsplit(filesInZip$Name,"/"), function(x) x[1])));
393 directories=directories[!(directories %in% c("__MACOSX")) & file.info(directories)$isdir]
394 directory = "."
395 if (length(directories) == 1) directory = directories
396
397 cat("files_root_directory\t",directory,"\n")
398
399 }
400 }