R/PoissonDistance.R

In PoiClaClu: Classification and Clustering of Sequencing Data Based on a Poisson Model

PoissonDistance <-
function(x,beta=1,type=c("mle","deseq","quantile"),transform=TRUE, alpha=NULL,perfeature=FALSE){
  type <- match.arg(type)
  if(!transform && !is.null(alpha)) stop("You have asked for NO transformation but have entered alpha.")
  if(transform && !is.null(alpha)){
    if(alpha>0 && alpha <= 1) x <- x^alpha
    if(alpha<=0 || alpha>1) stop("alpha must be between 0 and 1")
  }  
  if(transform && is.null(alpha)){
    alpha <- FindBestTransform(x)
    x <- x^alpha
  }
  dd <- matrix(0, nrow=nrow(x), ncol=nrow(x))
  ddd <- NULL
  if(perfeature) ddd <- array(0, dim=c(nrow(x), nrow(x), ncol(x)))
  for(i in 2:nrow(dd)){
    xi <- x[i,]
    for(j in 1:(i-1)){
      xj <- x[j,]
      n <- NullModel(x[c(i,j),],type=type)$n
      ni <- n[1,]
      nj <- n[2,]
      di <- (xi+beta)/(ni+beta)
      dj <- (xj+beta)/(nj+beta)
      dd[i,j] <- sum(ni+nj-ni*di-nj*dj+xi*log(di)+xj*log(dj))
      if(perfeature) ddd[j,i,] <- ddd[i,j,] <- ni+nj-ni*di-nj*dj+xi*log(di)+xj*log(dj)
    }
  }
  save <- list(dd=as.dist(dd+t(dd)), alpha=alpha, x=x, ddd=ddd, alpha=alpha, type=type)
  class(save) <- "poidist"
  return(save)
}

print.poidist <- function(x,...){
  if(!is.null(x$alpha)) cat("Value of alpha used to transform data: ", x$alpha, fill=TRUE)
  if(is.null(x$alpha)) cat("No transformation performed.",fill=TRUE)
  cat("This type of normalization was performed:", x$type, fill=TRUE)
  cat("Dissimilarity computed for ", nrow(x$x), " observations.", fill=TRUE)
}