R/report-utils.R

Defines functions .do.create.protein.info .create.or.load.my.protein.infos .create.or.load.quant.table .set .round.distr .create.or.load.ratiodistr .create.or.load.ptm.info .create.or.load.noise.model .create.or.load.ibspectra .load.property .check.property .exists.property property .get.property .create.or.load .get.or.load initialize.env load.properties .compile.tex create.reports testPdflatex

Documented in create.reports initialize.env load.properties property testPdflatex

testPdflatex <- function(pdflatex.cmd) {
  if (system(pdflatex.cmd) != 0) 
    stop("pdflatex does not seem to be installed. ",
         "Install LaTeX using the TeXLive or MiKTex distribution to generate ",
         "PDF reports with isobar. Set 'write.qc.report=FALSE' and 'write.report=FALSE', otherwise.")
}
    

create.reports <- function(properties.file="properties.R",
                           global.properties.file=system.file("report","properties.R",package="isobar"),
                           args=NULL,...,recreate.properties.env=TRUE,recreate.report.env=TRUE) {
  ow <- options("warn")
  if (!exists("properties.env") || recreate.properties.env) {
    properties.env <- load.properties(properties.file,
                                      global.properties.file,
                                      args=args,...)
    assign("properties.env",properties.env,envir=.GlobalEnv)
  }
  options(warn=properties.env[["warning.level"]])

  if (!exists("report.env") || recreate.report.env) {
    report.env <- .GlobalEnv
    initialize.env(report.env,properties.env)
  }

  zip.files <- c(properties.file)

  ## generate XLS report
  if(property('write.xls.report',properties.env)) {
    message("Writing isobar-analysis.xls")
    write.xls.report(properties.env,report.env)
    zip.files <- c(zip.files,"isobar-analysis.xls")
  }
 
  ## generate Latex/Sweave report
  if(property('write.qc.report',properties.env)) {
    message("Weaving isobar-qc report")
    Sweave(system.file("report","isobar-qc.Rnw",package="isobar"))
    if (property('use.name.for.report',properties.env)) {
      qc.name <- .sanitize.sh(sprintf("%s.qc",property('name',properties.env)))
    	file.rename("isobar-qc.tex",sprintf("%s.tex",qc.name))
    } else {
        qc.name <- "isobar-qc"
    }

    zip.files <- c(zip.files,sprintf("%s.tex",qc.name))
    if (properties.env[["compile"]]) 
      zip.files <- .compile.tex(qc.name,zip.files)
  }

  if(property('write.report',properties.env) && properties.env[["report.level"]] != 'peptide') {
    message("Weaving isobar-analysis report")
    name <- switch(properties.env[["report.level"]],
                   protein="isobar-analysis",
                   peptide="isobar-peptide-analysis",
                   stop(properties.env[["report.level"]]," report.level not known",
                        " - choose protein or peptide"))
    Sweave(system.file("report",paste(name,".Rnw",sep=""),package="isobar"))

    if (property('use.name.for.report',properties.env)) {
    	tex.name <- sprintf("%s.tex",name)
      name <- .sanitize.sh(sprintf("%s.quant",property('name',properties.env)))
    	file.rename(tex.name,sprintf("%s.tex",name))
    } else {
    }


    zip.files <- c(zip.files,sprintf("%s.tex",name))
    if (properties.env[["compile"]])
      zip.files <- .compile.tex(name,zip.files)
  }

  if (properties.env[["zip"]]) {
    zip.f <- sprintf("%s.zip",property('name',properties.env))
    zip(zip.f,zip.files)
    message("Created zip archive ",zip.f)
  }

  options(ow) 
  message("\nSUCCESSFULLY CREATED REPORTS\n")
}

  .compile.tex <- function(name,zip.files) {
    r.cmd = shQuote(file.path(R.home("bin"),"R"))
    testPdflatex(paste(r.cmd,"CMD","pdflatex -version"))
    
    dir <- tempdir()
    cat("compiling ",name,".tex ...  1",sep="")
    .call.cmd(sprintf("%s CMD pdflatex -halt-on-error -output-directory=%s %s.tex",r.cmd,dir,name),
              paste(dir,"/",basename(name),".stdout",sep=""))
    cat(" 2")
    .call.cmd(sprintf("%s CMD pdflatex -halt-on-error -output-directory=%s %s.tex",r.cmd,dir,name),
              paste(dir,"/",basename(name),".stdout",sep=""))
    cat(" done!\n\n")
    file.copy(file.path(dir,paste0(name,".pdf")),file.path(getwd(),paste0(name,".pdf")),overwrite=TRUE)
    file.remove(file.path(dir,paste0(name,".pdf")))
    c(zip.files,sprintf("%s.pdf",name))
  }
 
#- load properties

load.properties <- function(properties.file="properties.R",
                            global.properties.file=system.file("report","properties.R",package="isobar"),
                            args=NULL,...) {

  properties.env <- new.env()
  tmp.properties.env <- new.env()

  message("Loading global properties file ",global.properties.file," ...")
  tryCatch(sys.source(global.properties.file,properties.env),
           error=function(e) stop("\nCould not read properties file:\n", e) )

  message("Looking for local properties file ",properties.file," ...")
  if (!is.null(properties.file) && file.exists(properties.file)) {
    message("  Loading local properties file ...")
    ## load properties file
    tryCatch(sys.source(properties.file,tmp.properties.env),
             error=function(e) stop("\nCould not read properties file:\n", e) )
    .env.copy(properties.env,tmp.properties.env)
    properties.env[["properties.file"]] <- properties.file
    properties.env[["properties.file.content"]] <- readLines(properties.file)
  } else {
    message("  No local properties file.")
  }
  
  ## command argument parsing
  tmp.properties.env <- new.env()
  if (length(args) > 0)
    message("parsing command line arguments ...")
  for (arg in args) {
    if (grepl("^--",arg)) {

      arg.n.val <- strsplit(substring(arg,3),"=")[[1]]
      if (length(arg.n.val) == 1)
        tmp.properties.env[[arg.n.val]] <- TRUE
      else if (length(arg.n.val) == 2) 
        tmp.properties.env[[arg.n.val[1]]] <- switch(arg.n.val[2],'TRUE'=TRUE,'FALSE'=FALSE,arg.n.val[2])
      else
        stop("Could not parse command line argument ",arg)
    } else {
        stop("Could not parse command line argument ",arg)
    }
    .env.copy(properties.env,tmp.properties.env)
  }

  ## ... parsing
  dotargs <- list(...)
  if (length(dotargs) > 0) {
    message("parsing function arguments")
    .env.copy(properties.env,list2env(dotargs))
  }

  return(properties.env)
}

#- initialize environment
initialize.env <- function(env,properties.env) {
  ## get property and exists property convenience functions
  get.property <- function(name) .get.property(name,properties.env)

  if (file.exists(get.property('cachedir'))) {
    if (!file.info(get.property('cachedir'))$isdir)
      stop("Cannot write into cachedir [",get.property('cachedir'),"]",
           " - it is a file and no directory")
  } else {
    ret <- dir.create(get.property('cachedir'))
    if (!ret) stop("Error creating cachedir [",get.property('cachedir'),"]")
  }

  ib.name <- file.path(.get.property('cachedir',properties.env),"ibspectra.rda")
  if (file.exists(ib.name)) {
    load(ib.name)
    env[["ibspectra"]] <- get("ibspectra")
  } else {
    env[["ibspectra"]] <- .create.or.load.ibspectra(properties.env)
    save(list='ibspectra',envir=env,file=ib.name)
  }
  if ("site.probs" %in% colnames(fData(env[["ibspectra"]])) 
      && ! "pep.siteprobs" %in% colnames(fData(env[["ibspectra"]])))
    fData(env[["ibspectra"]])$pep.siteprobs <- 
      .convertPhosphoRSPepProb(fData(env[["ibspectra"]])[,'peptide'],fData(env[["ibspectra"]])[,'site.probs'],
                               round.to.frac=20)


  env[["noise.model"]] <- .create.or.load.noise.model(env,properties.env)
  env[["ratiodistr"]] <- .create.or.load.ratiodistr(env,properties.env)
  if (identical(properties.env[["report.level"]],"peptide") ) {
    env[["ptm.info"]]  <- .create.or.load.ptm.info(properties.env,proteinGroup(env[["ibspectra"]]))
    if (all(c('use.for.quant','pep.siteprobs') %in% colnames(fData(env[['ibspectra']]))) && !all(fData(env[['ibspectra']])[['use.for.quant']]))
      env[["quant.tbl.notlocalized"]] <- 
        .create.or.load.quant.table(env,properties.env,name="quant.tbl.notlocalized",type="other-sites")
  }
  env[["quant.tbl"]] <- .create.or.load.quant.table(env,properties.env)
  if (!"ac" %in% colnames(env[["quant.tbl"]]) && "protein" %in% colnames(env[["quant.tbl"]]))
    env[["quant.tbl"]][,'ac'] <- env[["quant.tbl"]]$protein

  ## required for TeX
  if (property('write.report',properties.env) && identical(properties.env[["report.level"]],"protein"))
    env[["my.protein.infos"]] <- .create.or.load.my.protein.infos(env,properties.env)

  if (property('write.xls.report',properties.env)) {
    xls.table.name <- paste0("xls.quant.tbl.",properties.env[["xls.report.format"]])
    env[["xls.quant.tbl"]] <- 
      .create.or.load.xls.quant.tbl(env,properties.env,xls.table.name,"quant.tbl")
    if (exists("quant.tbl.notlocalized",envir=env)) {
      xls.ntable.name <- paste0("xls.quant.tbl.notlocalized.",properties.env[["xls.report.format"]])
      env[["xls.quant.tbl.notlocalized"]] <- 
        .create.or.load.xls.quant.tbl(env,properties.env,xls.ntable.name,"quant.tbl.notlocalized")
    }
  }
}

#- property loading helper functions
.DOES.NOT.EXIST = "NOT AVAILABLE"

.get.or.load <- function(name,envir,msg.f=name,class=NULL,null.ok=FALSE,do.load=FALSE) {
  if (.exists.property(name,envir,null.ok=null.ok)) {
    o <- .get.property(name,envir)
    if (is.null(o) && null.ok) return(NULL)
    if (!is.null(class) && is(o,class)) {
      return(o)
    } else if (is(o,"character")) {
      file.name <- o
    } else {
      stop("property ",name," is neither 'character' nor of a specified class")
    }
  } else {
    file.name <- sprintf("%s/%s.rda",.get.property('cachedir',envir),name)
  }
  if (file.exists(file.name) && (!envir[["regen"]] || do.load)) {
    message(sprintf("  loading %s from %s ... ",msg.f, file.name),appendLF=FALSE)
    x <- .load.property(name,file.name)
    message("done")
    return(x)
  } else {
    stop("Cannot get or load property [",name,"] - would expect it in [",file.name,"], but file does not exist.")
  }
}

.create.or.load <- function(name,envir,f,msg.f=name,regenerate=FALSE,do.load=FALSE,class=NULL,error=stop,default.value=NULL,...) {
  x <- if ( !regenerate ) tryCatch(.get.or.load(name,envir,msg.f,class,do.load=do.load),error=function(e) .DOES.NOT.EXIST) else .DOES.NOT.EXIST
  if (identical(x,.DOES.NOT.EXIST)) {
    message(paste("  creating",msg.f,"... "),appendLF=FALSE)
    tryCatch({
      x <- f(...)
      assign(name,x)
      file.name <- sprintf("%s/%s.rda",.get.property('cachedir',envir),name)
      message("done ",appendLF=FALSE)
      save(list=c(name),file=file.name)
      message("& saved to ",file.name)

    },error=error)

    if (!exists(name,inherits=FALSE)) {
      x <- default.value
      assign(name,x)
    }
  }
  if (!is.null(class) && !is(x,class))
    stop("property [",name,"] should be of class [",class,"] but is of class [",class(x),"]")

  return(x)
}
 
.get.property <- function(x,envir) get(x,envir=envir,inherits=FALSE)

property <- function(x, envir, null.ok=TRUE,class=NULL) {
  if (!.exists.property(x,envir)) 
    stop("property ",x," does not exist!")

  o <- get(x,envir=envir,inherits=FALSE)

  if (!null.ok && is.null(o))
    stop("property ",x," may not be NULL, but it is!")

  if (!is.null(class) && !is(o,class))
    stop("property ",x," does not have class ",class,"!")

  return(o)
}

.exists.property <- function(x,envir,null.ok=TRUE) 
  exists(x,envir=envir,inherits=FALSE) &&
    (null.ok || !is.null(.get.property(x,envir)))

.check.property <- function(x,envir,inherits=FALSE,msg="",print=TRUE,valid=NULL,check.file=FALSE,...) {
  def <- grep(x,envir[["properties.file.content"]],value=TRUE)
  if (length(def) > 0)
    def <- paste("\n  Corresponding line in ",envir[["properties.file"]],":\n\t",
                 paste(def,collapse="\n\t"),"\n\n",sep="")

  p <- .get.property(x,envir=envir)
  does.not.exist <- !exists(x,envir=envir,inherits=inherits,...)
  is.null.p <- does.not.exist || is.null(.get.property(x,envir=envir))
  length.0 <- does.not.exist || length(.get.property(x,envir=envir))==0
  if (does.not.exist || is.null.p || length.0) 
    stop ("  property '",x,"' not defined or did not evaluate in the properties file, but it is required",
          ifelse(is.null(valid),"",paste(" to be one of \n\t",paste(valid,collapse="\n\t"),
          sep="")),def)
  
  if (check.file && !file.exists(.get.property(x,envir=envir)))
    stop("  property '",x,"' should be assigned to a file, but it is not:",def)

  if (!is.null(valid)) {
    isnt.valid <- !p %in% valid
    if (isnt.valid) stop(" property '",x,"' must be one of \n\t",paste(valid,collapse="\n\t"),",\n\n",
                         "  but it is ",p )
  }
  if (print) {
    message("  property '",x,"' defined",ifelse(is.character(p),paste0(": ",paste(p,collapse="; ")),"."))
  }
}

.load.property <- function(name,file) {
  tmp.env <- new.env()
  load(file,envir=tmp.env)
  if (exists(name,envir=tmp.env))
    return(get(name,envir=tmp.env))
  else if (length(ls(envir=tmp.env)) == 1)
    return(get(ls(envir=tmp.env)[1],envir=tmp.env))
  else
    stop("Could not load property ",name," from file ",file)
}

.create.or.load.ibspectra <- function(properties.env) {
  get.property <- function(name) .get.property(name,properties.env)
  check.property <- function(name,...) .check.property(name,properties.env,...)
  ## check that required properties are defined
  check.property('ibspectra')
  check.property('type',valid=IBSpectraTypes())

  readIBSpectra.args <- get.property('readIBSpectra.args')
  if (is.null(readIBSpectra.args)) readIBSpectra.args <- list()
  for (name in names(readIBSpectra.args)) {
      arg <- readIBSpectra.args[[name]]
      if (!is.null(arg)) {
          if (!is.null(names(arg)))
            arg <- paste(names(arg),arg,collapse="=")
          if (is.function(arg)) { message("    ",name,": function")
          } else {
          message("    ",name,": ",
                  paste(arg,collapse=ifelse(length(arg)>2,"\n\t",", ")))
          }
      }
  }
  readIBSpectra.args[["type"]]=get.property('type')
  readIBSpectra.args[["id.file"]]=get.property('ibspectra')
  readIBSpectra.args[["fragment.precision"]]=get.property('fragment.precision')
  readIBSpectra.args[["fragment.outlier.prob"]]=get.property('fragment.outlier.prob')
  readIBSpectra.args[["proteinGroupTemplate"]]=.get.or.load('protein.group.template',properties.env,"ProteinGroup",null.ok=TRUE)

  if (is.data.frame(get.property('ibspectra')) || all(file.exists(get.property('ibspectra')))) {
    if (is.data.frame(get.property('ibspectra'))) {
      readIBSpectra.args[["id.file"]] <- get.property('ibspectra')
      ibspectra <- do.call(readIBSpectra,readIBSpectra.args)
    } else if (grepl(".csv",get.property('ibspectra'))) {
        message("ibspectra ends on .csv:\n",
                sprintf('ibspectra <- readIBSpectra("%s",%s) ...',
                        get.property('type'),
                        ifelse(length(get.property('ibspectra'))==1,
                               paste0('"',get.property('ibspectra'),'"'),
                               paste0('c("',get.property('ibspectra'),'")',collapse='","')
                               )))
        ibspectra <- do.call(readIBSpectra,readIBSpectra.args)
    } else if (grepl(".rda",get.property("ibspectra"))) {
        message("ibspectra ends on .rda: ","loading")
        ibspectra <- .load.property("ibspectra",get.property("ibspectra"))
    } else {
        stop("weird naming of ibspectra: ",get.property("ibspectra"),". use XXX.csv or XXX.rda")
    }

  } else {
    # create ibspectra file
    message(sprintf('\n  file %s does not exist, creating ibspectra',get.property('ibspectra')))
    check.property('peaklist',print=FALSE,check.file=TRUE)
    check.property('identifications',print=FALSE,check.file=TRUE)
    message(sprintf('    id.files = %s%s\n    peaklist.files = %s%s\n',
          ifelse(length(get.property('identifications'))>1,"\n      ",""),
          paste(get.property('identifications'),collapse="\n      "),
          ifelse(length(get.property('peaklist'))>1,"\n      ",""),
          paste(get.property('peaklist'),collapse="\n      ")))

    readIBSpectra.args[["id.file"]]=get.property('identifications')
    readIBSpectra.args[["peaklist.file"]]=get.property('peaklist')

    ibspectra <- do.call(readIBSpectra,readIBSpectra.args)
    if (grepl(".csv",get.property('ibspectra'))) {
        write.table(as.data.frame(ibspectra),sep="\t",row.names=F,file=get.property('ibspectra'),quote=FALSE)
    } else if (grepl(".rda",get.property("ibspectra"))) {
        save(ibspectra,file=get.property("ibspectra"),compress=TRUE)
    }  else {
        stop("weird naming of ibspectra: ",get.property("ibspectra"),". use XXX.csv or XXX.rda")
    }
  }

  if (!is.null(get.property("isotope.impurities"))) {
    isotopeImpurities(ibspectra) <- get.property("isotope.impurities")
  }

  if (property('correct.isotope.impurities',properties.env))
    ibspectra <- correctIsotopeImpurities(ibspectra)
  if (property('normalize',properties.env)) {
    ibspectra <-
      normalize(ibspectra,
                use.protein=property('normalize.use.protein',properties.env),
                exclude.protein=property('normalize.exclude.protein',properties.env),
                peptide.specificity=property('peptide.specificity',properties.env),
                f=property('normalize.function',properties.env),
                channels=property('normalize.channels',properties.env),
                na.rm=property('normalize.na.rm',properties.env),
                normalize.factors=property('normalize.factors',properties.env))
  }

  class.labels <- LETTERS[1:length(reporterTagNames(ibspectra))]
  if (.exists.property('class.labels',properties.env,null.ok=FALSE))
    class.labels <- get.property('class.labels')
  if (!any(table(class.labels)>1) && property('summarize',properties.env)) {
    stop("When summarize=TRUE, the must be more then one channel per class")
  }
  if (!is.character(class.labels)) {
    stop("Please provide class.labels of class character!")
  }
  classLabels(ibspectra) <- class.labels

  if (!is.null(property('protein.info.f',properties.env)))
    tryCatch({
    proteinInfo(proteinGroup(ibspectra)) <- 
      .create.or.load("protein.info",envir=properties.env,
                      f=property('protein.info.f',properties.env),
                      x=proteinGroup(ibspectra),
                      do.load=TRUE, msg.f="protein.info",
                      error=warning,default.value=proteinInfo(proteinGroup(ibspectra)))
    },error=function(e) stop("Error creating proteinInfo using function defined in [protein.info.f]: ",e,"Set protein.info.f to NULL if you want to create a report anyway."))

  if ("site.probs" %in% colnames(fData(ibspectra)) 
      && ! "pep.siteprobs" %in% colnames(fData(ibspectra)))
     fData(ibspectra)$pep.siteprobs <- 
        .convertPhosphoRSPepProb(fData(ibspectra)[,'peptide'],fData(ibspectra)[,'site.probs'],
                                 round.to.frac=20)
 
  return(ibspectra)
}

.create.or.load.noise.model <- function(env,properties.env) {
  noise.model.channels <- .get.property("noise.model.channels",properties.env)
  if (is.null(noise.model.channels)) 
    noise.model.channels <- reporterTagNames(env[["ibspectra"]])[!is.na(classLabels(env[["ibspectra"]]))]
  is.one.to.one <- isTRUE(.get.property("noise.model.is.technicalreplicates",properties.env))

  noise.model <- .get.property("noise.model",properties.env)
  if (!is(noise.model,"NoiseModel")) {
    noise.model.f <- .get.property("noise.model",properties.env)
    if (!file.exists(noise.model.f)) {
      message("estimating noise model as ",ifelse(is.one.to.one,"","non"),
              " one-to-one from channels ",paste0(noise.model.channels,collapse=", ")," ...")
      noise.model <- new("ExponentialNoiseModel",env[["ibspectra"]],one.to.one=is.one.to.one,
                         reporterTagNames=noise.model.channels,
                         min.spectra=property('noise.model.minspectra',properties.env))
      save(noise.model,file=noise.model.f,compress=TRUE)
    } else {
      message(sprintf("  loading noise model from %s ...",noise.model.f))
      noise.model <- .load.property("noise.model",noise.model.f)
    }
  }
  return(noise.model)
}

.create.or.load.ptm.info <- function(properties.env,protein.group) {
  if (is.null(property('ptm.info.f',properties.env))) 
    properties.env[["ptm.info.f"]] <- getPtmInfoFromNextprot

  return(.create.or.load("ptm.info",envir=properties.env,
                         f=property('ptm.info.f',properties.env),
                         protein.group=protein.group))
}

.create.or.load.ratiodistr <- function(env,properties.env) {
  
  return(.create.or.load("ratiodistr",envir=properties.env,class="Distribution",
                         msg.f="biological variability ratio distribution",f=function(){

    cl <- classLabels(env[["ibspectra"]])
    if (!is.null(property('ratiodistr.class.labels',properties.env)))
      cl <- property('ratiodistr.class.labels',properties.env)

    ratios.for.distr.fitting <- .create.or.load.ratiodistr.ratios(env,properties.env,cl)

    if (all(is.na(ratios.for.distr.fitting[,'lratio']))) 
      stop("All ratios for distr fitting are NA - are the correct class labels used?")

    if (!is.function(property('ratiodistr.fitting.f',properties.env)))
      stop("ratiodistr.fitting.f must be set to a function [e.g. fitCauchy or fitTd]")

    ratiodistr <- property('ratiodistr.fitting.f',properties.env)(ratios.for.distr.fitting[,'lratio'])
    ratiodistr <- .round.distr(ratiodistr,digits=5)
    attr(ratiodistr,"combn.method") <- attr(ratios.for.distr.fitting,"combn.method")
    attr(ratiodistr,"cl") <- attr(ratios.for.distr.fitting,"cl")
    attr(ratiodistr,"tagNames") <- reporterTagNames(env[["ibspectra"]])
    ratiodistr
  }))
}

.round.distr <- function(distr,digits) {
  for (s in slotNames(param(distr))) 
    if (s != "name" && is.numeric(slot(param(distr),s))) 
      slot(distr@param,s) <- round(slot(param(distr),s),digits)
  distr
}

.create.or.load.ratiodistr.ratios <- function (env,properties.env,cl) {
  .create.or.load("ratios.for.distr.fitting",envir=properties.env,class="data.frame",
                  msg.f="ratios for biological variability distribution fitting",f=function() {

    message("  Using class labels ",paste(cl,collapse=", "))
    if (property('summarize',properties.env) || any(table(cl)>=2)) {
      method <- "intraclass"
    } else {
      message(" WARNING: ratiodistr will be computed based on global ratios")
      method <- "global"
    }
    if (any(table(properties.env[["class.labels"]])>2))
      do.summarize <- properties.env[["summarize"]]
    else
      do.summarize <- FALSE

    if (identical(properties.env[["report.level"]],"peptide"))
      ratios.for.distr.fitting <- peptideRatios(env[["ibspectra"]],noise.model=env[["noise.model"]],do.warn=FALSE,
                                  cl=cl,combn.method=method,symmetry=isTRUE(properties.env[["ratiodistr.symmetry"]]),summarize=do.summarize)
    else
      ratios.for.distr.fitting <- proteinRatios(env[["ibspectra"]],noise.model=env[["noise.model"]],do.warn=FALSE,
                                      cl=cl,combn.method=method,symmetry=isTRUE(properties.env[["ratiodistr.symmetry"]]),summarize=do.summarize)

    if (all(is.nan(ratios.for.distr.fitting[["lratio"]])))
      stop("Cannot compute protein ratio distribution - no ratios available.\n",
           "Probably due to missing reporter intensities.")

    attr(ratios.for.distr.fitting,"combn.method") <- method
    attr(ratios.for.distr.fitting,"cl") <- cl

    ratios.for.distr.fitting
  })
}


.set <- function(x,name,list) {
  if (name %in% names(list)) {
    stop(name," already assigned in list!")
  }
  x
}

.create.or.load.quant.table <- function(env,properties.env,name="quant.tbl",type='confident-sites') {
  protein.group <- proteinGroup(env[["ibspectra"]])
  protein.info <- proteinInfo(protein.group)
  isoforms <- protein.group@isoformToGeneProduct
  
  .create.or.load(name,envir=properties.env,class="data.frame",
                  msg.f=paste("table of ratios of",properties.env[["report.level"]],"as",name),f=function() {
    if (!is.null(property('ratios.opts',properties.env)$summarize)) {
        message("WARNING: ratio.opts$summarize will be overwritten,",
                " define it outside of ratio.opts!")
        warning("ratio.opts$summarize will be overwritten,",
                " define it outside of ratio.opts!")
    }
    ratios.opts <- property('ratios.opts',properties.env)
    set.ratioopts <- function(x,name=names(x)) {
      if (is.null(intersect(name,names(ratios.opts))))
          stop("property '",intersect(name,names(ratios.opts)),
               "' already assigned in ratios.opts list - check your properties file!")

      if (is.list(x))
        ratios.opts <<- c(ratios.opts,x)
      else
        ratios.opts[[name]] <<- x
    }

    set.ratioopts(name="ibspectra",env[["ibspectra"]])
    set.ratioopts(name="noise.model",env[["noise.model"]])
    set.ratioopts(name="ratiodistr",env[["ratiodistr"]])
    
    if(identical(properties.env[["report.level"]],"peptide")) {
      if (!"use.for.quant" %in% colnames(fData(env[["ibspectra"]])))
        fData(env[["ibspectra"]])[["use.for.quant"]] <- TRUE
      if (type=='confident-sites')
        pep.n.modif <- unique(apply(fData(env[["ibspectra"]])[fData(env[["ibspectra"]])[["use.for.quant"]],
                                    c("peptide","modif")],2,cbind))
      else {
        pep.n.modif <- unique(apply(fData(env[["ibspectra"]])[!fData(env[["ibspectra"]])[["use.for.quant"]],
                                    c("peptide","modif","pep.siteprobs")],2,cbind))
        set.ratioopts(list(use.for.quant.only=FALSE))
      }


      set.ratioopts(list(peptide=pep.n.modif,
                         proteins=NULL))

    } else if (identical(properties.env[["report.level"]],"protein")) {
      set.ratioopts(list(peptide=NULL,
                         proteins=reporterProteins(proteinGroup(env[["ibspectra"]])),
                         quant.w.grouppeptides=property('quant.w.grouppeptides',properties.env)))
    } else {
      stop("don't known level ",properties.env[["report.level"]])
    }

    if (is.null(property('cmbn',properties.env)) & 
	!is.null(property('vs.class',properties.env)))
      properties.env[["cmbn"]] <- combn.matrix(reporterTagNames(env[["ibspectra"]]),
					   "versus.class",
					   property('class.labels',properties.env),
					   vs=property('vs.class',properties.env))

    set.ratioopts(list(combn.method=property('combn.method',properties.env),
                       cl=classLabels(env[["ibspectra"]]),
                       summarize=property('summarize',properties.env),
                       cmbn=property('cmbn',properties.env),
                       use.na=property('use.na',properties.env),
                       zscore.threshold=properties.env[["zscore.threshold"]],
                       do.warn=FALSE))

    if (!is.null(property('correct.peptide.ratios.with',properties.env))) {
      protein.quant.tbl <- .get.or.load("correct.peptide.ratios.with",properties.env)
      correct.protein.group <- if ( !is.null(property('correct.peptide.ratios.with_protein.group',properties.env)) ) {
        .get.or.load("correct.peptide.ratios.with_protein.group",properties.env)
      } else {
        warning( '"correct.peptide.ratios.with_protein.group" property not found, ',
                 'using the existing protein group' )
        .get.or.load("protein.group.template",properties.env)
      }
      ratios.opts[["before.summarize.f"]] <- function(...)
        correct.peptide.ratios(..., protein.quant.tbl=protein.quant.tbl,
                               protein.group.combined = correct.protein.group,
                               correlation = property('peptide.protein.correlation',properties.env))
    }

    quant.tbl <- do.call("proteinRatios",ratios.opts)

    quant.tbl[,"sd"] <- sqrt(quant.tbl[,"variance"])
    
#    quant.tbl$sign.string <- "not significant"
#    quant.tbl$sign.string[quant.tbl$is.significant] <- "is significant"
    
#    if (length(property('preselected',properties.env)) > 0) {
#      preselected <- unique(ip[sub("-.*","",names(ip)) %in% property('preselected',properties.env)])
#      quant.tbl$is.preselected <- quant.tbl$protein %in% preselected

#      quant.tbl$sign.string <- "not significant"
#      quant.tbl$sign.string[is.sign & !quant.tbl$is.preselected] <- "is significant [ours]"
#      quant.tbl$sign.string[is.sign & quant.tbl$is.preselected] <- "is significant [both]"
#      quant.tbl$sign.string[!is.sign & quant.tbl$is.preselected] <- "is significant [theirs]"
#    }

    if (identical(properties.env[["report.level"]],"protein")) {
      quant.tbl[,"gene_names"] <- sapply(quant.tbl[,"ac"], function(x) {
        if (length(protein.info) == 0) return("")
        allreporter <- indistinguishableProteins(protein.group,protein.g=x)
        acs <- unique(isoforms[allreporter,"proteinac.wo.splicevariant"])
        paste(sort(unique(protein.info[protein.info[["accession"]] %in% acs,"gene_name"])),
              collapse=", ")
      })
      sort.genenames <- quant.tbl[,"gene_names"]
      sort.genenames[sort.genenames==""] <- quant.tbl[sort.genenames=="","ac"]
      quant.tbl <- quant.tbl[order(sort.genenames,quant.tbl[,"r1"],quant.tbl[,"r2"]),]
      quant.tbl[,"group"] <- as.numeric(factor(quant.tbl[,"ac"],levels=unique(quant.tbl[,"ac"])))
    }
  
    if (all(is.na(quant.tbl[["lratio"]])))
      stop("All ratios are NA")

    return(quant.tbl)
  })
}

.create.or.load.my.protein.infos <- function(env,properties.env) {
  .create.or.load("my.protein.infos",envir=properties.env,
                  f=function() {

    protein.groupnames <-unique(env[["quant.tbl"]][,"ac"])
    ## if (is.null(protein.info)) { stop("protein info is null!")}                
    my.protein.infos <- llply(protein.groupnames, .do.create.protein.info, protein.group=proteinGroup(env[["ibspectra"]]), 
                              .parallel=isTRUE(getOption('isobar.parallel')))
    #my.protein.infos <- lapply(protein.groupnames, .do.create.protein.info, protein.group=protein.group)
    names(my.protein.infos) <- protein.groupnames
    return(my.protein.infos)
  })
}
  
.do.create.protein.info <- function(x,protein.group) {
    protein.group.table <- proteinGroupTable(protein.group)
      allgroupmember <- indistinguishableProteins(protein.group, protein.g =
                                                  protein.group.table[protein.group.table[,'reporter.protein'] %in% x,"protein.g"])
     
      reporter.protein.info <- my.protein.info(protein.group,x)
      collapsed.gene_name <- human.protein.names(reporter.protein.info)
      peptides <- peptides(protein.group,protein=x,do.warn=FALSE)
      peptides.gs <- peptides(protein.group,protein=x,
                              specificity=GROUPSPECIFIC,do.warn=FALSE)
      peptides.rs <- peptides(protein.group,protein=x,
                              specificity=REPORTERSPECIFIC,do.warn=FALSE)
      n.spectra <- length(names(spectrumToPeptide(protein.group))[spectrumToPeptide(protein.group)%in%peptides])
      
      tbl.protein.name <- sort(collapsed.gene_name[,'protein_name'])[1];
      if (length(unique(collapsed.gene_name[,'protein_name'])) > 1)
        tbl.protein.name <- paste(tbl.protein.name,", ...",sep="")
      
      list(n.reporter = nrow(reporter.protein.info),
           n.groupmember = length(allgroupmember),
           reporter.protein.info = reporter.protein.info,
           n.peptides=length(peptides),
           n.spectra=n.spectra,
           collapsed.gene_name = collapsed.gene_name,
           table.name = ifelse(
             length(collapsed.gene_name[,'ac_link'])>3,
             paste(paste(collapsed.gene_name[1:3,'ac_link'],collapse=", "),
                   ", \\dots",sep=""),
             paste(collapsed.gene_name[,'ac_link'],collapse=", ")),
           section.name = sanitize(paste(unique(collapsed.gene_name[,'name_nolink']),
             collapse=", ")),
           table.protein.name = tbl.protein.name,
           gene.name = paste(sort(unique(reporter.protein.info[,'gene_name'])),collapse=", ")
           )
    }

## copys objects from env into parentenv.
## Those objects MUST exist in parentenv before.
.env.copy <- function (parentenv,env) {
  printval <- function(x) if(is.character(x)) capture.output(dput(x)) else toString(x)

  for (object in ls(envir=env)) {
    if (exists(object,envir=parentenv)) {
      ## assign object to parent env
      val <- get(object,envir=env)
      if (is.list(val)) {
        val.s <- paste0("list(",paste(names(val),sapply(val,printval),sep=": ",collapse=";"),")")
      } else {
        val.s <- tryCatch(printval(val),error=function(e) capture.output(print(val)))
      }
      message(sprintf("%30s = %s",object,paste(val.s,collapse=paste0("\n",paste(rep(" ",35),collapse="")))))
      assign(object, value=get(object,envir=env),envir=parentenv)
      ## remove(object,envir=env)
    } else {
      stop("Illegal property ",object,"!\n",
           "Call script with --help for usage details.\n\n",
           "Available properties:\n",
           "\t",paste(ls(envir=parentenv),collapse="\n\t"),"\n")
                      
      ## more verbose list:
      ##cat(paste(unlist(lapply(objects(),function(o) {sprintf("%s [%s]",o,class(get(o)))})),collapse="\n\t"))
    }
  }
}

Try the isobar package in your browser

Any scripts or data that you put into this service are public.

isobar documentation built on Nov. 8, 2020, 7:48 p.m.