--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/lib.r	Mon Jul 04 04:29:25 2016 -0400
@@ -0,0 +1,323 @@
+# lib.r ProbMetab version="1.0.0"
+# Author: Misharl Monsoor ABIMS TEAM mmonsoor@sb-roscoff.fr
+# Contributors: Yann Guitton and Jean-francois Martin
+
+
+##Main probmetab function launch by the Galaxy ProbMetab wrapper
+probmetab = function(xa, xaP, xaN, variableMetadata, variableMetadataP, variableMetadataN, listArguments){
+	##ONE MODE ACQUISITION##
+	if(listArguments[["mode_acquisition"]]=="one") {
+		comb=NULL
+
+		#Get the polarity from xa object
+		polarity=xa@polarity
+		#SNR option
+		if ("xsetnofill" %in% names(listArguments)) {
+			load(listArguments[["xsetnofill"]])
+			xsetnofill=xset
+		}
+		else{
+			xsetnofill=NULL
+		}
+		#Exclude samples
+		if ("toexclude" %in% names(listArguments)) {
+			toexclude=listArguments[["toexclude"]]
+		}
+		else {
+			toexclude=NULL
+		}
+		ionAnnot=get.annot(xa, polarity=polarity, allowMiss=listArguments[["allowMiss"]],xset=xsetnofill,toexclude=toexclude)
+		comb=NULL
+	}
+
+	##TWO MODES ACQUISITION##
+	#Mode annotatediffreport
+	else if(listArguments[["inputs_mode"]]=="two"){
+		##Prepare the objects that will be used for the get.annot function
+		comb=1
+		
+
+		xsetPnofill=NULL
+		xsetNnofill=NULL
+		# TODO: a reactiver		
+		#if ("xsetPnofill" %in% names(listArguments)) {
+		#	load(listArguments[["xsetPnofill"]])
+		#	xsetPnofill=xset
+		#}
+		#if ("xsetNnofill" %in% names(listArguments)) {
+		#	load(listArguments[["xsetNnofill"]])
+		#	xsetNnofill=xset
+		#}
+		# include CAMERA non-annotated compounds, and snr retrieval 
+		# comb 2+ - used on Table 1	
+		ionAnnotP2plus = get.annot(axP, allowMiss=listArguments[["allowMiss"]], xset=xsetPnofill,toexclude=listArguments[["toexclude"]]) 
+		ionAnnotN2plus = get.annot(axN, polarity="negative", allowMiss=listArguments[["allowMiss"]], xset=xsetNnofill,toexclude=listArguments[["toexclude"]])
+		ionAnnot = combineMolIon(ionAnnotP2plus, ionAnnotN2plus)
+		print(sum(ionAnnot$molIon[,3]==1))
+		print(sum(ionAnnot$molIon[,3]==0))
+		write.table(ionAnnot[1], sep="\t", quote=FALSE, row.names=FALSE, file="CombineMolIon.tsv")
+		#Merge variableMetadata Negative and positive acquisitions mode
+		
+
+		#Mode combinexsannos TODO bug avec tableau issus de combinexsannos
+		#else {
+			#load(listArguments[["image_combinexsannos"]])
+			#image_combinexsannos=cAnnot
+			##Prepare the objects that will be used for the combineMolIon function
+			#load(listArguments[["image_pos"]])
+			#image_pos=xa
+			#ionAnnot=combineMolIon(peaklist=cAnnot, cameraobj=image_pos, polarity="pos")
+		#}
+		
+	}
+
+	##DATABASE MATCHING##
+	if (listArguments[["kegg_db"]]=="KEGG"){		
+		DB=build.database.kegg(orgID = NULL)
+	}
+	else{	
+		table_list <<- NULL
+		ids=strsplit(listArguments[["kegg_db"]],",")
+		ids=ids[[1]]
+		if(length(ids)>1){
+			for(i in 1:length(ids)){
+				 table_list[[i]] <- build.database.kegg(ids[i])
+			}
+			db_table=do.call("rbind",table_list)
+			DB=unique(db_table)
+		}
+		else{
+			DB=build.database.kegg(listArguments[["kegg_db"]])
+		}
+	}	
+	#Matching des mass exactes mesurees avec les masses des compounds KEGG (pas M+H ou M-H)
+	reactionM = create.reactionM(DB, ionAnnot, ppm.tol=listArguments[["ppm_tol"]])
+	##PROBABILITY RANKING##
+	# number of masses with candidates inside the fixed mass window
+	# and masses with more than one candidate
+	length(unique(reactionM[reactionM[,"id"]!="unknown",1])) 
+	sum(table(reactionM[reactionM[,"id"]!="unknown",1])>1)
+	#if (listArguments[["useIso"]]){
+		#BUG TODO
+		# Calculate the ratio between observed and theoretical isotopic patterns.
+		# If you don't have an assessment of carbon offset to carbon number prediction 
+		# skip this step and use the reactionM as input to weigthM function. 
+		#isoPatt < incorporate.isotopes(comb2plus, reactionM, , samp=12:23, DB=DB)  
+		#  calculate   the   likelihood   of   each   mass   to   compound   assignment   using   mass   accuracy,and isotopic pattern, when present
+		#wl < weightM(isoPatt,intervals=seq(0,1000,by=500), offset=c(3.115712, 3.434146, 2.350798))
+		
+			#isoPatt=incorporate.isotopes(ionAnnot, reactionM,comb=comb,var=listArguments[["var"]],DB=DB)
+
+		#wl = weightM(reactionM, useIso=true)
+	#}
+	#else {
+		#wl = weightM(reactionM, useIso=FALSE)
+	#}
+	wl =weightM(reactionM, useIso=FALSE)	
+	w = design.connection(reactionM)
+	# Probability calculations
+	x = 1:ncol(wl$wm)
+	y = 1:nrow(wl$wm)
+	conn = gibbs.samp(x, y, 5000, w, wl$wm)
+	ansConn = export.class.table(conn, reactionM, ionAnnot, DB=DB,html=listArguments[["html"]],filename="AnalysisExample",prob=listArguments[["prob"]])
+	if(listArguments[["html"]]){
+		#Zip the EICS plot
+		system(paste('zip -r "Analysis_Report.zip" "AnalysisExample_fig"'))
+	}
+	
+	# calculate the correlations and partial correlations and cross reference then with reactions
+	mw=which(w==1,arr.ind=TRUE)
+	#reac2cor function : Use the intensity of putative molecules in repeated samples to calculate correlations and partial
+	#correlation in a user defined threshold of false discovery rate for significance testing. After the
+	#correlation test the function also overlay significant correlations with all putative reactions between
+	#two masses.
+	#It generates a list of estimated correlations and reactions.
+	corList=reac2cor(mw,ansConn$classTable,listArguments[["opt"]],listArguments[["corths"]],listArguments[["corprob"]],listArguments[["pcorprob"]])
+	ans=list("ansConn"=ansConn,"corList"=corList)
+	#Generate the siff table for CytoScape
+	cytoscape_output(corList,ansConn)
+
+
+	#Execute the merge_probmetab function to merge the variableMetadata table and annotations from ProbMetab results
+	
+	if(listArguments[["mode_acquisition"]]=="one") {
+		#Retrocompatibility with previous annotateDiffreport variableMetadata dataframe (must replace mzmed column by mz, and rtmed by rt)
+		names(variableMetadata)[names(variableMetadata)=="mzmed"] <- "mz"
+		names(variableMetadata)[names(variableMetadata)=="rtmed"] <- "rt"
+		variableM=merge_probmetab(variableMetadata, ansConn)
+		write.table(variableM, sep="\t", quote=FALSE, row.names=FALSE, file="variableMetadata.tsv")
+	} else if (listArguments[["mode_acquisition"]]=="two") {
+		#Retrocompatibility with previous annotateDiffreport variableMetadata dataframe (must replace mzmed column by mz, and rtmed by rt)
+		names(variableMetadataP)[names(variableMetadata)=="mzmed"] <- "mz"
+		names(variableMetadataP)[names(variableMetadata)=="rtmed"] <- "rt"
+		names(variableMetadataN)[names(variableMetadata)=="mzmed"] <- "mz"
+		names(variableMetadataN)[names(variableMetadata)=="rtmed"] <- "rt"
+		variableMP=merge_probmetab(variableMetadataP, ansConn)
+		write.table(variableMP, sep="\t", quote=FALSE, row.names=FALSE, file="variableMetadata_Positive.tsv")
+		variableMN=merge_probmetab(variableMetadataN, ansConn)
+		write.table(variableMN, sep="\t", quote=FALSE, row.names=FALSE, file="variableMetadata_Negative.tsv")
+	}
+
+	return(ans)
+
+}
+
+##Function that generates a siff table for CytoScape
+cytoscape_output=function(corList,ansConn){
+	signif_cor=as.data.frame(corList$signif.cor)
+	classTable=as.data.frame(ansConn$classTable)
+	#Siff table
+	siff_table=cbind(signif_cor["node1"],signif_cor["cor"],signif_cor["node2"])
+	#attribute table output for Cytoscape
+
+	## START  Code part from the export2cytoscape function of ProbMetab written by Ricardo R. Silva 
+	for (i in 1:nrow(classTable)) if (classTable[i, 1] == ""){
+		classTable[i, c(1, 4, 6, 7)] <- classTable[i - 1, c(1, 4, 6, 7)]
+	}
+ 	msel <- as.matrix(classTable[, 1:7])
+ 	msel <- cbind(msel[, 6], msel[,-6])
+ 	colnames(msel)[1] <- "Id"
+ 	msel[, 1] <- sub("^\\s+", "", msel[, 1])
+ 	colnames(msel)[1] <- "Id"
+	ids <- unique(msel[, 1])
+ 	attrMatrix <- matrix("", nrow = length(ids), ncol = ncol(msel)-1)
+	for (i in 1:length(ids)) {
+		    attrMatrix[i, 1] <- unique(msel[msel[, 1] == ids[i], 
+		        2])
+		    attrMatrix[i, 2] <- paste("[", paste(msel[msel[, 
+		        1] == ids[i], 3], collapse = ", "), "]", sep = "")
+		    attrMatrix[i, 3] <- paste("[", paste(msel[msel[, 
+		        1] == ids[i], 4], collapse = ", "), "]", sep = "")
+		    attrMatrix[i, 4] <- unique(msel[msel[, 1] == ids[i], 
+		        5])
+		    attrMatrix[i, 5] <- paste("[", paste(msel[msel[, 
+		        1] == ids[i], 6], collapse = ", "), "]", sep = "")
+		    attrMatrix[i, 6] <- unique(msel[msel[, 1] == ids[i], 
+		        7])
+        }
+	ids <- as.numeric(unique(msel[, 1]))
+	attrMatrix <- cbind(ids, attrMatrix)
+	colnames(attrMatrix) <- colnames(msel)
+	## END Code part from the export2cytoscape function of ProbMetab writieen by Ricardo R. Silva 
+	write.table(attrMatrix, sep="\t", quote=FALSE, row.names=FALSE, file="Analysis_Report.tsv")
+	write.table(siff_table, sep="\t", quote=FALSE, row.names=FALSE, file="sif.tsv")
+
+	return(attrMatrix)
+}
+
+##Functions written by Jean-Francois Martin
+
+deter_ioni <- function (aninfo, pm)
+{
+  # determine ionisation in ProbMetab result file, used in function merge_probmetab
+  # input : for 1 ion, aninfo = string with m/z rt and CAMERA annotation from ProbMetab result file
+  # if the difference between m/z and the probmetab proposed mass is ~1 we use the sign (positive or negative) of this diference 
+  # to define the type of ionisation
+  # If adduct or fragment was detected, therefore diff >>1 and so, we search for substring "+" ou "2+" ou "3+" ou "-"... 
+  # to define the type of ionisation
+  # aninfo : vecteur of character resulting of the parsing(sep="#") of the probmetab annotation
+  if (round(abs(as.numeric(aninfo[1]) - pm),0) ==1) {
+    if (as.numeric(aninfo[1]) - pm <0) {esi <- "n"} else {esi <- "p"}
+  } else 
+    if (!is.na(aninfo[4])) {
+      anstr <- aninfo[4]
+      # cat(anstr)
+      if ((grepl("]+",anstr,fixed=T)==T) || (grepl("]2+",anstr,fixed=T)==T) || (grepl("]3+",anstr,fixed=T)==T)) { esi <- "p"}
+      else 
+        if ((grepl("]-",anstr,fixed=T)==T) || (grepl("]2-",anstr,fixed=T)==T) || (grepl("]3-",anstr,fixed=T)==T)) { esi <- "n"}
+      # cat(" ioni ",esi,"\n")
+    } else
+    { esi <- "u"} 
+  
+  return(esi)
+}
+
+
+merge_probmetab <- function(metaVar,ansConn) {
+  ## Parse ProbMetab information result file and merge in variable_metaData initial file
+  ##  inputs : 
+  ##      metaVar : data.frame of metadataVariable input of probmetab function
+  ##     ansConn  : data.frame of ProbMetab result
+  ## output : dataframe with Probmetab results merge with variableMetadata
+  ## Constante
+  ## iannot : indice de la colonne annotation dans le resultat de probMetab
+  iannot <- 4
+  
+  ## definition of an unique identification of ions mz with 3 decimals and rt(sec) with 1 decimal to avoid 
+  ## duplicate ions name in the diffreport result file
+  ions <- paste ("M",round(metaVar$mz,3),"T",round(metaVar$rt,1),sep="")
+  metaVar <- data.frame(ions,metaVar)
+  
+  ###### Result data.frame from ProbMetab result list
+  an_ini <- ansConn$classTable
+  
+  ## Suppression of rows without  mz and rt or unknown and columns of intensities
+  ## COLUMNS SUBSCRIPTS HAVE TO BE CHECKED WITh DIFFERENT RESULTS FILES 
+  an <- an_ini[(an_ini[,2]!="unknown"),c(1,2,3,7)]
+  ## initialisation of vectors receiving the result of the parse of the column annotation (subscrip iannot) 
+  mz <- rep(0,dim(an)[1])
+  rt <- rep(0,dim(an)[1])
+  propmz <- rep(0,dim(an)[1])
+  ioni <- rep("u",dim(an)[1])
+  
+  ## parse the column annotation and define ionisation mode
+  for (i in 1:dim(an)[1]) {
+    if (an[i,1] != "") {
+      info_mzrt <- unlist(strsplit(an[i,iannot],"#"))
+      propmz[i] <- as.numeric(an[i,1])
+      mz[i] <- as.numeric(info_mzrt[1])
+      rt[i] <- as.numeric(info_mzrt[2])
+      ioni[i] <- deter_ioni(info_mzrt,as.numeric(an[i,1]))
+    }
+    else {
+      propmz[i] <- as.numeric(propmz[i-1])
+      mz[i] <- as.numeric(mz[i-1])
+      rt[i] <- as.numeric(rt[i-1])
+      ioni[i] <- ioni[i-1]
+    }
+  }
+  
+  ## definition of an unique identification of ions : mz with 3 decimals and rt(sec) with 1 decimal
+  ## The same as for the metadataVariable data.frame to match with.
+  ions <- paste ("M",round(mz,3),"T",round(rt,1),sep="")
+  an <- data.frame(ions,ioni,propmz,mz,rt,an)
+  
+  ## transposition of the different probmetab annotations which are in different rows in the initial result data.frame
+  ## on only 1 row separated with a ";"
+  li <- as.matrix(table(an$propmz))
+  li <- data.frame(dimnames(li)[1],li)
+  dimnames(li)[[2]][1] <- "propmz"
+  ions   <- rep("u",dim(li)[1])
+  propmz <- rep(0,dim(li)[1])
+  mpc    <- rep("c",dim(li)[1])
+  proba  <- rep("p",dim(li)[1])
+  c <- 0
+  while (c < dim(li)[1]) {
+    c <- c + 1
+    suban     <- an[an$propmz==li[c,1],]
+    ions[c]   <- as.character(suban[1,1])
+    propmz[c] <- as.numeric(suban[1,3])
+    mpc[c]    <- paste(suban[,7],collapse=";")
+    proba[c]  <- paste(as.character(suban[,8]),collapse=";") 
+  }
+  
+  ## Creation of the data.frame with 1 row per ions
+  anc <- data.frame(ions,propmz,mpc,proba)
+  anc <- anc[order(anc[,1]),]
+  
+  metaVarFinal <- merge(metaVar, anc, by.x=1, by.y=1, all.x=T, all.y=T)
+  metaVarFinal <- metaVarFinal[,-1]
+  #write.table(metaVarFinal,file="res.txt", sep="\t", row.names=F, quote=F)
+  
+  return (metaVarFinal)
+}
+
+# RETROCOMPATIBILITE avec ancienne version de annotate
+getVariableMetadata = function(xa) {
+	# --- variableMetadata ---
+	peakList=getPeaklist(xa)
+	peakList=cbind(groupnames(xa@xcmsSet),peakList); colnames(peakList)[1] = c("name");
+	variableMetadata=peakList[,!(colnames(peakList) %in% c(sampnames(xa@xcmsSet)))]
+	variableMetadata$name= paste("M",round(variableMetadata$mz),"T",round(variableMetadata$rt),sep="")
+	return (variableMetadata)
+}