R/consensus.R

In seqinr: Biological Sequences Retrieval and Analysis

consensus <- function(matali, method = c( "majority", "threshold", "IUPAC", "profile"), threshold = 0.60,
  warn.non.IUPAC = FALSE, type = c("DNA", "RNA")){

  if(inherits(matali, "alignment")) matali <- as.matrix(matali)
  if(!is.matrix(matali)) stop("matrix or alignment object expected")
  if(storage.mode(matali) != "character") stop("matrix of characters expected")
  
  method <- match.arg(method)
  
  if(method == "IUPAC"){
    type <- match.arg(type)
    res <- apply(matali, 2, bma, warn.non.IUPAC = warn.non.IUPAC, type = type)
    names(res) <- NULL
    return(res)
  }
  
  if(method == "majority"){
    majority <- function(x) names(which.max(table(x)))
    res <- apply(matali, 2, majority)
    names(res) <- NULL
    return(res)
  }
  
  if(method == "profile"){
    obsvalue <- levels(factor(matali))
    nrow <- length(obsvalue)
    row.names(matali)<-NULL
    res <- apply(matali, 2, function(x) table(factor(x, levels = obsvalue)))
    return(res)
  }
  
  if(method == "threshold"){
    profile <- consensus(matali, method = "profile")
    profile.rf <- apply(profile, 2, function(x) x/sum(x))
    res <- rownames(profile.rf)[apply(profile.rf, 2, which.max)]
    res <- ifelse(apply(profile.rf, 2, max) >= threshold, res, NA)
    names(res) <- NULL
    return(res)
  }
}

con <- consensus