diff lib.r @ 9:972152941964 draft

planemo upload commit a634879c0e651eb0eb7b435a107ee40cf30524fa

--- a/lib.r	Wed Feb 01 12:23:24 2017 -0500
+++ b/lib.r	Fri Feb 10 11:11:49 2017 -0500
@@ -1,188 +1,225 @@
 # lib.r version="2.2.1"
 
+#@author G. Le Corguille
 #The function create a pdf from the different png generated by diffreport
 diffreport_png2pdf <- function(filebase, new_file_path) {
 
-  pdfEicOutput = paste(new_file_path,filebase,"-eic_visible_pdf",sep="")
-  pdfBoxOutput = paste(new_file_path,filebase,"-box_visible_pdf",sep="")
+    pdfEicOutput = paste(new_file_path,filebase,"-eic_visible_pdf",sep="")
+    pdfBoxOutput = paste(new_file_path,filebase,"-box_visible_pdf",sep="")
+
+    system(paste("gm convert ",filebase,"_eic/*.png ",filebase,"_eic.pdf",sep=""))
+    system(paste("gm convert ",filebase,"_box/*.png ",filebase,"_box.pdf",sep=""))
+
+    file.copy(paste(filebase,"_eic.pdf",sep=""), pdfEicOutput)
+    file.copy(paste(filebase,"_box.pdf",sep=""), pdfBoxOutput)
+}
 
-  system(paste("gm convert ",filebase,"_eic/*.png ",filebase,"_eic.pdf",sep=""))
-  system(paste("gm convert ",filebase,"_box/*.png ",filebase,"_box.pdf",sep=""))
+#@author G. Le Corguille
+#This function convert if it is required the Retention Time in minutes
+RTSecondToMinute <- function(variableMetadata, convertRTMinute) {
+    if (convertRTMinute){
+        #converting the retention times (seconds) into minutes
+        print("converting the retention times into minutes in the variableMetadata")
+        variableMetadata[,"rt"]=variableMetadata[,"rt"]/60
+        variableMetadata[,"rtmin"]=variableMetadata[,"rtmin"]/60
+        variableMetadata[,"rtmax"]=variableMetadata[,"rtmax"]/60
+    }
+    return (variableMetadata)
+}
 
-  file.copy(paste(filebase,"_eic.pdf",sep=""), pdfEicOutput)
-  file.copy(paste(filebase,"_box.pdf",sep=""), pdfBoxOutput)
+#@author G. Le Corguille
+#This function format ions identifiers
+formatIonIdentifiers <- function(variableMetadata, numDigitsRT=0, numDigitsMZ=0) {
+    splitDeco = strsplit(as.character(variableMetadata$name),"_")
+    idsDeco = sapply(splitDeco, function(x) { deco=unlist(x)[2]; if (is.na(deco)) return ("") else return(paste0("_",deco)) })
+    namecustom = make.unique(paste0("M",round(variableMetadata[,"mz"],numDigitsMZ),"T",round(variableMetadata[,"rt"],numDigitsRT),idsDeco))
+    variableMetadata=cbind(name=variableMetadata$name, namecustom=namecustom, variableMetadata[,!(colnames(variableMetadata) %in% c("name"))])
+    return(variableMetadata)
 }
 
 #The function annotateDiffreport without the corr function which bugs
 annotatediff <- function(xset=xset, listArguments=listArguments, variableMetadataOutput="variableMetadata.tsv", dataMatrixOutput="dataMatrix.tsv",new_file_path=NULL) {
-  # Resolve the bug with x11, with the function png
-  options(bitmapType='cairo')
+    # Resolve the bug with x11, with the function png
+    options(bitmapType='cairo')
+
+    #Check if the fillpeaks step has been done previously, if it hasn't, there is an error message and the execution is stopped.
+    res=try(is.null(xset@filled))
 
-  #Check if the fillpeaks step has been done previously, if it hasn't, there is an error message and the execution is stopped.
-  res=try(is.null(xset@filled))
+    # ------ annot -------
+    listArguments[["calcCiS"]]=as.logical(listArguments[["calcCiS"]])
+    listArguments[["calcIso"]]=as.logical(listArguments[["calcIso"]])
+    listArguments[["calcCaS"]]=as.logical(listArguments[["calcCaS"]])
 
-  # ------ annot -------
-  listArguments[["calcCiS"]]=as.logical(listArguments[["calcCiS"]])
-  listArguments[["calcIso"]]=as.logical(listArguments[["calcIso"]])
-  listArguments[["calcCaS"]]=as.logical(listArguments[["calcCaS"]])
+    # common parameters
+    listArguments4annotate = list(object=xset,
+    nSlaves=listArguments[["nSlaves"]],sigma=listArguments[["sigma"]],perfwhm=listArguments[["perfwhm"]],
+    maxcharge=listArguments[["maxcharge"]],maxiso=listArguments[["maxiso"]],minfrac=listArguments[["minfrac"]],
+    ppm=listArguments[["ppm"]],mzabs=listArguments[["mzabs"]],quick=listArguments[["quick"]],
+    polarity=listArguments[["polarity"]],max_peaks=listArguments[["max_peaks"]],intval=listArguments[["intval"]])
 
-  #graphMethod parameter bugs where this parameter is not defined in quick=true
-  if(listArguments[["quick"]]==TRUE) {
-    xa= annotate(object=xset,nSlaves=1,sigma=listArguments[["sigma"]],perfwhm=listArguments[["perfwhm"]],maxcharge=listArguments[["maxcharge"]],maxiso=listArguments[["maxiso"]],minfrac=listArguments[["minfrac"]],ppm=listArguments[["ppm"]],mzabs=listArguments[["mzabs"]],quick=listArguments[["quick"]],polarity=listArguments[["polarity"]],max_peaks=listArguments[["max_peaks"]],intval=listArguments[["intval"]])
-  }
-  else {
-    xa= annotate(object=xset,nSlaves=1,sigma=listArguments[["sigma"]],perfwhm=listArguments[["perfwhm"]],graphMethod=listArguments[["graphMethod"]],cor_eic_th=listArguments[["cor_eic_th"]],pval=listArguments[["pval"]],calcCiS=listArguments[["calcCiS"]],calcIso=listArguments[["calcIso"]],calcCaS=listArguments[["calcCaS"]],multiplier=listArguments[["multiplier"]],maxcharge=listArguments[["maxcharge"]],maxiso=listArguments[["maxiso"]],minfrac=listArguments[["minfrac"]],ppm=listArguments[["ppm"]],mzabs=listArguments[["mzabs"]],quick=listArguments[["quick"]],polarity=listArguments[["polarity"]],max_peaks=listArguments[["max_peaks"]],intval=listArguments[["intval"]])
+    # quick == FALSE
+    if(listArguments[["quick"]]==FALSE) {
+        listArguments4annotate = append(listArguments4annotate,
+            list(graphMethod=listArguments[["graphMethod"]],cor_eic_th=listArguments[["cor_eic_th"]],pval=listArguments[["pval"]],
+            calcCiS=listArguments[["calcCiS"]],calcIso=listArguments[["calcIso"]],calcCaS=listArguments[["calcCaS"]]))
+        # no ruleset
+        if (!is.null(listArguments[["multiplier"]])) {
+            listArguments4annotate = append(listArguments4annotate,
+                list(multiplier=listArguments[["multiplier"]]))
+        }
+        # ruleset
+        else {
+            rulset=read.table(listArguments[["rules"]], h=T, sep=";")
+            if (ncol(rulset) < 4) rulset=read.table(listArguments[["rules"]], h=T, sep="\t")
+            if (ncol(rulset) < 4) rulset=read.table(listArguments[["rules"]], h=T, sep=",")
+            if (ncol(rulset) < 4) {
+                error_message="Your ruleset file seems not well formatted. The column separators accepted are ; , and tabulation"
+                print(error_message)
+                stop(error_message)
+            }
 
-  }  
-  peakList=getPeaklist(xa,intval=listArguments[["intval"]])
-  peakList=cbind(groupnames(xa@xcmsSet),peakList); colnames(peakList)[1] = c("name");
+            listArguments4annotate = append(listArguments4annotate,
+                list(rules=rulset))
+        }
+    }
 
 
-  # --- Multi condition : diffreport --- 
-  diffrep=NULL
-  if (!is.null(listArguments[["runDiffreport"]]) & nlevels(sampclass(xset))>=2) {
+    # launch annotate
+    xa = do.call("annotate", listArguments4annotate)
+    peakList=getPeaklist(xa,intval=listArguments[["intval"]])
+    peakList=cbind(groupnames(xa@xcmsSet),peakList); colnames(peakList)[1] = c("name");
+
+    # --- dataMatrix ---
+    dataMatrix = peakList[,(make.names(colnames(peakList)) %in% c("name", make.names(sampnames(xa@xcmsSet))))]
+    write.table(dataMatrix, sep="\t", quote=FALSE, row.names=FALSE, file=dataMatrixOutput)
+
+
+    # --- Multi condition : diffreport ---
+    diffrep=NULL
+    if (!is.null(listArguments[["runDiffreport"]]) & nlevels(sampclass(xset))>=2) {
     #Check if the fillpeaks step has been done previously, if it hasn't, there is an error message and the execution is stopped.
     res=try(is.null(xset@filled))
     classes=levels(sampclass(xset))
     x=1:(length(classes)-1)
     for (i in seq(along=x) ) {
-      y=1:(length(classes))
-      for (n in seq(along=y)){
-        if(i+n <= length(classes)){
-          filebase=paste(classes[i],class2=classes[i+n],sep="-vs-")
+        y=1:(length(classes))
+        for (n in seq(along=y)){
+            if(i+n <= length(classes)){
+                filebase=paste(classes[i],class2=classes[i+n],sep="-vs-")
+
+                diffrep=diffreport(object=xset,class1=classes[i],class2=classes[i+n],filebase=filebase,eicmax=listArguments[["eicmax"]],eicwidth=listArguments[["eicwidth"]],sortpval=TRUE,value=listArguments[["value"]],h=listArguments[["h"]],w=listArguments[["w"]],mzdec=listArguments[["mzdec"]])
 
-          diffrep=diffreport(object=xset,class1=classes[i],class2=classes[i+n],filebase=filebase,eicmax=listArguments[["eicmax"]],eicwidth=listArguments[["eicwidth"]],sortpval=TRUE,value=listArguments[["value"]],h=listArguments[["h"]],w=listArguments[["w"]],mzdec=listArguments[["mzdec"]])
-          #combines results
-          diffreportTSV=merge(peakList, diffrep[,c("name","fold","tstat","pvalue")], by.x="name", by.y="name", sort=F)
-          diffreportTSV=cbind(diffreportTSV[,!(colnames(diffreportTSV) %in% c(sampnames(xa@xcmsSet)))],diffreportTSV[,(colnames(diffreportTSV) %in% c(sampnames(xa@xcmsSet)))])
+                diffrepOri = diffrep
+
+                # renamming of the column rtmed to rt to fit with camera peaklist function output
+                colnames(diffrep)[colnames(diffrep)=="rtmed"] <- "rt"
+                colnames(diffrep)[colnames(diffrep)=="mzmed"] <- "mz"
+
+                # combines results and reorder columns
+                diffrep = merge(peakList, diffrep[,c("name","fold","tstat","pvalue")], by.x="name", by.y="name", sort=F)
+                diffrep = cbind(diffrep[,!(colnames(diffrep) %in% c(sampnames(xa@xcmsSet)))],diffrep[,(colnames(diffrep) %in% c(sampnames(xa@xcmsSet)))])
 
-          if(listArguments[["sortpval"]]){
-            diffreportTSV=diffreportTSV[order(diffreportTSV$pvalue), ]
-          }
-    
-          if (listArguments[["convert_param"]]){
-              #converting the retention times (seconds) into minutes
-              diffreportTSV$rt=diffreportTSV$rt/60;diffreportTSV$rtmin=diffreportTSV$rtmin/60; diffreportTSV$rtmax=diffreportTSV$rtmax/60;
-          }
-          write.table(diffreportTSV, sep="\t", quote=FALSE, row.names=FALSE, file=paste(new_file_path,filebase,"-tabular_visible_tabular",sep=""))
+                diffrep = RTSecondToMinute(diffrep, listArguments[["convertRTMinute"]])
+                diffrep = formatIonIdentifiers(diffrep, numDigitsRT=listArguments[["numDigitsRT"]], numDigitsMZ=listArguments[["numDigitsMZ"]])
+
+                if(listArguments[["sortpval"]]){
+                    diffrep=diffrep[order(diffrep$pvalue), ]
+                }
 
-          if (listArguments[["eicmax"]] != 0) {
-              diffreport_png2pdf(filebase, new_file_path)
-          }
+                write.table(diffrep, sep="\t", quote=FALSE, row.names=FALSE, file=paste(new_file_path,filebase,"-tabular_visible_tabular",sep=""))
+
+                if (listArguments[["eicmax"]] != 0) {
+                    diffreport_png2pdf(filebase, new_file_path)
+                }
+            }
         }
-      }
+        }
     }
-  }
-
-
 
 
-  # --- variableMetadata ---
-  variableMetadata=peakList[,!(make.names(colnames(peakList)) %in% c(make.names(sampnames(xa@xcmsSet))))]
-  # if we have 2 conditions, we keep stat of diffrep
-  if (!is.null(listArguments[["runDiffreport"]]) & nlevels(sampclass(xset))==2) {
-    variableMetadata = merge(variableMetadata, diffrep[,c("name","fold","tstat","pvalue")],by.x="name", by.y="name", sort=F)
-    if(exists("listArguments[[\"sortpval\"]]")){
-      variableMetadata=variableMetadata[order(variableMetadata$pvalue), ]
+    # --- variableMetadata ---
+    variableMetadata=peakList[,!(make.names(colnames(peakList)) %in% c(make.names(sampnames(xa@xcmsSet))))]
+    variableMetadata = RTSecondToMinute(variableMetadata, listArguments[["convertRTMinute"]])
+    variableMetadata = formatIonIdentifiers(variableMetadata, numDigitsRT=listArguments[["numDigitsRT"]], numDigitsMZ=listArguments[["numDigitsMZ"]])
+    # if we have 2 conditions, we keep stat of diffrep
+    if (!is.null(listArguments[["runDiffreport"]]) & nlevels(sampclass(xset))==2) {
+        variableMetadata = merge(variableMetadata, diffrep[,c("name","fold","tstat","pvalue")],by.x="name", by.y="name", sort=F)
+        if(exists("listArguments[[\"sortpval\"]]")){
+          variableMetadata=variableMetadata[order(variableMetadata$pvalue), ]
+        }
     }
-  }
 
-  variableMetadataOri=variableMetadata
-  if (listArguments[["convert_param"]]){
-    #converting the retention times (seconds) into minutes
-    print("converting the retention times into minutes in the variableMetadata")
-    variableMetadata$rt=variableMetadata$rt/60;variableMetadata$rtmin=variableMetadata$rtmin/60; variableMetadata$rtmax=variableMetadata$rtmax/60;
-  }
-  #Transform metabolites name
-  variableMetadata$name= paste("M",round(variableMetadata$mz,digits=listArguments[["num_digits"]]),"T",round(variableMetadata$rt),sep="")    
-  write.table(variableMetadata, sep="\t", quote=FALSE, row.names=FALSE, file=variableMetadataOutput)
+    variableMetadataOri=variableMetadata
+    write.table(variableMetadata, sep="\t", quote=FALSE, row.names=FALSE, file=variableMetadataOutput)
 
-  # --- dataMatrix ---
-  dataMatrix = peakList[,(make.names(colnames(peakList)) %in% c(make.names(sampnames(xa@xcmsSet))))]
-  dataMatrix=cbind(peakList$name,dataMatrix); colnames(dataMatrix)[1] = c("name");
-
-  if (listArguments[["convert_param"]]){
-    #converting the retention times (seconds) into minutes
-    print("converting the retention times into minutes in the dataMatrix ids")
-    peakList$rt=peakList$rt/60
-  }
-  dataMatrix$name= paste("M",round(peakList$mz,digits=listArguments[["num_digits"]]),"T",round(peakList$rt),sep="")
-  write.table(dataMatrix, sep="\t", quote=FALSE, row.names=FALSE, file=dataMatrixOutput)
-  
-  return(list("xa"=xa,"diffrep"=diffrep,"variableMetadata"=variableMetadataOri));
+    return(list("xa"=xa,"diffrep"=diffrepOri,"variableMetadata"=variableMetadataOri));
 
 }
 
 
-combinexsAnnos_function <- function(xaP, xaN, listOFlistArgumentsP,listOFlistArgumentsN, diffrepP=NULL,diffrepN=NULL,convert_param=FALSE,pos=TRUE,tol=2,ruleset=NULL,keep_meta=TRUE, variableMetadataOutput="variableMetadata.tsv"){
+combinexsAnnos_function <- function(xaP, xaN, listOFlistArgumentsP,listOFlistArgumentsN, diffrepP=NULL,diffrepN=NULL,pos=TRUE,tol=2,ruleset=NULL,keep_meta=TRUE, intval="into", convertRTMinute=F, numDigitsMZ=0, numDigitsRT=0, variableMetadataOutput="variableMetadata.tsv"){
 
-  #Load the two Rdata to extract the xset objects from positive and negative mode
-  cat("\tObject xset from positive mode\n")
-  print(xaP)
-  cat("\n")  
+    #Load the two Rdata to extract the xset objects from positive and negative mode
+    cat("\tObject xset from positive mode\n")
+    print(xaP)
+    cat("\n")
+
+    cat("\tObject xset from negative mode\n")
+    print(xaN)
+    cat("\n")
 
-  cat("\tObject xset from negative mode\n")
-  print(xaN)
-  cat("\n")
-    
-  cat("\n")
-  cat("\tCombining...\n")
-  #Convert the string to numeric for creating matrix
-  row=as.numeric(strsplit(ruleset,",")[[1]][1])
-  column=as.numeric(strsplit(ruleset,",")[[1]][2])
-  ruleset=cbind(row,column)
-  #Test if the file comes from an older version tool
-  if ((!is.null(xaP)) & (!is.null(xaN))) {
-    #Launch the combinexsannos function from CAMERA
-    cAnnot=combinexsAnnos(xaP, xaN,pos=pos,tol=tol,ruleset=ruleset)
-  } else {
-    stop("You must relauch the CAMERA.annotate step with the lastest version.")
-  }
-
-  
+    cat("\n")
+    cat("\tCombining...\n")
+    #Convert the string to numeric for creating matrix
+    row=as.numeric(strsplit(ruleset,",")[[1]][1])
+    column=as.numeric(strsplit(ruleset,",")[[1]][2])
+    ruleset=cbind(row,column)
+    #Test if the file comes from an older version tool
+    if ((!is.null(xaP)) & (!is.null(xaN))) {
+        #Launch the combinexsannos function from CAMERA
+        cAnnot=combinexsAnnos(xaP, xaN,pos=pos,tol=tol,ruleset=ruleset)
+    } else {
+        stop("You must relauch the CAMERA.annotate step with the lastest version.")
+    }
 
-  if(pos){
-    xa=xaP
-    listOFlistArgumentsP=listOFlistArguments
-    mode="neg. Mode"
-  } else {
-    xa=xaN
-    listOFlistArgumentsN=listOFlistArguments
-    mode="pos. Mode"
-  }
-  intval = "into"; for (steps in names(listOFlistArguments)) { if (!is.null(listOFlistArguments[[steps]]$intval)) intval = listOFlistArguments[[steps]]$intval }
-  peakList=getPeaklist(xa,intval=intval)
-  peakList=cbind(groupnames(xa@xcmsSet),peakList); colnames(peakList)[1] = c("name");
-  variableMetadata=cbind(peakList, cAnnot[, c("isotopes", "adduct", "pcgroup",mode)]);
-  variableMetadata=variableMetadata[,!(colnames(variableMetadata) %in% c(sampnames(xa@xcmsSet)))]
+    if(pos){
+        xa=xaP
+        listOFlistArgumentsP=listOFlistArguments
+        mode="neg. Mode"
+    } else {
+        xa=xaN
+        listOFlistArgumentsN=listOFlistArguments
+        mode="pos. Mode"
+    }
+
+    peakList=getPeaklist(xa,intval=intval)
+    peakList=cbind(groupnames(xa@xcmsSet),peakList); colnames(peakList)[1] = c("name");
+    variableMetadata=cbind(peakList, cAnnot[, c("isotopes", "adduct", "pcgroup",mode)]);
+    variableMetadata=variableMetadata[,!(colnames(variableMetadata) %in% c(sampnames(xa@xcmsSet)))]
+
+    #Test if there are more than two classes (conditions)
+    if ( nlevels(sampclass(xaP@xcmsSet))==2 & (!is.null(diffrepN)) & (!is.null(diffrepP))) {
+        diffrepP = diffrepP[,c("name","fold","tstat","pvalue")]; colnames(diffrepP) = paste("P.",colnames(diffrepP),sep="")
+        diffrepN = diffrepN[,c("name","fold","tstat","pvalue")]; colnames(diffrepN) = paste("N.",colnames(diffrepN),sep="")
 
-  #Test if there are more than two classes (conditions)
-  if ( nlevels(sampclass(xaP@xcmsSet))==2 & (!is.null(diffrepN)) & (!is.null(diffrepP))) {
-    diffrepP = diffrepP[,c("name","fold","tstat","pvalue")]; colnames(diffrepP) = paste("P.",colnames(diffrepP),sep="")
-    diffrepN = diffrepN[,c("name","fold","tstat","pvalue")]; colnames(diffrepN) = paste("N.",colnames(diffrepN),sep="")
-    
-    variableMetadata = merge(variableMetadata, diffrepP, by.x="name", by.y="P.name")
-    variableMetadata = merge(variableMetadata, diffrepN, by.x="name", by.y="N.name")
-  } 
-    
-  rownames(variableMetadata) = NULL
-  #TODO: checker
-  #colnames(variableMetadata)[1:2] = c("name","mz/rt");
+        variableMetadata = merge(variableMetadata, diffrepP, by.x="name", by.y="P.name")
+        variableMetadata = merge(variableMetadata, diffrepN, by.x="name", by.y="N.name")
+    }
+
+    rownames(variableMetadata) = NULL
+    #TODO: checker
+    #colnames(variableMetadata)[1:2] = c("name","mz/rt");
 
-  #If the user want to convert the retention times (seconds) into minutes.
-  if (listArguments[["convert_param"]]){
-    #converting the retention times (seconds) into minutes
-    cat("\tConverting the retention times into minutes\n")
-    variableMetadata$rtmed=cAnnot$rt/60; variableMetadata$rtmin=cAnnot$rtmin/60; variableMetadata$rtmax=cAnnot$rtmax/60;
-  }
+    variableMetadata = RTSecondToMinute(variableMetadata, convertRTMinute)
+    variableMetadata = formatIonIdentifiers(variableMetadata, numDigitsRT=numDigitsRT, numDigitsMZ=numDigitsMZ)
 
-  #If the user want to keep only the metabolites which match a difference
-  if(keep_meta){
-    variableMetadata=variableMetadata[variableMetadata[,c(mode)]!="",]
-  }
-      
-  #Write the output into a tsv file
-  write.table(variableMetadata, sep="\t", quote=FALSE, row.names=FALSE, file=variableMetadataOutput)
-  return(variableMetadata);
+    #If the user want to keep only the metabolites which match a difference
+    if(keep_meta){
+        variableMetadata=variableMetadata[variableMetadata[,c(mode)]!="",]
+    }
+
+    #Write the output into a tsv file
+    write.table(variableMetadata, sep="\t", quote=FALSE, row.names=FALSE, file=variableMetadataOutput)
+    return(variableMetadata);
 
 }