R/evalScoring.R

In macat: MicroArray Chromosome Analysis Tool

evalScoring <- function(data,class,chromosome,nperms=1000,permute="labels", pcompute="empirical", subset=NULL, newlabels=NULL,kernel=rbf,kernelparams=NULL,cross.validate=TRUE,paramMultipliers=2^(-4:4),ncross=10,step.width=100000,memory.limit=TRUE, verbose=TRUE){
  # evaluate scores by permutations
  # data: list generated by Loader.R
  # class: tumor class to evaluate scoring for, has to occur in 'labels'
  # chromosome: chromosome to investigate
  # nperms: number of permutations
  # permute: permutation mode for null model, one of "labels","locations"
  # subset: if only subset of samples is used, enter the column-indices of these in the expression matrix here
  # step.width: number of base pairs between two stops while computing sliding scores
  # memory.limit: is memory < 2GB ?

  ### check arguments: ###
  stopifnot(!is.null(data$expr),!is.null(data$labels),
            !is.null(data$geneName),!is.null(data$geneLocation),
            !is.null(data$chromosome),!is.null(data$chip))
  if(is.null(subset))
    subset <- seq(ncol(data$expr))
  if(!is.null(newlabels)){
    if(length(newlabels)!=length(subset))
      stop("Length of new labels vector does not match number of samples!\nConsider using the 'subset' argument!\n")
    labels <- newlabels
  } else labels <- data$labels[subset]
  if (!(class %in% names(table(labels))))
    stop(paste(class,"is no class in labels vector!\n"))
  stopifnot(chromosome %in% unique(data$chromosome))
  permute = match.arg(permute, c("locations","labels"))

  ### interpret arguments: ### 
  areOfClass <- as.numeric(labels == class) #binary labels
  if (verbose)
    cat(paste("Investigating",sum(areOfClass),"samples of class",class,"...\n"))
  
  n.possible.permutations <- choose(length(areOfClass),sum(areOfClass))
  if ((permute=="labels") && (n.possible.permutations < nperms))
    warning(paste("For these class labels, the number of possible permutations is",n.possible.permutations,", less than the selected number of permuations,",nperms,".\n Consider fewer permutations or select 'permute=\"locations\"'\n"))
  
  # Assess list components:
  geneID <- data$geneName
  geneLoc <- data$geneLocation
  chrom <- data$chromosome
  X <- data$expr[,subset,drop=FALSE]
  if (is.null(kernelparams)){ # set default kernel parameters based on chosen kernel function
    kernelparams <- fitkernelparams(data,chromosome,kernel)
  } #then (is.null(kernelparams))
  
  # compute scores for all genes on chip:  
  if (permute=="labels"){
    samResult <- scoring(X,areOfClass,method="SAM", pcompute=pcompute, nperms=nperms, verbose=verbose)
  } else {
    samResult <- scoring(X,areOfClass,method="SAM",pcompute="none", verbose=verbose)
  } #else
  classScores <- matrix(samResult$observed,ncol=1)
  rownames(classScores) <- rownames(X)
  # now focus only on one chromosome:
  onChromIndex <- which(chrom==chromosome)
  onChromGenes <- geneID[onChromIndex]
  onChromLocs <- geneLoc[onChromIndex]
  absOnChromLocs <- abs(onChromLocs)
  onChromChrom <- chrom[onChromIndex]
  onChromScores <- classScores[onChromGenes,,drop=FALSE]
  onChromPvalues <- samResult$pvalues[onChromGenes]

  if (cross.validate){
    if (verbose)
      cat("Performing cross-validation to obtain optimal parameters...\n")
    evp <- evalParams(locations=absOnChromLocs,scores=onChromScores,kernel=kernel,kernelparams=kernelparams,paramMultipliers=paramMultipliers,ncross=10,verbose=verbose)
    kernelparams <- evp$best
  } # if(cross.validate

  if (verbose)
    cat("Computing sliding values for scores...\n")
  steps <- getsteps(absOnChromLocs,step.width)
  kernelweights <- kernelmatrix(steps,absOnChromLocs,kernel,kernelparams)
  classSlideScores <- kernelize(onChromScores,kernelweights)
    
  ### assess scores by permutations (two methods): ###
  if (permute=="locations"){
    # do permutations:
    if (verbose)
      cat("Building permutation matrix...\n")
    permMatrix <- matrix(nrow=nperms,ncol=length(absOnChromLocs))
    for (i in 1:nperms)
      permMatrix[i,] <- sample(absOnChromLocs)
    # auxiliary function:
    computePermSlideScores <- function(permLocs, verbose=FALSE){
      if (verbose) cat(".")
      permM <- list(geneName=onChromGenes,geneLocation=permLocs,
                    chromosome=onChromChrom,expr=onChromScores)
      permS <- compute.sliding(permM,chrom=chromosome,sample=1,kernel,kernelparams,step.width=step.width)
      return(permS[,2])
    }#computePermSlideScores
    if (verbose)
      cat(paste("Assessing",nperms,"permutations:\n"))
    permSlideScores <- apply(permMatrix,1,computePermSlideScores, verbose=verbose)
    if (verbose)
      cat("Computing expected borders...\n")
    expected.borders <- apply(permSlideScores,1,quantile,probs=c(0.025,0.975))
    lowerPermSlideScores <- expected.borders[1,]
    upperPermSlideScores <- expected.borders[2,]
  } # then (permute=="locations")
  
  if (permute=="labels"){
    if (verbose)
      cat("Compute sliding values for permutations...\n")
    lowerPermScores <- matrix(samResult$expected.lower,ncol=1)
    rownames(lowerPermScores) <- rownames(X)
    onChromPermScores <- lowerPermScores[onChromGenes,,drop=FALSE]
    lowerPermSlideScores <- kernelize(onChromPermScores,kernelweights)
    upperPermScores <- matrix(samResult$expected.upper,ncol=1)
    rownames(upperPermScores) <- rownames(X)
    onChromPermScores <- upperPermScores[onChromGenes,,drop=FALSE]
    upperPermSlideScores <- kernelize(onChromPermScores,kernelweights)
  } # then (permute=="labels")
  
  sortedLocs <- sort(absOnChromLocs,method="quick",index.return=TRUE)
  sortedScores <- onChromScores[sortedLocs$ix]
  sortedPvalues <- onChromPvalues[sortedLocs$ix]
  sortedGenes <- onChromGenes[sortedLocs$ix]
  result <- list(original.geneid=sortedGenes,original.loc=sortedLocs$x,
                 original.score=sortedScores,original.pvalue=sortedPvalues,
                 steps=steps,sliding.value=classSlideScores,
                 lower.permuted.border=lowerPermSlideScores,
                 upper.permuted.border=upperPermSlideScores,
                 chromosome=chromosome,class=class,chip=data$chip)
  class(result) <- "MACATevalScoring"
  if (verbose) cat("All done.\n")
  return(result)
} #evalScoring
  
plot.MACATevalScoring <- function(x, output="x11", HTMLfilename=NULL, mytitle=NULL ,new.device=TRUE,...){
  # x : result from MACAT-function 'evalScoring'
  # output: one of "x11","html"
  # HTMLfilename: filname for HTMLpage, default: 'results.html'
  # mytitle: tiltle of HTMLpage, default: 'Results'
  #attach(x); on.exit(detach(x))
  this.call <- as.list(match.call())
  if (!interactive())
    output <- "none"
  output <- match.arg(output, c("x11","html", "none"))
  require(gsub("\\.db\\.db$",".db",paste(x$chip,"db",sep=".")),character.only=TRUE)
  chromlength <- eval(as.symbol(paste(gsub("\\.db$","",x$chip),"CHRLENGTHS",sep="")))[x$chromosome]
  lowestpos <- 0 # old: min(min(original.loc),min(steps))
  highestpos <- chromlength # old: max(max(original.loc),max(steps))
  lowestscore <- min(min(x$original.score),min(x$sliding.value),min(x$lower.permuted.border))
  highestscore <- max(max(x$original.score),max(x$sliding.value),max(x$upper.permuted.border))
  if ((new.device) & interactive()){
    if (output=="x11"){
      if (capabilities()["X11"])
        x11(width=12,height=6)
    } else if (output=="html"){
      if (is.null(HTMLfilename))
        HTMLfilename=paste("Results",x$chromosome, "_" ,x$class ,".html",sep="")
      if (is.null(mytitle))
        mytitle=paste("Results for class ", x$class, " on chromosome ", x$chromosome, sep="")
      Slidingpic=paste("scoreplot", x$chromosome,"_",x$class, ".png", sep="")
      png(Slidingpic,width=900,height=480) ####
    }
    if (output=="x11")
      par(mar=c(5,5,4,1))
    else
      par(mar=c(1,5,4,1))
  } # if (new.device)

  if (!("xlim" %in% names(this.call)) || output=="html")
    plot(lowestpos,lowestscore,type="n",xaxt="n",xlab=NA,ylab="Score",
         xlim= c(lowestpos,highestpos), ylim=c(lowestscore,highestscore),
         frame.plot=FALSE,...)
  else
    plot(lowestpos,lowestscore,type="n",xaxt="n",xlab=NA,ylab="Score",
         ylim=c(lowestscore,highestscore), frame.plot=FALSE,...)

  points(x$original.loc,x$original.score,pch=20,cex=0.7,col="black")
  lines(x$steps,x$sliding.value,col="red",lwd=3)
  issig <- ((x$sliding.value>x$upper.permuted.border)|(x$sliding.value<x$lower.permuted.border))
  points(x$steps[issig],x$sliding.value[issig],col="gold",pch=20,cex=0.7,lwd=2)
  lines(x$steps,x$lower.permuted.border,col="gray",lwd=3)
  lines(x$steps,x$upper.permuted.border,col="gray",lwd=3)
  title(main=paste("Class",x$class,", Scores for Chromosome",x$chromosome))
  if ((output=="x11")&(new.device)){
    axis(1)
    mtext("Coordinate",1,line=3)
  } # if ((output=="x11")&(new.device))
  if (output=="html"){
    dev.off()
    ####
    getHtml(x, Slidingpic, HTMLfilename, mytitle)
  }
} # plot.MACATevalScoring