@@ -308,3 +308,7 @@ is decoded and scanned

  * fixed bug in readIdatFiles()

  * updated crlmmIllumina.Rnw and crlmmIllumina.pdf

+

+2009-11-14 B Carvalho - committed version 1.5.2

+

+ * code cleanup - removed unneeded C code

@@ -71,7 +71,6 @@ crlmmGT <- function(A, B, SNR, mixtureParams, cdfName, row.names=NULL,

   if(verbose) message("Loading annotations.")

   loader("genotypeStuff.rda", .crlmmPkgEnv, pkgname)

   loader("mixtureStuff.rda", .crlmmPkgEnv, pkgname)

-##  data(genotypeStuff, mixtureStuff, package=pkgname, envir=.crlmmPkgEnv)

   ## this is toget rid of the 'no visible binding' notes

   ## variable definitions

@@ -85,13 +84,13 @@ crlmmGT <- function(A, B, SNR, mixtureParams, cdfName, row.names=NULL,

   ##IF gender not provide, we predict

   if(is.null(gender)){

-	  if(verbose) message("Determining gender.")

-	  XMedian <- apply(log2(A[XIndex,, drop=FALSE])+log2(B[XIndex,, drop=FALSE]), 2, median)/2

-	  if(sum(SNR>SNRMin)==1){

-		  gender <- which.min(c(abs(XMedian-8.9), abs(XMedian-9.5)))

-	  }else{

-		  gender <- kmeans(XMedian, c(min(XMedian[SNR>SNRMin]), max(XMedian[SNR>SNRMin])))[["cluster"]]

-	  }

+    if(verbose) message("Determining gender.")

+    XMedian <- apply(log2(A[XIndex,, drop=FALSE])+log2(B[XIndex,, drop=FALSE]), 2, median)/2

+    if(sum(SNR>SNRMin)==1){

+      gender <- which.min(c(abs(XMedian-8.9), abs(XMedian-9.5)))

+    }else{

+      gender <- kmeans(XMedian, c(min(XMedian[SNR>SNRMin]), max(XMedian[SNR>SNRMin])))[["cluster"]]

+    }

   }

   Indexes <- list(autosomeIndex, XIndex, YIndex)

@@ -206,15 +205,7 @@ crlmmGT <- function(A, B, SNR, mixtureParams, cdfName, row.names=NULL,

    tmp4batchQC <- DD[autosomeIndex,]/(params[["N"]][autosomeIndex,]+1)

    tmpSnpQc <- dev[autosomeIndex]

    SS <- cov(tmp4batchQC[tmpSnpQc < badSNP,])

-####     tmp4batchQC <- DD[snps2keep,]/(params[["N"]][snps2keep,]+1)

-####     tmpSnpQc <- dev[snps2keep]

-####     SS <- cov(DD[snps2keep,])

-

-##     DD <- sweep(GG[snps2keep, ], 2, colMeans(DD[snps2keep, ]))

-##     tmpSnpQc <- dev[snps2keep]

-##     SS <- cov(DD[tmpSnpQc < badSNP, ])

-####    batchQC <- sqrt(sum(diag(cov(tmp4batchQC[tmpSnpQc < badSNP,]))))*sum(params[["N"]][snps2keep[1],])

-    batchQC <- mean(diag(SS))

+   batchQC <- mean(diag(SS))

   }else{

     SS <- matrix(0, 3, 3)

     batchQC <- Inf

@@ -271,663 +262,3 @@ gtypeCallerR2 <- function(A, B, fIndex, mIndex, theCenters, theScales,

         as.integer(noTraining), as.integer(noInfo), PACKAGE="crlmm")

 }

-

-

-

-###############################

-####### THIS IS TEMPORARY NOT OFFICIALLY USED

-#####################################

-crlmmNM <- function(filenames, row.names=TRUE, col.names=TRUE,

-                  probs=c(1/3, 1/3, 1/3), DF=6, SNRMin=5, gender=NULL,

-                  save.it=FALSE, load.it=FALSE,

-                  intensityFile="tmpcrlmmintensities.rda",

-                  desctrucitve=FALSE, mixtureSampleSize=10^5, eps=0.1,

-                  verbose=TRUE, cdfName, sns, recallMin=10,

-                  recallRegMin=1000, returnParams=FALSE){

-  

-  if (missing(sns)) sns <- basename(filenames)

-  if (load.it & !file.exists(intensityFile)){

-    load.it <- FALSE

-    message("File ", intensityFile, " does not exist.")

-    message("Not loading it, but running SNPRMA from scratch.")

-  }

-  if (!load.it){

-    res <- snprma(filenames, fitMixture=TRUE,

-                  mixtureSampleSize=mixtureSampleSize, verbose=verbose,

-                  eps=eps, cdfName=cdfName, sns=sns)

-    if(save.it){

-      t0 <- proc.time()

-      save(res, file=intensityFile)

-      t0 <- proc.time()-t0

-      message("Used ", t0[3], " seconds to save ", intensityFile, ".")

-    }

-  }else{

-    message("Loading ", intensityFile, ".")

-    obj <- load(intensityFile)

-    message("Done.")

-    if (obj != "res")

-      stop("Object in ", intensityFile, " seems to be invalid.")

-  }

-  if(row.names) row.names=res$gns else row.names=NULL

-  if(col.names) col.names=res$sns else col.names=NULL

-

-  res2 <- crlmmGTnm(res[["A"]], res[["B"]], res[["SNR"]],

-                  res[["mixtureParams"]], res[["cdfName"]],

-                  gender=gender, row.names=row.names,

-                  col.names=col.names, recallMin=recallMin,

-                  recallRegMin=1000, SNRMin=SNRMin,

-                  returnParams=returnParams)

-

-  res2[["SNR"]] <- res[["SNR"]]

-  

-  return(res2)

-}

-

-crlmmTNoN <- function(filenames, row.names=TRUE, col.names=TRUE,

-                  probs=c(1/3, 1/3, 1/3), DF=6, SNRMin=5, gender=NULL,

-                  save.it=FALSE, load.it=FALSE,

-                  intensityFile="tmpcrlmmintensities.rda",

-                  desctrucitve=FALSE, mixtureSampleSize=10^5, eps=0.1,

-                  verbose=TRUE){

-  if (load.it & !file.exists(intensityFile)){

-    load.it <- FALSE

-    message("File ", intensityFile, " does not exist.")

-    message("Not loading it, but running SNPRMA from scratch.")

-  }

-  if (!load.it){

-    res <- snprma(filenames, fitMixture=TRUE,

-                  mixtureSampleSize=mixtureSampleSize, verbose=verbose,

-                  eps=eps)

-    if(save.it) save(res, file=intensityFile)

-  }else{

-    message("Loading ", intensityFile, ".")

-    obj <- load(intensityFile)

-    message("Done.")

-    if (obj != "res")

-      stop("Object in ", intensityFile, " seems to be invalid.")

-  }

-  if(row.names) row.names=res$gns else row.names=NULL

-  if(col.names) col.names=res$sns else col.names=NULL

-  res2 <- crlmmGTTNoN(res[["A"]], res[["B"]], res[["SNR"]],

-                  res[["mixtureParams"]], res[["cdfName"]],

-                  gender=gender, row.names=row.names,

-                  col.names=col.names)

-  res2[["SNR"]] <- res[["SNR"]]

-  return(res2)

-}

-

-crlmmNormalNoN <- function(filenames, row.names=TRUE, col.names=TRUE,

-                  probs=c(1/3, 1/3, 1/3), DF=6, SNRMin=5, gender=NULL,

-                  save.it=FALSE, load.it=FALSE,

-                  intensityFile="tmpcrlmmintensities.rda",

-                  desctrucitve=FALSE, mixtureSampleSize=10^5, eps=0.1,

-                  verbose=TRUE){

-  if (load.it & !file.exists(intensityFile)){

-    load.it <- FALSE

-    message("File ", intensityFile, " does not exist.")

-    message("Not loading it, but running SNPRMA from scratch.")

-  }

-  if (!load.it){

-    res <- snprma(filenames, fitMixture=TRUE,

-                  mixtureSampleSize=mixtureSampleSize, verbose=verbose,

-                  eps=eps)

-    if(save.it) save(res, file=intensityFile)

-  }else{

-    message("Loading ", intensityFile, ".")

-    obj <- load(intensityFile)

-    message("Done.")

-    if (obj != "res")

-      stop("Object in ", intensityFile, " seems to be invalid.")

-  }

-  if(row.names) row.names=res$gns else row.names=NULL

-  if(col.names) col.names=res$sns else col.names=NULL

-  res2 <- crlmmGTNormalNoN(res[["A"]], res[["B"]], res[["SNR"]],

-                  res[["mixtureParams"]], res[["cdfName"]],

-                  gender=gender, row.names=row.names,

-                  col.names=col.names)

-  res2[["SNR"]] <- res[["SNR"]]

-  return(res2)

-}

-

-crlmmGTTNoN <- function(A, B, SNR, mixtureParams, cdfName,

-                         row.names=NULL, col.names=NULL, probs=c(1/3,

-                         1/3, 1/3), DF=6, SNRMin=6, gender=NULL,

-                         desctrucitve=FALSE, verbose=TRUE){

-  keepIndex <- which(SNR>SNRMin)

-  if(length(keepIndex)==0) stop("No arrays above quality threshold!")

-  

-  NC <- ncol(A)

-  NR <- nrow(A)

-  pkgname <- getCrlmmAnnotationName(cdfName)

-  

-  if(verbose) message("Loading annotations.")

-  loader("genotypeStuff.rda", .crlmmPkgEnv, pkgname)

-  loader("mixtureStuff.rda", .crlmmPkgEnv, pkgname)

-##  data(genotypeStuff, mixtureStuff, package=pkgname, envir=.crlmmPkgEnv)

-

-  ## this is toget rid of the 'no visible binding' notes

-  ## variable definitions

-  XIndex <- getVarInEnv("XIndex")

-  autosomeIndex <- getVarInEnv("autosomeIndex")

-  YIndex <- getVarInEnv("YIndex")

-  SMEDIAN <- getVarInEnv("SMEDIAN")

-  theKnots <- getVarInEnv("theKnots")

-  regionInfo <- getVarInEnv("regionInfo")

-  

-  ##IF gender not provide, we predict

-  if(is.null(gender)){

-    if(verbose) message("Determining gender.")

-    XMedian <- apply(log2(A[XIndex,, drop=FALSE])+log2(B[XIndex,, drop=FALSE]), 2, median)/2

-    if(sum(SNR>SNRMin)==1) gender <- which.min(c(abs(XMedian-8.9), abs(XMedian-9.5))) else  gender <- kmeans(XMedian, c(min(XMedian[SNR>SNRMin]), max(XMedian[SNR>SNRMin])))[["cluster"]]

-  }

-  

-  Indexes <- list(autosomeIndex, XIndex, YIndex)

-  cIndexes <- list(keepIndex, 

-                   keepIndex[which(gender[keepIndex]==2)], 

-                   keepIndex[which(gender[keepIndex]==1)])

-  

-  if(verbose) cat("Calling", NR, "SNPs for recalibration... ")

-

-  ## call C

-  fIndex <- which(gender==2)

-  mIndex <- which(gender==1)

-  t0 <- proc.time()

-  newparams <- gtypeCallerRTNoN(A, B, fIndex, mIndex,

-                            params[["centers"]], params[["scales"]], params[["N"]],

-                            Indexes, cIndexes,

-                            sapply(Indexes, length), sapply(cIndexes, length),

-                            SMEDIAN, theKnots,

-                            mixtureParams, DF, probs, 0.025)

-  t0 <- proc.time()-t0

-  message("Part 1 took ", t0[3], " seconds.")

-  names(newparams) <- c("centers", "scales", "N")

-  

-  if(verbose) message("Done.")

-  if(verbose) message("Estimating recalibration parameters.")

-  d <- newparams[["centers"]] - params$centers

-

-  ##regression 

-  MIN <- 10

-  Index <- intersect(which(pmin(newparams[["N"]][, 1], newparams[["N"]][, 2], newparams[["N"]][, 3])>MIN & !apply(regionInfo, 1, any)), autosomeIndex)

-  data4reg <- as.data.frame(newparams[["centers"]][Index,])

-  names(data4reg) <- c("AA", "AB", "BB")

-  regParams <- cbind(  coef(lm(AA~AB*BB, data=data4reg)),

-                     c(coef(lm(AB~AA+BB, data=data4reg)), 0), 

-                       coef(lm(BB~AA*AB, data=data4reg)))

-  rownames(regParams) <- c("intercept", "X", "Y", "XY")

-  rm(data4reg)

-  

-  minN <- 3

-  newparams[["centers"]][newparams[["N"]]<minN] <- NA

-  Index <- setdiff(which(rowSums(is.na(newparams[["centers"]]))==1), YIndex)

-  if(verbose) cat("Filling out empty centers")

-  for(i in Index){

-    if(verbose) if(i%%10000==0)cat(".")

-    mu <- newparams[["centers"]][i, ]

-    j <- which(is.na(mu))

-    newparams[["centers"]][i, j] <- c(1, mu[-j], prod(mu[-j]))%*%regParams[, j]

-  }

-  ##remaing NAs are made like originals

-  if(length(YIndex)>0){

-    noMoveIndex <- union(setdiff(which(rowSums(is.na(newparams[["centers"]]))>0), YIndex), 

-                         YIndex[rowSums(is.na(newparams[["centers"]][YIndex, ])>1)])

-  }

-  newparams[["centers"]][is.na(newparams[["centers"]])] <- params[["centers"]][is.na(newparams[["centers"]])]

-  if(verbose) cat("\n")

-  

-  if(verbose) message("Calculating and standardizing size of shift.")

-  DD <- newparams[["centers"]] - params[["centers"]]

-  MM <- colMeans(DD[autosomeIndex, ])

-  DD <- sweep(DD, 2, MM)

-  SS <- cov(DD[autosomeIndex, ])

-  SSI <- solve(SS)

-  dev <- vector("numeric", nrow(DD))

-  if(length(YIndex)){

-    dev[-YIndex] <- apply(DD[-YIndex, ], 1, function(x) x%*%SSI%*%x)

-    dev[-YIndex] <- 1/sqrt( (2*pi)^3*det(SS))*exp(-0.5*dev[-YIndex])

-    ##Now Y (only two params)

-    SSY <- SS[c(1, 3), c(1, 3)]

-    SSI <- solve(SSY) 

-    dev[YIndex] <- apply(DD[YIndex, c(1, 3)], 1, function(x) x%*%SSI%*%x)

-    dev[YIndex] <- 1/sqrt( (2*pi)^2*det(SSY))*exp(-0.5*dev[YIndex])

-  } else {

-    dev=apply(DD,1,function(x) x%*%SSI%*%x)

-    dev=1/sqrt( (2*pi)^3*det(SS))*exp(-0.5*dev)

-  } 

-

-  ## BC: must keep SD

-  params[-2] <- newparams[-2]

-  rm(newparams);gc(verbose=FALSE)  

-  if(verbose) cat("Calling", NR, "SNPs... ")

-  ## ###################

-  ## ## MOVE TO C#######

-  t0 <- proc.time()

-  ImNull <- gtypeCallerR2TNoN(A, B, fIndex, mIndex, params[["centers"]],

-                          params[["scales"]], params[["N"]], Indexes,

-                          cIndexes, sapply(Indexes, length),

-                          sapply(cIndexes, length), SMEDIAN, theKnots,

-                          mixtureParams, DF, probs, 0.025,

-                          which(regionInfo[,2]),

-                          which(regionInfo[,1]))

-  t0 <- proc.time()-t0

-  message("Part 2 took ", t0[3], " seconds.")

-  ##  END MOVE TO C#######

-  ## ##################

-  

-  dev <- dev/(dev+1/383)

-  if(!is.null(row.names)){ rownames(A) <- rownames(B) <- names(dev) <- row.names}

-  if(!is.null(col.names)){ colnames(A) <- colnames(B) <- col.names}

-  

-  if(verbose) message("Done.")

-  return(list(calls=A, confs=B, SNPQC=dev, batchQC=mean(diag(SS))))

-}

-

-

-gtypeCallerRTNoN <- function(A, B, fIndex, mIndex, theCenters, theScales,

-                         theNs, Indexes, cIndexes, nIndexes,

-                         ncIndexes, SMEDIAN, knots, params, dft,

-                         probs, trim){

-

-  stopifnot(!missing(A), !missing(B), dim(A)==dim(B),

-            nrow(A)==nrow(theCenters), nrow(A)==nrow(theScales),

-            nrow(A) == nrow(theNs), length(knots)==3, nrow(params)==4,

-            ncol(params)==ncol(A), length(trim)==1, length(probs)==3)

-  

-  .Call("gtypeCallerPart1TNoN", A, B, fIndex, mIndex, theCenters,

-        theScales, theNs, Indexes, cIndexes, nIndexes, ncIndexes,

-        SMEDIAN, knots, params, dft, probs, trim, PACKAGE="crlmm")

-  

-}

-

-gtypeCallerR2TNoN <- function(A, B, fIndex, mIndex, theCenters, theScales,

-                         theNs, Indexes, cIndexes, nIndexes,

-                         ncIndexes, SMEDIAN, knots, params, dft,

-                         probs, trim, noTraining, noInfo){

-

-  stopifnot(!missing(A), !missing(B), dim(A)==dim(B),

-            nrow(A)==nrow(theCenters), nrow(A)==nrow(theScales),

-            nrow(A) == nrow(theNs), length(knots)==3, nrow(params)==4,

-            ncol(params)==ncol(A), length(trim)==1, length(probs)==3)

-  

-  .Call("gtypeCallerPart2TNoN", A, B, fIndex, mIndex, theCenters,

-        theScales, theNs, Indexes, cIndexes, nIndexes, ncIndexes,

-        SMEDIAN, knots, params, dft, probs, trim, noTraining, noInfo, PACKAGE="crlmm")

-  

-}

-

-crlmmGTNormalNoN <- function(A, B, SNR, mixtureParams, cdfName,

-                         row.names=NULL, col.names=NULL, probs=c(1/3,

-                         1/3, 1/3), DF=6, SNRMin=6, gender=NULL,

-                         desctrucitve=FALSE, verbose=TRUE){

-  keepIndex <- which(SNR>SNRMin)

-  if(length(keepIndex)==0) stop("No arrays above quality threshold!")

-  

-  NC <- ncol(A)

-  NR <- nrow(A)

-  pkgname <- getCrlmmAnnotationName(cdfName)

-

-  if(verbose) message("Loading annotations.")

-  loader("genotypeStuff.rda", .crlmmPkgEnv, pkgname)

-  loader("mixtureStuff.rda", .crlmmPkgEnv, pkgname)

-##  data(genotypeStuff, mixtureStuff, package=pkgname, envir=.crlmmPkgEnv)

-

-  ## this is toget rid of the 'no visible binding' notes

-  ## variable definitions

-  XIndex <- getVarInEnv("XIndex")

-  autosomeIndex <- getVarInEnv("autosomeIndex")

-  YIndex <- getVarInEnv("YIndex")

-  SMEDIAN <- getVarInEnv("SMEDIAN")

-  theKnots <- getVarInEnv("theKnots")

-  regionInfo <- getVarInEnv("regionInfo")

-  

-  ##IF gender not provide, we predict

-  if(is.null(gender)){

-    if(verbose) message("Determining gender.")

-    XMedian <- apply(log2(A[XIndex,, drop=FALSE])+log2(B[XIndex,, drop=FALSE]), 2, median)/2

-    if(sum(SNR>SNRMin)==1) gender <- which.min(c(abs(XMedian-8.9), abs(XMedian-9.5))) else  gender <- kmeans(XMedian, c(min(XMedian[SNR>SNRMin]), max(XMedian[SNR>SNRMin])))[["cluster"]]

-  }

-  

-  Indexes <- list(autosomeIndex, XIndex, YIndex)

-  cIndexes <- list(keepIndex, 

-                   keepIndex[which(gender[keepIndex]==2)], 

-                   keepIndex[which(gender[keepIndex]==1)])

-  

-  if(verbose) cat("Calling", NR, "SNPs for recalibration... ")

-

-  ## call C

-  fIndex <- which(gender==2)

-  mIndex <- which(gender==1)

-  t0 <- proc.time()

-  newparams <- gtypeCallerRNormalNoN(A, B, fIndex, mIndex,

-                            params[["centers"]], params[["scales"]], params[["N"]],

-                            Indexes, cIndexes,

-                            sapply(Indexes, length), sapply(cIndexes, length),

-                            SMEDIAN, theKnots,

-                            mixtureParams, DF, probs, 0.025)

-  t0 <- proc.time()-t0

-  message("Part 1 took ", t0[3], " seconds.")

-  names(newparams) <- c("centers", "scales", "N")

-  

-  if(verbose) message("Done.")

-  if(verbose) message("Estimating recalibration parameters.")

-  d <- newparams[["centers"]] - params$centers

-

-  ##regression 

-  MIN <- 10

-  Index <- intersect(which(pmin(newparams[["N"]][, 1], newparams[["N"]][, 2], newparams[["N"]][, 3])>MIN & !apply(regionInfo, 1, any)), autosomeIndex)

-  data4reg <- as.data.frame(newparams[["centers"]][Index,])

-  names(data4reg) <- c("AA", "AB", "BB")

-  regParams <- cbind(  coef(lm(AA~AB*BB, data=data4reg)),

-                     c(coef(lm(AB~AA+BB, data=data4reg)), 0), 

-                       coef(lm(BB~AA*AB, data=data4reg)))

-  rownames(regParams) <- c("intercept", "X", "Y", "XY")

-  rm(data4reg)

-  

-  minN <- 3

-  newparams[["centers"]][newparams[["N"]]<minN] <- NA

-  Index <- setdiff(which(rowSums(is.na(newparams[["centers"]]))==1), YIndex)

-  if(verbose) cat("Filling out empty centers")

-  for(i in Index){

-    if(verbose) if(i%%10000==0)cat(".")

-    mu <- newparams[["centers"]][i, ]

-    j <- which(is.na(mu))

-    newparams[["centers"]][i, j] <- c(1, mu[-j], prod(mu[-j]))%*%regParams[, j]

-  }

-  ##remaing NAs are made like originals

-  if(length(YIndex)>0){

-    noMoveIndex <- union(setdiff(which(rowSums(is.na(newparams[["centers"]]))>0), YIndex), 

-                         YIndex[rowSums(is.na(newparams[["centers"]][YIndex, ])>1)])

-  }

-  newparams[["centers"]][is.na(newparams[["centers"]])] <- params[["centers"]][is.na(newparams[["centers"]])]

-  if(verbose) cat("\n")

-  

-  if(verbose) message("Calculating and standardizing size of shift.")

-  DD <- newparams[["centers"]] - params[["centers"]]

-  MM <- colMeans(DD[autosomeIndex, ])

-  DD <- sweep(DD, 2, MM)

-  SS <- cov(DD[autosomeIndex, ])

-  SSI <- solve(SS)

-  dev <- vector("numeric", nrow(DD))

-  if(length(YIndex)){

-    dev[-YIndex] <- apply(DD[-YIndex, ], 1, function(x) x%*%SSI%*%x)

-    dev[-YIndex] <- 1/sqrt( (2*pi)^3*det(SS))*exp(-0.5*dev[-YIndex])

-    ##Now Y (only two params)

-    SSY <- SS[c(1, 3), c(1, 3)]

-    SSI <- solve(SSY) 

-    dev[YIndex] <- apply(DD[YIndex, c(1, 3)], 1, function(x) x%*%SSI%*%x)

-    dev[YIndex] <- 1/sqrt( (2*pi)^2*det(SSY))*exp(-0.5*dev[YIndex])

-  } else {

-    dev=apply(DD,1,function(x) x%*%SSI%*%x)

-    dev=1/sqrt( (2*pi)^3*det(SS))*exp(-0.5*dev)

-  } 

-

-  ## BC: must keep SD

-  params[-2] <- newparams[-2]

-  rm(newparams);gc(verbose=FALSE)  

-  if(verbose) cat("Calling", NR, "SNPs... ")

-  ## ###################

-  ## ## MOVE TO C#######

-  t0 <- proc.time()

-  ImNull <- gtypeCallerR2NormalNoN(A, B, fIndex, mIndex, params[["centers"]],

-                          params[["scales"]], params[["N"]], Indexes,

-                          cIndexes, sapply(Indexes, length),

-                          sapply(cIndexes, length), SMEDIAN, theKnots,

-                          mixtureParams, DF, probs, 0.025,

-                          which(regionInfo[,2]),

-                          which(regionInfo[,1]))

-  t0 <- proc.time()-t0

-  message("Part 2 took ", t0[3], " seconds.")

-  ##  END MOVE TO C#######

-  ## ##################

-  

-  dev <- dev/(dev+1/383)

-  if(!is.null(row.names)){ rownames(A) <- rownames(B) <- names(dev) <- row.names}

-  if(!is.null(col.names)){ colnames(A) <- colnames(B) <- col.names}

-  

-  if(verbose) message("Done.")

-  return(list(calls=A, confs=B, SNPQC=dev, batchQC=mean(diag(SS))))

-}

-

-

-gtypeCallerRNormalNoN <- function(A, B, fIndex, mIndex, theCenters, theScales,

-                         theNs, Indexes, cIndexes, nIndexes,

-                         ncIndexes, SMEDIAN, knots, params, dft,

-                         probs, trim){

-

-  stopifnot(!missing(A), !missing(B), dim(A)==dim(B),

-            nrow(A)==nrow(theCenters), nrow(A)==nrow(theScales),

-            nrow(A) == nrow(theNs), length(knots)==3, nrow(params)==4,

-            ncol(params)==ncol(A), length(trim)==1, length(probs)==3)

-  

-  .Call("gtypeCallerPart1NormalNoN", A, B, fIndex, mIndex, theCenters,

-        theScales, theNs, Indexes, cIndexes, nIndexes, ncIndexes,

-        SMEDIAN, knots, params, dft, probs, trim, PACKAGE="crlmm")

-  

-}

-

-gtypeCallerR2NormalNoN <- function(A, B, fIndex, mIndex, theCenters, theScales,

-                         theNs, Indexes, cIndexes, nIndexes,

-                         ncIndexes, SMEDIAN, knots, params, dft,

-                         probs, trim, noTraining, noInfo){

-

-  stopifnot(!missing(A), !missing(B), dim(A)==dim(B),

-            nrow(A)==nrow(theCenters), nrow(A)==nrow(theScales),

-            nrow(A) == nrow(theNs), length(knots)==3, nrow(params)==4,

-            ncol(params)==ncol(A), length(trim)==1, length(probs)==3)

-  

-  .Call("gtypeCallerPart2NormalNoN", A, B, fIndex, mIndex, theCenters,

-        theScales, theNs, Indexes, cIndexes, nIndexes, ncIndexes,

-        SMEDIAN, knots, params, dft, probs, trim, noTraining, noInfo, PACKAGE="crlmm")

-  

-}

-

-

-crlmmGTnm <- function(A, B, SNR, mixtureParams, cdfName, row.names=NULL,

-                    col.names=NULL, probs=c(1/3, 1/3, 1/3), DF=6,

-                    SNRMin=5, recallMin=10, recallRegMin=1000,

-                    gender=NULL, desctrucitve=FALSE, verbose=TRUE,

-                    returnParams=FALSE){

-  

-  keepIndex <- which(SNR>SNRMin)

-  if(length(keepIndex)==0) stop("No arrays above quality threshold!")

-  

-  NC <- ncol(A)

-  NR <- nrow(A)

-  

-  pkgname <- getCrlmmAnnotationName(cdfName)

-  if(!require(pkgname, character.only=TRUE)){

-    suggCall <- paste("library(", pkgname, ", lib.loc='/Altern/Lib/Loc')", sep="")

-    msg <- paste("If", pkgname, "is installed on an alternative location, please load it manually by using", suggCall)

-    message(strwrap(msg))

-    stop("Package ", pkgname, " could not be found.")

-    rm(suggCall, msg)

-  }

-

-  if(verbose) message("Loading annotations.")

-  loader("genotypeStuff.rda", .crlmmPkgEnv, pkgname)

-  loader("mixtureStuff.rda", .crlmmPkgEnv, pkgname)

-##  data(genotypeStuff, mixtureStuff, package=pkgname, envir=.crlmmPkgEnv)

-

-  ## this is toget rid of the 'no visible binding' notes

-  ## variable definitions

-  XIndex <- getVarInEnv("XIndex")

-  autosomeIndex <- getVarInEnv("autosomeIndex")

-  YIndex <- getVarInEnv("YIndex")

-  SMEDIAN <- getVarInEnv("SMEDIAN")

-  theKnots <- getVarInEnv("theKnots")

-  regionInfo <- getVarInEnv("regionInfo")

-  

-  ##IF gender not provide, we predict

-  if(is.null(gender)){

-    if(verbose) message("Determining gender.")

-    XMedian <- apply(log2(A[XIndex,, drop=FALSE])+log2(B[XIndex,, drop=FALSE]), 2, median)/2

-    if(sum(SNR>SNRMin)==1){

-      gender <- which.min(c(abs(XMedian-8.9), abs(XMedian-9.5)))

-    }else{

-      gender <- kmeans(XMedian, c(min(XMedian[SNR>SNRMin]), max(XMedian[SNR>SNRMin])))[["cluster"]]

-    }

-  }

-  

-  Indexes <- list(autosomeIndex, XIndex, YIndex)

-  cIndexes <- list(keepIndex, 

-                   keepIndex[which(gender[keepIndex]==2)], 

-                   keepIndex[which(gender[keepIndex]==1)])

-  

-  if(verbose) cat("Calling", NR, "SNPs for recalibration... ")

-

-  ## call C

-  fIndex <- which(gender==2)

-  mIndex <- which(gender==1)

-  t0 <- proc.time()

-  newparams <- gtypeCallerRnm(A, B, fIndex, mIndex,

-                            params[["centers"]], params[["scales"]], params[["N"]],

-                            Indexes, cIndexes,

-                            sapply(Indexes, length), sapply(cIndexes, length),

-                            SMEDIAN, theKnots,

-                            mixtureParams, DF, probs, 0.025)

-  t0 <- proc.time()-t0

-  message("Part 1 took ", t0[3], " seconds.")

-  gc(verbose=FALSE)

-  names(newparams) <- c("centers", "scales", "N")

-  

-  if(verbose) message("Done.")

-  if(verbose) message("Estimating recalibration parameters.")

-  d <- newparams[["centers"]] - params$centers

-

-  ##regression 

-  Index <- intersect(which(pmin(newparams[["N"]][, 1],

-                                newparams[["N"]][, 2],

-                                newparams[["N"]][, 3]) > recallMin &

-                                !apply(regionInfo, 1, any)),

-                                autosomeIndex)

-  

-  if(length(Index) < recallRegMin){

-    warning("Recallibration not possible.")

-    newparams <- params

-    dev <- vector("numeric", nrow(newparams[["centers"]]))

-    SS <- matrix(Inf, 3, 3)

-  }else{

-    data4reg <- as.data.frame(newparams[["centers"]][Index,])

-    names(data4reg) <- c("AA", "AB", "BB")

-    regParams <- cbind(  coef(lm(AA~AB*BB, data=data4reg)),

-                       c(coef(lm(AB~AA+BB, data=data4reg)), 0), 

-                       coef(lm(BB~AA*AB, data=data4reg)))

-    rownames(regParams) <- c("intercept", "X", "Y", "XY")

-    rm(data4reg)

-  

-    minN <- 3

-    newparams[["centers"]][newparams[["N"]] < minN] <- NA

-    Index <- setdiff(which(rowSums(is.na(newparams[["centers"]]))==1), YIndex)

-    if(verbose) cat("Filling out empty centers")

-    for(i in Index){

-      if(verbose) if(i%%10000==0)cat(".")

-      mu <- newparams[["centers"]][i, ]

-      j <- which(is.na(mu))

-      newparams[["centers"]][i, j] <- c(1, mu[-j], prod(mu[-j]))%*%regParams[, j]

-    }

-    

-    ##remaing NAs are made like originals

-    if(length(YIndex)>0){

-      noMoveIndex <- union(setdiff(which(rowSums(is.na(newparams[["centers"]]))>0), YIndex), 

-                           YIndex[rowSums(is.na(newparams[["centers"]][YIndex, ])>1)])

-    }

-    newparams[["centers"]][is.na(newparams[["centers"]])] <- params[["centers"]][is.na(newparams[["centers"]])]

-    if(verbose) cat("\n")

-  

-    if(verbose) message("Calculating and standardizing size of shift.")

-    DD <- newparams[["centers"]] - params[["centers"]]

-    DD <- sweep(DD, 2, colMeans(DD[autosomeIndex, ]))

-    SS <- cov(DD[autosomeIndex, ])

-    SSI <- solve(SS)

-    dev <- vector("numeric", nrow(DD))

-    if(length(YIndex)){

-      dev[-YIndex] <- apply(DD[-YIndex, ], 1, function(x) x%*%SSI%*%x)

-      dev[-YIndex] <- 1/sqrt( (2*pi)^3*det(SS))*exp(-0.5*dev[-YIndex])

-      ##Now Y (only two params)

-      SSY <- SS[c(1, 3), c(1, 3)]

-      SSI <- solve(SSY) 

-      dev[YIndex] <- apply(DD[YIndex, c(1, 3)], 1, function(x) x%*%SSI%*%x)

-      dev[YIndex] <- 1/sqrt( (2*pi)^2*det(SSY))*exp(-0.5*dev[YIndex])

-    } else {

-      dev=apply(DD,1,function(x) x%*%SSI%*%x)

-      dev=1/sqrt( (2*pi)^3*det(SS))*exp(-0.5*dev)

-    }

-  }

-    

-  ## BC: must keep SD

-  params[-2] <- newparams[-2]

-  

-  rm(newparams);gc(verbose=FALSE)  

-  if(verbose) cat("Calling", NR, "SNPs... ")

-  ## ###################

-  ## ## MOVE TO C#######

-  t0 <- proc.time()

-  ImNull <- gtypeCallerR2nm(A, B, fIndex, mIndex, params[["centers"]],

-                          params[["scales"]], params[["N"]], Indexes,

-                          cIndexes, sapply(Indexes, length),

-                          sapply(cIndexes, length), SMEDIAN, theKnots,

-                          mixtureParams, DF, probs, 0.025,

-                          which(regionInfo[,2]),

-                          which(regionInfo[,1]))

-  t0 <- proc.time()-t0

-  gc(verbose=FALSE)

-  message("\n Part 2 took ", t0[3], " seconds.")

-  ##  END MOVE TO C#######

-  ## ##################

-  

-  dev <- dev/(dev+1/383)

-  if(!is.null(row.names)){ rownames(A) <- rownames(B) <- names(dev) <- row.names}

-  if(!is.null(col.names)){ colnames(A) <- colnames(B) <- col.names}

-  

-  if(verbose) message("Done.")

-  if (returnParams){

-    return(list(calls=A, confs=B, SNPQC=dev, batchQC=mean(diag(SS)), params=params))

-  }else{

-    return(list(calls=A, confs=B, SNPQC=dev, batchQC=mean(diag(SS))))

-  }

-}

-

-

-gtypeCallerRnm <- function(A, B, fIndex, mIndex, theCenters, theScales,

-                         theNs, Indexes, cIndexes, nIndexes,

-                         ncIndexes, SMEDIAN, knots, params, dft,

-                         probs, trim){

-

-  stopifnot(!missing(A), !missing(B), dim(A)==dim(B),

-            nrow(A)==nrow(theCenters), nrow(A)==nrow(theScales),

-            nrow(A) == nrow(theNs), length(knots)==3, nrow(params)==4,

-            ncol(params)==ncol(A), length(trim)==1, length(probs)==3)

-

-  ## make code robust

-  ## check types before passing to C

-  

-  .Call("gtypeCallerPart1nm", A, B,

-        as.integer(fIndex), as.integer(mIndex),

-        as.numeric(theCenters), as.numeric(theScales),

-        as.integer(theNs), lapply(Indexes, as.integer), lapply(cIndexes, as.integer), as.integer(nIndexes), as.integer(ncIndexes),

-        as.numeric(SMEDIAN), as.numeric(knots), as.numeric(params),

-        as.integer(dft), as.numeric(probs), as.numeric(trim),

-        PACKAGE="crlmm")

-  

-}

-

-gtypeCallerR2nm <- function(A, B, fIndex, mIndex, theCenters, theScales,

-                         theNs, Indexes, cIndexes, nIndexes,

-                         ncIndexes, SMEDIAN, knots, params, dft,

-                         probs, trim, noTraining, noInfo){

-

-  stopifnot(!missing(A), !missing(B), dim(A)==dim(B),

-            nrow(A)==nrow(theCenters), nrow(A)==nrow(theScales),

-            nrow(A) == nrow(theNs), length(knots)==3, nrow(params)==4,

-            ncol(params)==ncol(A), length(trim)==1, length(probs)==3)

-

-  .Call("gtypeCallerPart2nm", A, B,

-        as.integer(fIndex), as.integer(mIndex),

-        as.numeric(theCenters), as.numeric(theScales),

-        as.integer(theNs), Indexes, cIndexes, nIndexes, ncIndexes,

-        as.numeric(SMEDIAN), as.numeric(knots), as.numeric(params),

-        as.integer(dft), as.numeric(probs), as.numeric(trim),

-        as.integer(noTraining), as.integer(noInfo), PACKAGE="crlmm")

-  

-}

-

@@ -1,85 +1,85 @@

 snprma <- function(filenames, mixtureSampleSize=10^5, fitMixture=FALSE, eps=0.1, verbose=TRUE, seed=1, cdfName, sns){

-	if (missing(sns)) sns <- basename(filenames)

-	##ADD CHECK TO SEE IF LOADED

-	if (missing(cdfName))

-		cdfName <- read.celfile.header(filenames[1])$cdfName

-	##  stuffDir <- changeToCrlmmAnnotationName(cdfName)

-	pkgname <- getCrlmmAnnotationName(cdfName)

-	if(!require(pkgname, character.only=TRUE, quietly=!verbose)){

-		suggCall <- paste("library(", pkgname, ", lib.loc='/Altern/Lib/Loc')", sep="")

-		msg <- paste("If", pkgname, "is installed on an alternative location, please load it manually by using", suggCall)

-		message(strwrap(msg))

-		stop("Package ", pkgname, " could not be found.")

-		rm(suggCall, msg)

-	}

-	if(verbose) message("Loading annotations and mixture model parameters.")

-	loader("preprocStuff.rda", .crlmmPkgEnv, pkgname)

-	loader("genotypeStuff.rda", .crlmmPkgEnv, pkgname)

-	loader("mixtureStuff.rda", .crlmmPkgEnv, pkgname)

-	autosomeIndex <- getVarInEnv("autosomeIndex")

-	pnsa <- getVarInEnv("pnsa")

-	pnsb <- getVarInEnv("pnsb")

-	fid <- getVarInEnv("fid")

-	reference <- getVarInEnv("reference")

-	aIndex <- getVarInEnv("aIndex")

-	bIndex <- getVarInEnv("bIndex")

-	SMEDIAN <- getVarInEnv("SMEDIAN")

-	theKnots <- getVarInEnv("theKnots")

-	gns <- getVarInEnv("gns")

-

-	##We will read each cel file, summarize, and run EM one by one

-	##We will save parameters of EM to use later

-	mixtureParams <- matrix(0, 4, length(filenames))

-	SNR <- vector("numeric", length(filenames))

-	SKW <- vector("numeric", length(filenames))

-

-	## This is the sample for the fitting of splines

-	## BC: I like better the idea of the user passing the seed,

-	##     because this might intefere with other analyses

-	##     (like what happened to GCRMA)

-	set.seed(seed)

-	idx <- sort(sample(autosomeIndex, mixtureSampleSize))

-	##S will hold (A+B)/2 and M will hold A-B

-	##NOTE: We actually dont need to save S. Only for pics etc...

-	##f is the correction. we save to avoid recomputing

-	A <- matrix(as.integer(0), length(pnsa), length(filenames))

-	B <- matrix(as.integer(0), length(pnsb), length(filenames))

+  if (missing(sns)) sns <- basename(filenames)

+  ##ADD CHECK TO SEE IF LOADED

+  if (missing(cdfName))

+    cdfName <- read.celfile.header(filenames[1])$cdfName

+  ##  stuffDir <- changeToCrlmmAnnotationName(cdfName)

+  pkgname <- getCrlmmAnnotationName(cdfName)

+  if(!require(pkgname, character.only=TRUE, quietly=!verbose)){

+    suggCall <- paste("library(", pkgname, ", lib.loc='/Altern/Lib/Loc')", sep="")

+    msg <- paste("If", pkgname, "is installed on an alternative location, please load it manually by using", suggCall)

+    message(strwrap(msg))

+    stop("Package ", pkgname, " could not be found.")

+    rm(suggCall, msg)

+  }

+  if(verbose) message("Loading annotations and mixture model parameters.")

+  loader("preprocStuff.rda", .crlmmPkgEnv, pkgname)

+  loader("genotypeStuff.rda", .crlmmPkgEnv, pkgname)

+  loader("mixtureStuff.rda", .crlmmPkgEnv, pkgname)

+  autosomeIndex <- getVarInEnv("autosomeIndex")

+  pnsa <- getVarInEnv("pnsa")

+  pnsb <- getVarInEnv("pnsb")

+  fid <- getVarInEnv("fid")

+  reference <- getVarInEnv("reference")

+  aIndex <- getVarInEnv("aIndex")

+  bIndex <- getVarInEnv("bIndex")

+  SMEDIAN <- getVarInEnv("SMEDIAN")

+  theKnots <- getVarInEnv("theKnots")

+  gns <- getVarInEnv("gns")

+  

+  ##We will read each cel file, summarize, and run EM one by one

+  ##We will save parameters of EM to use later

+  mixtureParams <- matrix(0, 4, length(filenames))

+  SNR <- vector("numeric", length(filenames))

+  SKW <- vector("numeric", length(filenames))

-	if(verbose){

-		message("Processing ", length(filenames), " files.")

-		if (getRversion() > '2.7.0') pb <- txtProgressBar(min=0, max=length(filenames), style=3)

-	}

-	##We start looping throug cel files

-	idx2 <- sample(length(fid), 10^5) ##for skewness. no need to do everything

-	for(i in seq(along=filenames)){

-		y <- as.matrix(read.celfile(filenames[i], intensity.means.only=TRUE)[["INTENSITY"]][["MEAN"]][fid])

-		x <- log2(y[idx2])

-		SKW[i] <- mean((x-mean(x))^3)/(sd(x)^3)

-		rm(x)

-		y <- normalize.quantiles.use.target(y, target=reference)

-		A[, i] <- intMedianSummaries(y[aIndex, 1, drop=FALSE], pnsa)

-		B[, i] <- intMedianSummaries(y[bIndex, 1, drop=FALSE], pnsb)

-		##Now to fit the EM

-		if(fitMixture){

-			S <- (log2(A[idx, i])+log2(B[idx, i]))/2 - SMEDIAN

-			M <- log2(A[idx, i])-log2(B[idx, i])

+  ## This is the sample for the fitting of splines

+  ## BC: I like better the idea of the user passing the seed,

+  ##     because this might intefere with other analyses

+  ##     (like what happened to GCRMA)

+  set.seed(seed)

+  idx <- sort(sample(autosomeIndex, mixtureSampleSize))

+  ##S will hold (A+B)/2 and M will hold A-B

+  ##NOTE: We actually dont need to save S. Only for pics etc...

+  ##f is the correction. we save to avoid recomputing

+  A <- matrix(as.integer(0), length(pnsa), length(filenames))

+  B <- matrix(as.integer(0), length(pnsb), length(filenames))

+  

+  if(verbose){

+    message("Processing ", length(filenames), " files.")

+    if (getRversion() > '2.7.0') pb <- txtProgressBar(min=0, max=length(filenames), style=3)

+  }

+  ##We start looping throug cel files

+  idx2 <- sample(length(fid), 10^5) ##for skewness. no need to do everything

+  for(i in seq(along=filenames)){

+    y <- as.matrix(read.celfile(filenames[i], intensity.means.only=TRUE)[["INTENSITY"]][["MEAN"]][fid])

+    x <- log2(y[idx2])

+    SKW[i] <- mean((x-mean(x))^3)/(sd(x)^3)

+    rm(x)

+    y <- normalize.quantiles.use.target(y, target=reference)

+    A[, i] <- intMedianSummaries(y[aIndex, 1, drop=FALSE], pnsa)

+    B[, i] <- intMedianSummaries(y[bIndex, 1, drop=FALSE], pnsb)

+    ##Now to fit the EM

+    if(fitMixture){

+      S <- (log2(A[idx, i])+log2(B[idx, i]))/2 - SMEDIAN

+      M <- log2(A[idx, i])-log2(B[idx, i])

-			##we need to test the choice of eps.. it is not the max diff between funcs

-			tmp <- fitAffySnpMixture56(S, M, theKnots, eps=eps)

-			

-			mixtureParams[, i] <- tmp[["coef"]]

-			SNR[i] <- tmp[["medF1"]]^2/(tmp[["sigma1"]]^2+tmp[["sigma2"]]^2)

-		}

-		if (verbose)

-			if (getRversion() > '2.7.0') setTxtProgressBar(pb, i)

-			else cat(".")

-	}

-	if (verbose)

-		if (getRversion() > '2.7.0') close(pb)

-		else cat("\n")

-	if (!fitMixture) SNR <- mixtureParams <- NA

-	## gns comes from preprocStuff.rda

-	list(A=A, B=B, sns=sns, gns=gns, SNR=SNR, SKW=SKW, mixtureParams=mixtureParams, cdfName=cdfName)

+      ##we need to test the choice of eps.. it is not the max diff between funcs

+      tmp <- fitAffySnpMixture56(S, M, theKnots, eps=eps)

+      

+      mixtureParams[, i] <- tmp[["coef"]]

+      SNR[i] <- tmp[["medF1"]]^2/(tmp[["sigma1"]]^2+tmp[["sigma2"]]^2)

+    }

+    if (verbose)

+      if (getRversion() > '2.7.0') setTxtProgressBar(pb, i)

+      else cat(".")

+  }

+  if (verbose)

+    if (getRversion() > '2.7.0') close(pb)

+    else cat("\n")

+  if (!fitMixture) SNR <- mixtureParams <- NA

+  ## gns comes from preprocStuff.rda