Mercurial > repos > mmonsoor > probmetab

comparison lib.r @ 0:e13ec2c3fabe draft

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
planemo upload commit 25fd6a739741295e3f434e0be0286dee61e06825
author mmonsoor
date Mon, 04 Jul 2016 04:29:25 -0400
parents
children abcfa1648b66
comparison
equal deleted inserted replaced
-1:000000000000 0:e13ec2c3fabe
1 # lib.r ProbMetab version="1.0.0"
2 # Author: Misharl Monsoor ABIMS TEAM mmonsoor@sb-roscoff.fr
3 # Contributors: Yann Guitton and Jean-francois Martin
4
5
6 ##Main probmetab function launch by the Galaxy ProbMetab wrapper
7 probmetab = function(xa, xaP, xaN, variableMetadata, variableMetadataP, variableMetadataN, listArguments){
8 ##ONE MODE ACQUISITION##
9 if(listArguments[["mode_acquisition"]]=="one") {
10 comb=NULL
11
12 #Get the polarity from xa object
13 polarity=xa@polarity
14 #SNR option
15 if ("xsetnofill" %in% names(listArguments)) {
16 load(listArguments[["xsetnofill"]])
17 xsetnofill=xset
18 }
19 else{
20 xsetnofill=NULL
21 }
22 #Exclude samples
23 if ("toexclude" %in% names(listArguments)) {
24 toexclude=listArguments[["toexclude"]]
25 }
26 else {
27 toexclude=NULL
28 }
29 ionAnnot=get.annot(xa, polarity=polarity, allowMiss=listArguments[["allowMiss"]],xset=xsetnofill,toexclude=toexclude)
30 comb=NULL
31 }
32
33 ##TWO MODES ACQUISITION##
34 #Mode annotatediffreport
35 else if(listArguments[["inputs_mode"]]=="two"){
36 ##Prepare the objects that will be used for the get.annot function
37 comb=1
38
39
40 xsetPnofill=NULL
41 xsetNnofill=NULL
42 # TODO: a reactiver
43 #if ("xsetPnofill" %in% names(listArguments)) {
44 # load(listArguments[["xsetPnofill"]])
45 # xsetPnofill=xset
46 #}
47 #if ("xsetNnofill" %in% names(listArguments)) {
48 # load(listArguments[["xsetNnofill"]])
49 # xsetNnofill=xset
50 #}
51 # include CAMERA non-annotated compounds, and snr retrieval
52 # comb 2+ - used on Table 1
53 ionAnnotP2plus = get.annot(axP, allowMiss=listArguments[["allowMiss"]], xset=xsetPnofill,toexclude=listArguments[["toexclude"]])
54 ionAnnotN2plus = get.annot(axN, polarity="negative", allowMiss=listArguments[["allowMiss"]], xset=xsetNnofill,toexclude=listArguments[["toexclude"]])
55 ionAnnot = combineMolIon(ionAnnotP2plus, ionAnnotN2plus)
56 print(sum(ionAnnot$molIon[,3]==1))
57 print(sum(ionAnnot$molIon[,3]==0))
58 write.table(ionAnnot[1], sep="\t", quote=FALSE, row.names=FALSE, file="CombineMolIon.tsv")
59 #Merge variableMetadata Negative and positive acquisitions mode
60
61
62 #Mode combinexsannos TODO bug avec tableau issus de combinexsannos
63 #else {
64 #load(listArguments[["image_combinexsannos"]])
65 #image_combinexsannos=cAnnot
66 ##Prepare the objects that will be used for the combineMolIon function
67 #load(listArguments[["image_pos"]])
68 #image_pos=xa
69 #ionAnnot=combineMolIon(peaklist=cAnnot, cameraobj=image_pos, polarity="pos")
70 #}
71
72 }
73
74 ##DATABASE MATCHING##
75 if (listArguments[["kegg_db"]]=="KEGG"){
76 DB=build.database.kegg(orgID = NULL)
77 }
78 else{
79 table_list <<- NULL
80 ids=strsplit(listArguments[["kegg_db"]],",")
81 ids=ids[[1]]
82 if(length(ids)>1){
83 for(i in 1:length(ids)){
84 table_list[[i]] <- build.database.kegg(ids[i])
85 }
86 db_table=do.call("rbind",table_list)
87 DB=unique(db_table)
88 }
89 else{
90 DB=build.database.kegg(listArguments[["kegg_db"]])
91 }
92 }
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"]])
95 ##PROBABILITY RANKING##
96 # number of masses with candidates inside the fixed mass window
97 # and masses with more than one candidate
98 length(unique(reactionM[reactionM[,"id"]!="unknown",1]))
99 sum(table(reactionM[reactionM[,"id"]!="unknown",1])>1)
100 #if (listArguments[["useIso"]]){
101 #BUG TODO
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
104 # skip this step and use the reactionM as input to weigthM function.
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
107 #wl < weightM(isoPatt,intervals=seq(0,1000,by=500), offset=c(3.115712, 3.434146, 2.350798))
108
109 #isoPatt=incorporate.isotopes(ionAnnot, reactionM,comb=comb,var=listArguments[["var"]],DB=DB)
110
111 #wl = weightM(reactionM, useIso=true)
112 #}
113 #else {
114 #wl = weightM(reactionM, useIso=FALSE)
115 #}
116 wl =weightM(reactionM, useIso=FALSE)
117 w = design.connection(reactionM)
118 # Probability calculations
119 x = 1:ncol(wl$wm)
120 y = 1:nrow(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"]])
123 if(listArguments[["html"]]){
124 #Zip the EICS plot
125 system(paste('zip -r "Analysis_Report.zip" "AnalysisExample_fig"'))
126 }
127
128 # calculate the correlations and partial correlations and cross reference then with reactions
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
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
133 #two masses.
134 #It generates a list of estimated correlations and reactions.
135 corList=reac2cor(mw,ansConn$classTable,listArguments[["opt"]],listArguments[["corths"]],listArguments[["corprob"]],listArguments[["pcorprob"]])
136 ans=list("ansConn"=ansConn,"corList"=corList)
137 #Generate the siff table for CytoScape
138 cytoscape_output(corList,ansConn)
139
140
141 #Execute the merge_probmetab function to merge the variableMetadata table and annotations from ProbMetab results
142
143 if(listArguments[["mode_acquisition"]]=="one") {
144 #Retrocompatibility with previous annotateDiffreport variableMetadata dataframe (must replace mzmed column by mz, and rtmed by rt)
145 names(variableMetadata)[names(variableMetadata)=="mzmed"] <- "mz"
146 names(variableMetadata)[names(variableMetadata)=="rtmed"] <- "rt"
147 variableM=merge_probmetab(variableMetadata, ansConn)
148 write.table(variableM, sep="\t", quote=FALSE, row.names=FALSE, file="variableMetadata.tsv")
149 } else if (listArguments[["mode_acquisition"]]=="two") {
150 #Retrocompatibility with previous annotateDiffreport variableMetadata dataframe (must replace mzmed column by mz, and rtmed by rt)
151 names(variableMetadataP)[names(variableMetadata)=="mzmed"] <- "mz"
152 names(variableMetadataP)[names(variableMetadata)=="rtmed"] <- "rt"
153 names(variableMetadataN)[names(variableMetadata)=="mzmed"] <- "mz"
154 names(variableMetadataN)[names(variableMetadata)=="rtmed"] <- "rt"
155 variableMP=merge_probmetab(variableMetadataP, ansConn)
156 write.table(variableMP, sep="\t", quote=FALSE, row.names=FALSE, file="variableMetadata_Positive.tsv")
157 variableMN=merge_probmetab(variableMetadataN, ansConn)
158 write.table(variableMN, sep="\t", quote=FALSE, row.names=FALSE, file="variableMetadata_Negative.tsv")
159 }
160
161 return(ans)
162
163 }
164
165 ##Function that generates a siff table for CytoScape
166 cytoscape_output=function(corList,ansConn){
167 signif_cor=as.data.frame(corList$signif.cor)
168 classTable=as.data.frame(ansConn$classTable)
169 #Siff table
170 siff_table=cbind(signif_cor["node1"],signif_cor["cor"],signif_cor["node2"])
171 #attribute table output for Cytoscape
172
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] == ""){
175 classTable[i, c(1, 4, 6, 7)] <- classTable[i - 1, c(1, 4, 6, 7)]
176 }
177 msel <- as.matrix(classTable[, 1:7])
178 msel <- cbind(msel[, 6], msel[,-6])
179 colnames(msel)[1] <- "Id"
180 msel[, 1] <- sub("^\\s+", "", msel[, 1])
181 colnames(msel)[1] <- "Id"
182 ids <- unique(msel[, 1])
183 attrMatrix <- matrix("", nrow = length(ids), ncol = ncol(msel)-1)
184 for (i in 1:length(ids)) {
185 attrMatrix[i, 1] <- unique(msel[msel[, 1] == ids[i],
186 2])
187 attrMatrix[i, 2] <- paste("[", paste(msel[msel[,
188 1] == ids[i], 3], collapse = ", "), "]", sep = "")
189 attrMatrix[i, 3] <- paste("[", paste(msel[msel[,
190 1] == ids[i], 4], collapse = ", "), "]", sep = "")
191 attrMatrix[i, 4] <- unique(msel[msel[, 1] == ids[i],
192 5])
193 attrMatrix[i, 5] <- paste("[", paste(msel[msel[,
194 1] == ids[i], 6], collapse = ", "), "]", sep = "")
195 attrMatrix[i, 6] <- unique(msel[msel[, 1] == ids[i],
196 7])
197 }
198 ids <- as.numeric(unique(msel[, 1]))
199 attrMatrix <- cbind(ids, attrMatrix)
200 colnames(attrMatrix) <- colnames(msel)
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")
203 write.table(siff_table, sep="\t", quote=FALSE, row.names=FALSE, file="sif.tsv")
204
205 return(attrMatrix)
206 }
207
208 ##Functions written by Jean-Francois Martin
209
210 deter_ioni <- function (aninfo, pm)
211 {
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
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
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
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) {
220 if (as.numeric(aninfo[1]) - pm <0) {esi <- "n"} else {esi <- "p"}
221 } else
222 if (!is.na(aninfo[4])) {
223 anstr <- aninfo[4]
224 # cat(anstr)
225 if ((grepl("]+",anstr,fixed=T)==T) || (grepl("]2+",anstr,fixed=T)==T) || (grepl("]3+",anstr,fixed=T)==T)) { esi <- "p"}
226 else
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")
229 } else
230 { esi <- "u"}
231
232 return(esi)
233 }
234
235
236 merge_probmetab <- function(metaVar,ansConn) {
237 ## Parse ProbMetab information result file and merge in variable_metaData initial file
238 ## inputs :
239 ## metaVar : data.frame of metadataVariable input of probmetab function
240 ## ansConn : data.frame of ProbMetab result
241 ## output : dataframe with Probmetab results merge with variableMetadata
242 ## Constante
243 ## iannot : indice de la colonne annotation dans le resultat de probMetab
244 iannot <- 4
245
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
248 ions <- paste ("M",round(metaVar$mz,3),"T",round(metaVar$rt,1),sep="")
249 metaVar <- data.frame(ions,metaVar)
250
251 ###### Result data.frame from ProbMetab result list
252 an_ini <- ansConn$classTable
253
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
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)
258 mz <- rep(0,dim(an)[1])
259 rt <- rep(0,dim(an)[1])
260 propmz <- rep(0,dim(an)[1])
261 ioni <- rep("u",dim(an)[1])
262
263 ## parse the column annotation and define ionisation mode
264 for (i in 1:dim(an)[1]) {
265 if (an[i,1] != "") {
266 info_mzrt <- unlist(strsplit(an[i,iannot],"#"))
267 propmz[i] <- as.numeric(an[i,1])
268 mz[i] <- as.numeric(info_mzrt[1])
269 rt[i] <- as.numeric(info_mzrt[2])
270 ioni[i] <- deter_ioni(info_mzrt,as.numeric(an[i,1]))
271 }
272 else {
273 propmz[i] <- as.numeric(propmz[i-1])
274 mz[i] <- as.numeric(mz[i-1])
275 rt[i] <- as.numeric(rt[i-1])
276 ioni[i] <- ioni[i-1]
277 }
278 }
279
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.
282 ions <- paste ("M",round(mz,3),"T",round(rt,1),sep="")
283 an <- data.frame(ions,ioni,propmz,mz,rt,an)
284
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 ";"
287 li <- as.matrix(table(an$propmz))
288 li <- data.frame(dimnames(li)[1],li)
289 dimnames(li)[[2]][1] <- "propmz"
290 ions <- rep("u",dim(li)[1])
291 propmz <- rep(0,dim(li)[1])
292 mpc <- rep("c",dim(li)[1])
293 proba <- rep("p",dim(li)[1])
294 c <- 0
295 while (c < dim(li)[1]) {
296 c <- c + 1
297 suban <- an[an$propmz==li[c,1],]
298 ions[c] <- as.character(suban[1,1])
299 propmz[c] <- as.numeric(suban[1,3])
300 mpc[c] <- paste(suban[,7],collapse=";")
301 proba[c] <- paste(as.character(suban[,8]),collapse=";")
302 }
303
304 ## Creation of the data.frame with 1 row per ions
305 anc <- data.frame(ions,propmz,mpc,proba)
306 anc <- anc[order(anc[,1]),]
307
308 metaVarFinal <- merge(metaVar, anc, by.x=1, by.y=1, all.x=T, all.y=T)
309 metaVarFinal <- metaVarFinal[,-1]
310 #write.table(metaVarFinal,file="res.txt", sep="\t", row.names=F, quote=F)
311
312 return (metaVarFinal)
313 }
314
315 # RETROCOMPATIBILITE avec ancienne version de annotate
316 getVariableMetadata = function(xa) {
317 # --- variableMetadata ---
318 peakList=getPeaklist(xa)
319 peakList=cbind(groupnames(xa@xcmsSet),peakList); colnames(peakList)[1] = c("name");
320 variableMetadata=peakList[,!(colnames(peakList) %in% c(sampnames(xa@xcmsSet)))]
321 variableMetadata$name= paste("M",round(variableMetadata$mz),"T",round(variableMetadata$rt),sep="")
322 return (variableMetadata)
323 }