comparison lib.r @ 25:230f0bc9e792 draft

planemo upload commit a634879c0e651eb0eb7b435a107ee40cf30524fa

24:10176a940ec6

# lib.r version="2.2.1"

#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="")

  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)
}

#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')

  #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"]])

  #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"]])

  }  
  peakList=getPeaklist(xa,intval=listArguments[["intval"]])
  peakList=cbind(groupnames(xa@xcmsSet),peakList); colnames(peakList)[1] = c("name");


  # --- 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-")

          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)))])

          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=""))

          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), ]
    }
  }

  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)

  # --- 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));

}


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"){

  #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("\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)))]

  #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");

  #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;
  }

  #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);

}

25:230f0bc9e792

# 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="")

    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)
}

#@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)
}

#@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')

    #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"]])

    # 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"]])

    # 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)
            }

            listArguments4annotate = append(listArguments4annotate,
                list(rules=rulset))
        }
    }


    # 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-")

                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"]])

                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)))])

                diffrep = RTSecondToMinute(diffrep, listArguments[["convertRTMinute"]])
                diffrep = formatIonIdentifiers(diffrep, numDigitsRT=listArguments[["numDigitsRT"]], numDigitsMZ=listArguments[["numDigitsMZ"]])

                if(listArguments[["sortpval"]]){
                    diffrep=diffrep[order(diffrep$pvalue), ]
                }

                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))))]
    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
    write.table(variableMetadata, sep="\t", quote=FALSE, row.names=FALSE, file=variableMetadataOutput)

    return(list("xa"=xa,"diffrep"=diffrepOri,"variableMetadata"=variableMetadataOri));

}


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")

    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.")
    }

    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="")

        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 = 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);

}