R/utils.R

In benilton/crlmmOld: Genotype Calling (CRLMM) and Copy Number Analysis tool for Affymetrix SNP 5.0 and 6.0 and Illumina arrays.

stupidSplineBasis <- function(x,knots){
  x <- pmin(x,knots[3])
  x <- pmax(x,knots[1])
  cbind(1, x, x^2, pmax(0, (x-knots[2]))^2)
}

changeToCrlmmAnnotationName <- function(x){
  pkgDir <- system.file(package=THISPKG)
  wanted <- paste(tolower(gsub("_", "", x)), "crlmm.stuff", sep=".")
  file.path(pkgDir, "extdata", wanted)
}

getCrlmmAnnotationName <- function(x){
	paste(tolower(gsub("_", "", x)), "Crlmm", sep="")
}

## medianSummaries <- function(mat, grps)
##   .Call("R_subColSummarize_median", mat, grps, PACKAGE = "preprocessCore")

medianSummaries <- function(mat, grps)
.Call("subColSummarizeMedianPP", mat, grps)

intMedianSummaries <- function(mat, grps)
  as.integer(medianSummaries(mat, grps))



## .crlmmPkgEnv is an enviroment that will
## store all the variables used by the pkg.
## it's meant to not overwrite user's variables
## and get rid of the NOTES generated by
## R CMD check

isLoaded <- function(dataset, environ=.crlmmPkgEnv)
	exists(dataset, envir=environ)

getVarInEnv <- function(dataset, environ=.crlmmPkgEnv){
	if (!isLoaded(dataset, environ=environ))
		stop("Variable ", dataset, " not found in supplied environment")
	environ[[dataset]]
}

list2SnpSet <- function(x, returnParams=FALSE){
  pd <- data.frame(SNR=x[["SNR"]][], gender=x[["gender"]],
                   batchQC=rep(x[["batchQC"]], ncol(x[["calls"]])),
                   row.names=colnames(x[["calls"]]))
  pdv <- data.frame(labelDescription=c("Signal-to-noise Ratio",
                      "Gender: Male (1) and Female (2)",
                      "Quality score for batch"),
                    row.names=c("SNR", "gender", "batchQC"))

  recall <- length(x[["DD"]]) > 1
  if (returnParams){
    if (recall){
      fd <- data.frame(SNPQC=x[["SNPQC"]],
                       cAA=x[["params"]][["centers"]][,1],
                       cAB=x[["params"]][["centers"]][,2],
                       cBB=x[["params"]][["centers"]][,3],
                       sAA=x[["params"]][["scales"]][,1],
                       sAB=x[["params"]][["scales"]][,2],
                       sBB=x[["params"]][["scales"]][,3],
                       nAA=x[["params"]][["N"]][,1],
                       nAB=x[["params"]][["N"]][,2],
                       nBB=x[["params"]][["N"]][,3],
                       spAA=x[["DD"]][,1],
                       spAB=x[["DD"]][,2],
                       spBB=x[["DD"]][,3],
                       row.names=rownames(x[["calls"]]))
      fdv <- data.frame(labelDescription=c(
                          "SNP Quality Score",
                          "Center AA", "Center AB", "Center BB",
                          "Scale AA", "Scale AB", "Scale BB",
                          "N AA", "N AB", "N BB",
                          "Shift in parameters AA",
                          "Shift in parameters AB",
                          "Shift in parameters BB"),
                        row.names=c(
                          "SNPQC",
                          "cAA", "cAB", "cBB",
                          "sAA", "sAB", "sBB",
                          "nAA", "nAB", "nBB",
                          "spAA", "spAB", "spBB"))
    }else{
      fd <- data.frame(SNPQC=x[["SNPQC"]],
                       cAA=x[["params"]][["centers"]][,1],
                       cAB=x[["params"]][["centers"]][,2],
                       cBB=x[["params"]][["centers"]][,3],
                       sAA=x[["params"]][["scales"]][,1],
                       sAB=x[["params"]][["scales"]][,2],
                       sBB=x[["params"]][["scales"]][,3],
                       nAA=x[["params"]][["N"]][,1],
                       nAB=x[["params"]][["N"]][,2],
                       nBB=x[["params"]][["N"]][,3],
                       row.names=rownames(x[["calls"]]))
      fdv <- data.frame(labelDescription=c(
                          "SNP Quality Score",
                          "Center AA", "Center AB", "Center BB",
                          "Scale AA", "Scale AB", "Scale BB",
                          "N AA", "N AB", "N BB"),
                        row.names=c(
                          "SNPQC",
                          "cAA", "cAB", "cBB",
                          "sAA", "sAB", "sBB",
                          "nAA", "nAB", "nBB"))
    }
  }else{
    if (recall){
      fd <- data.frame(SNPQC=x[["SNPQC"]],
                       spAA=x[["DD"]][,1],
                       spAB=x[["DD"]][,2],
                       spBB=x[["DD"]][,3],
                       row.names=rownames(x[["calls"]]))
      fdv <- data.frame(labelDescription=c("SNP Quality Score",
                          "Shift in parameters AA",
                          "Shift in parameters AB",
                          "Shift in parameters BB"),
                        row.names=c("SNPQC", "spAA", "spAB", "spBB"))
    }else{
      fd <- data.frame(SNPQC=x[["SNPQC"]],
                       row.names=rownames(x[["calls"]]))
      fdv <- data.frame(labelDescription=c("SNP Quality Score"),
                        row.names=c("SNPQC"))
    }
  }
  new("SnpSet",
      assayData=assayDataNew("lockedEnvironment",
        call=x[["calls"]], callProbability=x[["confs"]]),
      phenoData=new("AnnotatedDataFrame",
        data=pd, varMetadata=pdv),
      featureData=new("AnnotatedDataFrame",
        data=fd, varMetadata=fdv),
      annotation=x[["pkgname"]])
}

loader <- function(theFile, envir, pkgname){
	theFile <- file.path(system.file(package=pkgname),
			     "extdata", theFile)
	if (!file.exists(theFile))
		stop("File ", theFile, " does not exist in ", pkgname)
	load(theFile, envir=envir)
}

celDates <- function(celfiles){
	if(!all(file.exists(celfiles))) stop("1 or more cel file does not exist")
	celdates <- celtimes <- vector("character", length(celfiles))
	for(i in seq(along=celfiles)){
		if(i %% 100 == 0) cat(".")
		tmp <- read.celfile.header(celfiles[i], info="full")$DatHeader
		tmp <- strsplit(tmp, "\ +")
		celdates[i] <- tmp[[1]][6]
		celtimes[i] <- tmp[[1]][7]
	}
	tmp <- paste(celdates, celtimes)
	celdts <- strptime(tmp, "%m/%d/%y %H:%M:%S")
	return(celdts)
}

illuminaCdfNames <- function(){
	c("human1mv1c",# 1M
	  "human370v1c",            # 370CNV
	  "human650v3a",            # 650Y
	  "human610quadv1b",        # 610 quad
	  "human660quadv1a",        # 660 quad
	  "human370quadv3c",        # 370CNV quad
	  "human550v3b",            # 550K
	  "human1mduov3b",          # 1M Duo
	  "humanomni1quadv1b",      # Omni1 quad
	  "humanomni25quadv1b",     # Omni2.5 quad
	  "humanomni258v1a",        # Omni2.5 8 v1 A
          "humanomni258v1p1b",      # Omni2.5 8 v1.1 B
          "humanomni5quadv1b",      # Omni5 quad          
	  "humanomniexpress12v1b",  # Omni express 12
	  "humanimmuno12v1b",       # Immuno chip 12
	  "humancytosnp12v2p1h",    # CytoSNP 12
          "humanexome12v1p2a",      # Exome 12 v1.2 A
          "humanomniexpexome8v1p1b")
}

affyCdfNames <- function(){
	c("genomewidesnp6",
	  "genomewidesnp5")
}

validCdfNames <- function(){
	c(affyCdfNames(),
	  illuminaCdfNames())
}

cleancdfname <- function(x) strsplit(x, "Crlmm")[[1]][[1]]

isValidCdfName <- function(cdfName){
	chipList <- validCdfNames()
	match.arg(cleancdfname(cdfName), chipList)
	return(TRUE)
}

isPackageLoaded <- function(pkg){
	stopifnot(is.character(pkg))
	pkg <- paste("package:", pkg, sep="")
	pkg %in% search()
}

paramNames <- function(){
	c("tau2A",
	  "tau2B", "sig2A", "sig2B",
	  "nuA", ##"nuA.se",
	  "nuB",  ##"nuB.se",
	  "phiA",
	  "phiB",
	  "phiPrimeA",
	  "phiPrimeB",
	  ##"phiA.se", "phiB.se",
	  "corrAB",
	  "corrBB",
	  "corrAA")
}

loadObject <- function(filename, load.it){
	fname <- paste(filename, ".rda", sep="")
	if(load.it & file.exists(file.path(ldPath(), fname))){
		message("load.it is TRUE, loading previously saved ff object")
		return(TRUE)
	} else return(FALSE)
}



setMethod("annotatedDataFrameFrom", "ff_matrix", Biobase:::annotatedDataFrameFromMatrix)
setMethod("annotatedDataFrameFrom", "ffdf", Biobase:::annotatedDataFrameFromMatrix)

## Document this...
getBAF <- function(theta, canonicalTheta)
    .Call('normalizeBAF', theta, canonicalTheta)


validCEL <- function(celfiles){
	for(i in seq_along(celfiles)){
		res <- tryCatch(read.celfile(celfiles[i], intensity.means.only=TRUE), error=function(e) NULL)
		if(is.null(res)) {
			msg <- message("Problem reading ", celfiles[i])
			stop(msg)
		}
	}
	return("Successfully read all cel files")
}