R/methods-CNSet.R
c75c54c0
 setMethod("posteriorMean", signature(object="CNSet"), function(object) assayDataElement(object, "posteriorMean"))
 setReplaceMethod("posteriorMean", signature(object="CNSet", value="matrix"), function(object, value) assayDataElementReplace(object, "posteriorMean", value))
3ad9a188
 ##setReplaceMethod("mixtureParams", signature(object="CNSet", value="ANY"), function(object, value) object@mixtureParams <- mixtureParams)
c75c54c0
 
97e84a80
 linearParamElementReplace <- function(obj, elt, value) {
6846d583
     storage.mode <- storageMode(batchStatistics(obj))
97e84a80
     switch(storage.mode,
            "lockedEnvironment" = {
6846d583
                aData <- copyEnv(batchStatistics(obj))
97e84a80
                if (is.null(value)) rm(list=elt, envir=aData)
                else aData[[elt]] <- value
                Biobase:::assayDataEnvLock(aData)
6846d583
                batchStatistics(obj) <- aData
97e84a80
            },
            "environment" = {
6846d583
                if (is.null(value)) rm(list=elt, envir=batchStatistics(obj))
                else batchStatistics(obj)[[elt]] <- value
97e84a80
            },
6846d583
            list = batchStatistics(obj)[[elt]] <- value)
97e84a80
     obj
 }
 
6846d583
 ## parameters
0198a9ad
 ## allele A
 ##   autosome SNPs
 ##   autosome NPs
 ##   chromosome X NPs for women
 C1 <- function(object, marker.index, batch.index, sample.index){
a268cdc3
 	acn <- matrix(NA, nrow=length(marker.index), ncol=length(sample.index))
0198a9ad
 	for(k in seq_along(batch.index)){
 		l <- batch.index[k]
453e688a
 		## calculate cn for all the samples that are in this batch
0198a9ad
 		jj <- sample.index[as.character(batch(object))[sample.index] == batchNames(object)[l]]
 		bg <- nuA(object)[marker.index, l]
 		slope <- phiA(object)[marker.index, l]
7e59dced
 		I <- as.matrix(A(object)[marker.index, jj])
a268cdc3
 		acn[, match(jj, sample.index)] <- 1/slope * (I - bg)
0198a9ad
 	}
9568a401
 	return(as.matrix(acn))
0198a9ad
 }
 
197c7911
 ## allele B  (treated allele 'A' for chromosome X NPs)
0198a9ad
 ##   autosome SNPs
 ##   chromosome X for male nonpolymorphic markers
 C2 <- function(object, marker.index, batch.index, sample.index, NP.X=FALSE){
a268cdc3
 	acn <- matrix(NA, nrow=length(marker.index), ncol=length(sample.index))
0198a9ad
 	for(k in seq_along(batch.index)){
 		l <- batch.index[k]
 		jj <- sample.index[as.character(batch(object))[sample.index] == batchNames(object)[l]]
 		bg <- nuB(object)[marker.index, l]
 		slope <- phiB(object)[marker.index, l]
 		if(!NP.X){
7e59dced
 			I <- as.matrix(B(object)[marker.index, jj])
 		} else I <- as.matrix(A(object)[marker.index, jj])
a268cdc3
 		acn[, match(jj, sample.index)] <- 1/slope * (I - bg)
0198a9ad
 	}
a268cdc3
 ##	if(length(acn) > 1){
 ##		acn <- do.call("cbind", acn)
 ##	} else acn <- acn[[1]]
9568a401
 	return(as.matrix(acn))
0198a9ad
 }
 
 ## Chromosome X SNPs
 C3 <- function(object, allele, marker.index, batch.index, sample.index){
a268cdc3
 ##	acn <- vector("list", length(batch.index))
 	acn <- matrix(NA, nrow=length(marker.index), ncol=length(sample.index))
0198a9ad
 	for(k in seq_along(batch.index)){
 		l <- batch.index[k]
 		##j <- which(as.character(batch(object))[sample.index] == batchNames(object)[l])
 		jj <- sample.index[as.character(batch(object))[sample.index] == batchNames(object)[l]]
1813938e
 		##phiA2 and phiB2 are not always estimable  -- need both men and women
0198a9ad
 		phiA2 <- phiPrimeA(object)[marker.index, l]
 		phiB2 <- phiPrimeB(object)[marker.index, l]
 		phiA <- phiA(object)[marker.index, l]
 		phiB <- phiB(object)[marker.index, l]
 		nuA <- nuA(object)[marker.index, l]
89c5fe95
 		nuB <- nuB(object)[marker.index, l]
453e688a
 		phiA <- phiA(object)[marker.index, l]
7e59dced
 		IA <- as.matrix(A(object)[marker.index, jj])
 		IB <- as.matrix(B(object)[marker.index, jj])
453e688a
 		## I = nu + acn * phi
 		## acn = 1/phi* (I-nu)
89c5fe95
 		phistar <- phiB2/phiA
 		tmp <- (IB - nuB - phistar*IA + phistar*nuA)/phiB
 		CB <- tmp/(1-phistar*phiA2/phiB)
1813938e
 		index <- which(is.na(phiB2) | is.na(phiA2))
 		if(length(index) > 0){
 			cb <- 1/phiB[index] * (IB[index, ] - nuB[index])
 			CB[index, ] <- cb
 		}
0198a9ad
 		if(allele == "B"){
a268cdc3
 			acn[, match(jj, sample.index)] <- CB
 			##acn[[k]] <- CB
0198a9ad
 		}
 		if(allele == "A"){
1813938e
 			CA <- (IA-nuA-phiA2*CB)/phiA
 			if(length(index) > 0){
 				ca <- 1/phiA[index] * (IA[index, ] - nuA[index])
  				CA[index, ] <- ca
 			}
 			acn[, match(jj, sample.index)] <- CA
0198a9ad
 		}
 	}
a268cdc3
 ##	if(length(acn) > 1){
 ##		acn <- do.call("cbind", acn)
 ##	} else acn <- acn[[1]]
9568a401
 	return(as.matrix(acn))
0198a9ad
 }
 
8e2d9355
 
0198a9ad
 
518a2908
 
805ca33b
 ACN <- function(object, allele, i , j){
6846d583
 	if(missing(i) & missing(j)) stop("must specify rows (i) or columns (j)")
9568a401
 	is.ff <- is(calls(object), "ff") | is(calls(object), "ffdf")
0198a9ad
 	missing.i <- missing(i)
 	missing.j <- missing(j)
 	if(!missing.i){
 		is.ann <- !is.na(chromosome(object)[i])
 		is.X <- chromosome(object)[i]==23 & is.ann
 		is.auto <- chromosome(object)[i] < 23 & is.ann
 		is.snp <- isSnp(object)[i] & is.ann
 	} else{
 		is.ann <- !is.na(chromosome(object))
 		is.X <- chromosome(object)==23 & is.ann
 		is.auto <- chromosome(object) < 23 & is.ann
 		is.snp <- isSnp(object) & is.ann
 		i <- 1:nrow(object)
 	}
 	## Define batch.index and sample.index
 	if(!missing.j) {
063b3d14
 		batches <- unique(as.character(batch(object)[j]))
a268cdc3
 		##batches <- as.character(batch(object)[j])
0198a9ad
 		batch.index <- match(batches, batchNames(object))
 	} else {
 		batch.index <- seq_along(batchNames(object))
 		j <- 1:ncol(object)
 	}
 	nr <- length(i)
 	nc <- length(j)
 	acn <- matrix(NA, nr, nc)
a268cdc3
 	dimnames(acn) <- list(featureNames(object)[i],
 			      sampleNames(object)[j])
0198a9ad
 	if(allele == "A"){
 		if(is.ff){
 			open(nuA(object))
 			open(phiA(object))
 			open(A(object))
 		}
 		## --
 		## 4 types of markers for allele A
 		##--
 		## 1. autosomal SNPs or autosomal NPs
 		if(any(is.auto)){
 			auto.index <- which(is.auto)
 			marker.index <- i[is.auto]
 			acn[auto.index, ] <- C1(object, marker.index, batch.index, j)
 		}
 		if(any(is.X)){
 			##2. CHR X SNPs (men and women)
 			if(any(is.snp)){
 				if(is.ff) {
 					open(phiPrimeA(object))
 					open(phiPrimeB(object))
 					open(phiB(object))
 					open(nuB(object))
 					open(B(object))
 				}
 				marker.index <- i[is.X & is.snp]
 				acn.index <- which(is.X & is.snp)
 				acn[acn.index, ] <- C3(object, allele="A", marker.index, batch.index, j)
 				if(is.ff) {
 					close(phiPrimeA(object))
 					close(phiPrimeB(object))
 					close(phiB(object))
 					close(nuB(object))
 					close(B(object))
 				}
 			}
a268cdc3
 			if(any(!is.snp)){
0198a9ad
 				marker.index <- i[is.X & !is.snp]
 				acn.index <- which(is.X & !is.snp)
 				acn[acn.index, ] <- NA
4191b5f9
 				female.index <- j[object$gender[j] == 2]
 				## 3. CHR X NPs: women
 				if(length(female.index) > 0){
 					female.batch.index <- match(unique(as.character(batch(object))[female.index]), batchNames(object))
 					jj <- which(object$gender[j] == 2)
 					acn[acn.index, jj] <- C1(object, marker.index, female.batch.index, female.index)
 				}
 				male.index <- j[object$gender[j] == 1]
 				if(length(male.index) > 0){
 					if(is.ff){
 						open(nuB(object))
 						open(phiB(object))
 					}
 					male.batch.index <- match(unique(as.character(batch(object))[male.index]), batchNames(object))
 					jj <- which(object$gender[j] == 1)
 					acn[acn.index, jj] <- C2(object, marker.index, male.batch.index, male.index, NP.X=TRUE)
 					if(is.ff){
 						close(nuB(object))
 						close(phiB(object))
 					}
 				}
6846d583
 			}
 		}
0198a9ad
 		if(is.ff){
 			close(nuA(object))
 			close(phiA(object))
 			close(A(object))
805ca33b
 		}
 	}
0198a9ad
 	if(allele == "B"){
 		if(is.ff){
 			open(nuB(object))
 			open(phiB(object))
 			open(B(object))
805ca33b
 		}
a268cdc3
 		if(any(!is.snp)){
 			acn.index <- which(!is.snp)
 			acn[acn.index, ] <- 0
 		}
89c5fe95
 		if(any(is.auto)){
5133e591
 			auto.index <- which(is.auto & is.snp)
 			if(length(auto.index) > 0){
c7175779
 				marker.index <- i[auto.index]
5133e591
 				acn[auto.index, ] <- C2(object, marker.index, batch.index, j)
 			}
d5b39e21
 		}
89c5fe95
 		if(any(is.X)){
5cfbfadf
 			if(is.ff){
 				open(phiPrimeA(object))
 				open(phiPrimeB(object))
 				open(phiA(object))
 				open(nuA(object))
 				open(A(object))
 			}
 			marker.index <- i[is.X & is.snp]
 			acn.index <- which(is.X & is.snp)
 			acn[acn.index, ] <- C3(object, allele="B", marker.index, batch.index, j)
 			if(is.ff){
 				close(phiPrimeA(object))
 				close(phiPrimeB(object))
 				close(phiA(object))
 				close(nuA(object))
 				close(A(object))
0198a9ad
 			}
89c5fe95
 			if(any(!is.snp)){
 				acn.index <- which(!is.snp)
 				marker.index <- i[!is.snp]
 				acn[acn.index, ] <- 0
 			}
0198a9ad
 		}
f3738959
 	}
 	if(is.ff){
 		close(nuB(object))
 		close(phiB(object))
 		close(B(object))
805ca33b
 	}
 	return(acn)
 }
ddc94a58
 
69ec8644
 setMethod("CA", signature=signature(object="CNSet"),
 	  function(object, ...){
 		  ca <- ACN(object, allele="A", ...)
0198a9ad
 		  ca[ca < 0] <- 0
66900fea
 		  ca[ca > 5] <- 5
ddc94a58
 		  return(ca)
2c958cfa
 	  })
69ec8644
 setMethod("CB", signature=signature(object="CNSet"),
 	  function(object, ...) {
 		  cb <- ACN(object, allele="B", ...)
0198a9ad
 		  cb[cb < 0] <- 0
66900fea
 		  cb[cb > 5] <- 5
805ca33b
 		  return(cb)
 	  })
752632c4
 
6b75ec77
 setMethod("totalCopynumber", signature=signature(object="CNSet"),
 	  function(object, ...){
fda11070
 		  ca <- CA(object, ...)
 		  cb <- CB(object, ...)
 		  return(ca+cb)
 	  })
70878f91
 rawCopynumber <- totalCopynumber
69ec8644
 
e516fadb
 setMethod("posteriorProbability", signature(object="CNSet"),
96fb65cc
 	  function(object, predictRegion, copyNumber=0:4, w){
 		  if(missing(w)) w <- rep(1/length(copyNumber),length(copyNumber))
 		  stopifnot(sum(w)==1)
e516fadb
 		  logI <- array(NA, dim=c(nrow(object), ncol(object), 2), dimnames=list(NULL, NULL, LETTERS[1:2]))
 		  getIntensity <- function(object){
 			  logI[, , 1] <- log2(A(object))
 			  logI[, , 2] <- log2(B(object))
 			  return(logI)
 		  }
 		  logI <- getIntensity(object)
 		  ##gts <- lapply(as.list(0:4), genotypes)
 		  prob <- array(NA, dim=c(nrow(object), ncol(object), length(copyNumber)),
 				dimnames=list(NULL,NULL, paste("copynumber",copyNumber, sep="")))
96fb65cc
 		  bns <- batchNames(object)
b398d6f8
 		  snp.index <- which(isSnp(object))
 		  if(length(snp.index) > 0){
 			  for(i in seq_along(copyNumber)){
 				  G <- genotypes(copyNumber[i])
 				  P <- array(NA, dim=c(length(snp.index), ncol(object), length(G)),
 					     dimnames=list(NULL,NULL,G))
 				  for(g in seq_along(G)){
 					  gt <- G[g]
 					  for(j in seq_along(bns)){
 						  this.batch <- bns[j]
 						  sample.index <- which(batch(object) == this.batch)
 						  mu <- predictRegion[[gt]]$mu[snp.index, , this.batch, drop=FALSE]
 						  dim(mu) <- dim(mu)[1:2]
 						  if(all(is.na(mu))){
 							  P[, sample.index, gt] <- NA
 						  } else {
 							  Sigma <- predictRegion[[gt]]$cov[snp.index, , this.batch, drop=FALSE]
 							  dim(Sigma) <- dim(Sigma)[1:2]
 							  P[, sample.index, gt] <- dbvn(x=logI[snp.index, sample.index, , drop=FALSE], mu=mu, Sigma=Sigma)
 						  }
96fb65cc
 					  }
 				  }
b398d6f8
 				  ## the marginal probability for the total copy number
 				  ##  -- integrate out the probability for the different genotypes
 				  for(j in 1:ncol(P)){
 					  PP <- P[, j, , drop=FALSE]
 					  dim(PP) <- dim(PP)[c(1, 3)]
 					  prob[snp.index, j, i] <- rowSums(PP, na.rm=TRUE)
 				  }
e516fadb
 			  }
b398d6f8
 		  } ## length(snp.index) > 0
 		  ##---------------------------------------------------------------------------
 		  ##
 		  ##  nonpolymorphic markers
 		  ##
 		  ##---------------------------------------------------------------------------
 		  np.index <- which(!isSnp(object))
 		  if(length(np.index) > 0){
 			  for(i in seq_along(copyNumber)){
 				  G <- genotypes(copyNumber[i], is.snp=FALSE)
 				  P <- array(NA, dim=c(length(np.index), ncol(object), length(G)),
 					     dimnames=list(NULL,NULL,G))
 				  for(g in seq_along(G)){
 					  gt <- G[g]
 					  for(j in seq_along(bns)){
 						  this.batch <- bns[j]
 						  sample.index <- which(batch(object) == this.batch)
 						  mu <- predictRegion[[gt]]$mu[np.index, 1, this.batch, drop=FALSE]
 						  dim(mu) <- dim(mu)[1]
 						  if(all(is.na(mu))){
 							  P[, sample.index, gt] <- NA
 						  } else {
 							  Sigma <- predictRegion[[gt]]$cov[np.index, 1, this.batch, drop=FALSE]
 							  dim(Sigma) <- dim(Sigma)[1]
 							  P[, sample.index, gt] <- dnorm(logI[np.index, sample.index, 1], mean=mu, sd=Sigma)
 						  }
 					  } ## j in seq_along(bns)
 				  } ## g in seq_along(G)
 				  ## the marginal probability for the total copy number
 				  ##  -- integrate out the probability for the different genotypes
 				  prob[np.index, , i] <- P
 			  } ## seq_along(copyNumber)
 		  } ## length(np.index) > 0
96fb65cc
 		  ##constant.  Probs across copy number states must
 		  ##sum to one.  nc <- apply(prob, c(1,3), sum,
 		  ##na.rm=TRUE)
 		  wm <- matrix(w, nrow(object), length(copyNumber), byrow=TRUE)
e516fadb
 		  nc <- matrix(NA, nrow(object), ncol(object))
 		  for(j in seq_len(ncol(object))){
96fb65cc
 			  pp <- prob[, j, ] * wm
 			  nc <- rowSums(pp, na.rm=TRUE)
 			  prob[, j, ] <- pp/nc
e516fadb
 		  }
 		  return(prob)
 	  })
3e542411
 
e516fadb
 setMethod("calculatePosteriorMean", signature(object="CNSet"),
96fb65cc
 	  function(object, posteriorProb, copyNumber=0:4, ...){
e516fadb
 		  stopifnot(dim(posteriorProb)[[3]] == length(copyNumber))
 		  pm <- matrix(0, nrow(object), ncol(object))
 		  for(i in seq_along(copyNumber)){
96fb65cc
 			  pm <- pm + posteriorProb[, , i] * copyNumber[i]
e516fadb
 		  }
 		  return(pm)
 	  })
3e542411
 
ad3cc706
 ##.bivariateCenter <- function(nu, phi){
 ##			  ##  lexical scope for mus, CA, CB
 ##			  if(CA <= 2 & CB <= 2 & (CA+CB) < 4){
 ##				  mus[,1, ] <- log2(nu[, 1, ] + CA *
 ##						    phi[, 1, ])
 ##				  mus[,2, ] <- log2(nu[, 2, ] + CB *
 ##						    phi[, 2, ])
 ##			  } else { ## CA > 2
 ##				  if(CA > 2){
 ##					  theta <- pi/4*Sigma[,2,]
 ##					  shiftA <- CA/4*phi[, 1, ] * cos(theta)
 ##					  shiftB <- CA/4*phi[, 1, ] * sin(theta)
 ##					  mus[, 1, ] <- log2(nu[, 1, ] + 2 * phi[,1,]+shiftA)
 ##					  mus[, 2, ] <- log2(nu[, 2, ] + CB *phi[,2,]+shiftB)
 ##				  }
 ##				  if(CB > 2){
 ##					  ## CB > 2
 ##					  theta <- pi/2-pi/4*Sigma[,2,]
 ##					  shiftA <- CB/4*phi[, 2, ] * cos(theta)
 ##					  shiftB <- CB/4*phi[, 2, ] * sin(theta)
 ##					  mus[, 1, ] <- log2(nu[, 1, ] + CA*phi[,1,]+shiftA)
 ##					  mus[, 2, ] <- log2(nu[, 2, ]+ 2*phi[,2,]+shiftB)
 ##				  }
 ##				  if(CA == 2 & CB == 2){
 ##					  mus[, 1, ] <- log2(nu[, 1, ] + 1/2*CA*phi[,1,])
 ##					  mus[, 2, ] <- log2(nu[, 2, ]+ 1/2*CB*phi[,2,])
 ##				  }
 ##			  }
 ##			  mus
 ##		  }
3e542411
 
e516fadb
 ## for a given copy number, return a named list of bivariate normal prediction regions
 ##   - elements of list are named by genotype
 ##   'AAA'  list should have 'mu' and 'Sigma'
 ##                mu is a R x 2 matrix, R is number of features
 ##                Sigma is a R x 4 matrix.  Each row can be coerced to a 2 x 2 matrix
 ##          For nonpolymorphic markers, 2nd column of mu is NA and elements 2-4 of Sigma are NA
 ##
 setMethod("predictionRegion", signature(object="CNSet", copyNumber="integer"),
 	  function(object, copyNumber=0:4){
 		  ## would it be more efficient to store as an array?
 		  ## mu: features x genotype x allele
 		  ## Sigma: features x genotype x covariance
 		  stopifnot(all(copyNumber %in% 0:4))
0a5dbc06
 		  getNu <- function(object){
 			  nu[, 1, ] <- nuA(object)
 			  nu[, 2, ] <- nuB(object)
 			  nu
 		  }
 		  getPhi <- function(object){
 			  phi[,1,] <- phiA(object)
 			  phi[,2,] <- phiB(object)
 			  phi
 		  }
e516fadb
 		  gts <- lapply(as.list(copyNumber), genotypes)
 		  nms <- unlist(gts)
 		  res <- vector("list", length(nms))
 		  ##names(res) <- paste("copyNumber", copyNumber, sep="")
 		  names(res) <- nms
0a5dbc06
 		  bnms <- batchNames(object)
 		  nu <- array(NA, dim=c(nrow(object), 2, length(bnms)))
 		  phi <- array(NA, dim=c(nrow(object), 2, length(bnms)))
e516fadb
 		  nu <- getNu(object)
 		  phi <- getPhi(object)
0a5dbc06
 		  ##mus <- matrix(NA, nrow(nu), 2, dimnames=list(NULL, LETTERS[1:2]))
 		  mus <- array(NA, dim=c(nrow(nu), 2, length(bnms)),
 			       dimnames=list(NULL, LETTERS[1:2],
 			       bnms))
 		  ## Sigma <- matrix(NA, nrow(nu), 3)
 		  Sigma <- array(NA, dim=c(nrow(nu), 3, length(bnms)),
 				 dimnames=list(NULL, c("varA", "cor", "varB"),
 				 bnms))
e516fadb
 		  bivariateCenter <- function(nu, phi){
 			  ##  lexical scope for mus, CA, CB
97a2ac47
 			  mus[,1, ] <- log2(nu[, 1, ] + CA *
 					    phi[, 1, ])
 			  mus[,2, ] <- log2(nu[, 2, ] + CB *
 					    phi[, 2, ])
 			  return(mus)
e516fadb
 		  }
b398d6f8
 		  np.index <- which(!isSnp(object))
e516fadb
 		  for(i in seq_along(copyNumber)){
 			  G <- genotypes(copyNumber[i])
 			  tmp <- vector("list", length(G))
 			  names(tmp) <- G
 			  CN <- copyNumber[i]
 			  for(g in seq_along(G)){
 				  gt <- G[g]
 				  CB <- g-1
 				  CA <- CN-CB
 				  gt.corr <- genotypeCorrelation(gt)
 				  nma <- ifelse(CA == 0, "tau2A.BB", "tau2A.AA")
 				  nmb <- ifelse(CB == 0, "tau2B.AA", "tau2B.BB")
0a5dbc06
 				  Sigma[, 1, ] <- getVar(object, nma)
 				  Sigma[, 3, ] <- getVar(object, nmb)
 				  Sigma[, 2, ] <- getCor(object, gt.corr)
b398d6f8
 				  if(length(np.index) > 0){
 					  Sigma[, 1, ] <- getVar(object, "tau2A.AA")
 					  Sigma[np.index, 2, ] <- NA
 				  }
97a2ac47
 				  res[[gt]]$mu <- bivariateCenter(nu,phi)
 				  ## adjust correlation
 ##				  if(CA == 0 | CB == 0){
 ##					  Sigma[,2,] <- 0
 ##				  }
e516fadb
 				  res[[gt]]$cov <- Sigma
 			  }
 		  }
97a2ac47
 		  res <- as(res, "PredictionRegion")
e516fadb
 		  return(res)
 	  })
0a5dbc06
 
 
97a2ac47
 
0a5dbc06
 setMethod("xyplotcrlmm", signature(x="formula", data="CNSet", predictRegion="list"),
 	  function(x, data, predictRegion, ...){
96fb65cc
 		  fns <- featureNames(data)
 		  fns <- matrix(fns, nrow(data), ncol(data), byrow=FALSE)
 		  fns <- as.character(fns)
 		  df <- list(A=as.numeric(log2(A(data))),
 			     B=as.numeric(log2(B(data))),
 			     gt=as.integer(calls(data)),
 			     gt.conf=as.numeric(confs(data)),
 			     snpid=fns)#, snp=snpId)
 		  df <- as.data.frame(df)
99377638
 		  df$snpid <- factor(df$snpid, levels=unique(df$snpid), ordered=TRUE)
96fb65cc
 		  bns <- batchNames(data)
 		  predictRegion <- lapply(predictRegion, function(x, bns){
 			  batch.index <- match(bns, dimnames(x$mu)[[3]])
 			  x$mu <- x$mu[, , batch.index, drop=FALSE]
 			  x$cov <- x$cov[, , batch.index, drop=FALSE]
 			  return(x)
 		  }, bns=bns)
97a2ac47
 		  ##predictRegion is an argument of ABpanel
0a5dbc06
 		  xyplot(x, df, predictRegion=predictRegion, ...)
 	  })
32b8047e
 
0a5dbc06
 setMethod("xyplot", signature(x="formula", data="CNSet"),
 	  function(x, data, ...){
8be0a109
 		  if("predictRegion" %in% names(list(...))){
 			  xyplotcrlmm(x, data, ...)
 		  } else{
 			  callNextMethod()
 		  }
0a5dbc06
 })
4c4e9303
 
 setMethod("calculateRBaf", signature(object="CNSet"),
063b3d14
 	  function(object, batch.name, chrom){
 		  calculateRBafCNSet(object, batch.name, chrom)
 	  })
4c4e9303
 
063b3d14
 calculateRBafCNSet <- function(object, batch.name, chrom){
2bc29627
 	if(missing(batch.name)) {
 		batch.name <- batchNames(object)
 		if("grandMean" %in% batch.name)
 			batch.name <- batch.name[-length(batch.name)]
 	}
063b3d14
 	if(missing(chrom)) chrom <- unique(chromosome(object))
 	if(!(all(batch.name %in% batchNames(object)))) stop("batch.name must be belong to batchNames(object)")
 	chr <- chromosome(object)
 	valid.chrs <- chr <= 23 & chr %in% chrom
 	index <- which(valid.chrs)
 	indexlist <- split(index, chr[index])
 	J <- which(batch(object) %in% batch.name)
 	sns <- sampleNames(object)[J]
a3b625d4
 	sampleindex <- split(J, factor(batch(object)[J], levels=unique(batch(object)[J])))
063b3d14
 	if(!all(valid.chrs)) warning("Only computing log R ratios and BAFs for autosomes and chr X")
 	## if ff package is loaded, these will be ff objects
 	chr <- names(indexlist)
 	rlist <- blist <- vector("list", length(indexlist))
2bc29627
 	path <- ldPath()
 	processByChromosome <- function(i, chr, path){
 		ldPath(path)
 		nr <- length(i)
 		##CHR <- names(i)
 		bafname <- paste("baf_chr", chr, sep="")
 		rname <- paste("lrr_chr", chr, sep="")
063b3d14
 		bmatrix <- initializeBigMatrix(bafname, nr=nr, nc=length(sns), vmode="integer")
 		rmatrix <- initializeBigMatrix(rname, nr=nr, nc=length(sns), vmode="integer")
 		colnames(rmatrix) <- colnames(bmatrix) <- sns
 		## put rownames in order of physical position
2bc29627
 		ix <- order(position(object)[i])
 		i <- i[ix]
 		rownames(rmatrix) <- rownames(bmatrix) <- featureNames(object)[i]
063b3d14
 		for(j in seq_along(sampleindex)){
 			bname <- batch.name[j]
 			J <- sampleindex[[j]]
2bc29627
 			res <- crlmm:::calculateRTheta(object=object,
 						       batch.name=bname,
 						       feature.index=i)
063b3d14
 			k <- match(sampleNames(object)[J], sns)
 			bmatrix[, k] <- res[["baf"]]
 			rmatrix[, k] <- res[["lrr"]]
 		}
2bc29627
 		list(bmatrix, rmatrix)
 	}
 	## calcualte R BAF by chromosome
 	if(isPackageLoaded("ff")){
 		pkgs <- c("oligoClasses", "ff", "Biobase", "crlmm")
 	} else pkgs <- c("oligoClasses", "Biobase", "crlmm")
5b0470a6
 	i <- NULL
2bc29627
 	res <- foreach(i=indexlist, chr=names(indexlist), .packages=pkgs) %dopar% {
 		processByChromosome(i, chr, path)
063b3d14
 	}
2bc29627
 	blist <- lapply(res, "[[", 1)
 	rlist <- lapply(res, "[[", 2)
 	res <- list(baf=blist, lrr=rlist)
063b3d14
 	return(res)
 }
4c4e9303
 
063b3d14
 ##setMethod("calculateRBaf", signature(object="CNSet"),
 ##	  function(object, batch.name){
 ##		  all.autosomes <- all(chromosome(object) < 23)
 ##		  if(!all.autosomes){
 ##			  stop("method currently only defined for chromosomes 1-22")
 ##		  }
 ##		  if(missing(batch.name)) batch.name <- batchNames(object)[1]
 ##		  stopifnot(batch.name %in% batchNames(object))
 ##		  if(length(batch.name) > 1){
 ##			  warning("only the first batch in batch.name processed")
 ##			  batch.name <- batch.name[1]
 ##		  }
 ##		  RTheta.aa <- calculateRTheta(object, "AA", batch.name)
 ##		  RTheta.ab <- calculateRTheta(object, "AB", batch.name)
 ##		  RTheta.bb <- calculateRTheta(object, "BB", batch.name)
 ##
 ##		  J <- which(batch(object) == batch.name)
 ##
 ##		  theta.aa <- matrix(RTheta.aa[, "theta"], nrow(object), length(J), byrow=FALSE)
 ##		  theta.ab <- matrix(RTheta.ab[, "theta"], nrow(object), length(J), byrow=FALSE)
 ##		  theta.bb <- matrix(RTheta.bb[, "theta"], nrow(object), length(J), byrow=FALSE)
 ##
 ##		  a <- A(object)[, J, drop=FALSE]
 ##		  b <- B(object)[, J, drop=FALSE] ## NA's for b where nonpolymorphic
 ##		  is.np <- !isSnp(object)
 ##		  ##b[is.np, ] <- a[is.np, ]
 ##		  b[is.np, ] <- 0L
 ##		  dns <- dimnames(a)
 ##		  dimnames(a) <- dimnames(b) <- NULL
 ##		  obs.theta <- atan2(b, a)*2/pi
 ##
 ##		  lessAA <- obs.theta < theta.aa
 ##		  lessAB <- obs.theta < theta.ab
 ##		  lessBB <- obs.theta < theta.bb
 ##		  grAA <- !lessAA
 ##		  grAB <- !lessAB
 ##		  grBB <- !lessBB
 ##		  not.na <- !is.na(theta.aa)
 ##		  I1 <- grAA & lessAB & not.na
 ##		  I2 <- grAB & lessBB & not.na
 ##
 ##		  bf <- matrix(NA, nrow(object), ncol(a))
 ##		  bf[I1] <- 0.5 * ((obs.theta-theta.aa)/(theta.ab-theta.aa))[I1]
 ##		  bf[I2] <- (.5 * (obs.theta - theta.ab) / (theta.bb - theta.ab))[I2] + 0.5
 ##		  bf[lessAA] <- 0
 ##		  bf[grBB] <- 1
 ##
 ##		  r.expected <- matrix(NA, nrow(object), ncol(a))
 ##		  r.aa <- matrix(RTheta.aa[, "R"], nrow(object), length(J), byrow=FALSE)
 ##		  r.ab <- matrix(RTheta.ab[, "R"], nrow(object), length(J), byrow=FALSE)
 ##		  r.bb <- matrix(RTheta.bb[, "R"], nrow(object), length(J), byrow=FALSE)
 ##		  rm(RTheta.aa, RTheta.ab, RTheta.bb); gc()
 ##		  obs.r <- a+b
 ##
 ##		  lessAA <- lessAA & not.na
 ##		  grBB <- grBB & not.na
 ##		  tmp <- ((obs.theta - theta.aa) * (r.ab-r.aa)/(theta.ab-theta.aa))[I1] + r.aa[I1]
 ##		  r.expected[I1] <- tmp
 ##		  tmp <- ((obs.theta - theta.ab) * (r.bb - r.ab)/(theta.bb-theta.ab))[I2] + r.ab[I2]
 ##		  r.expected[I2] <- tmp
 ##		  r.expected[lessAA] <- r.aa[lessAA]
 ##		  r.expected[grBB] <- r.bb[grBB]
 ##		  index.np <- which(is.np)
 ##		  if(length(index.np) > 0){
 ##			  a.np <- A(object)[index.np, J, drop=FALSE]
 ##			  meds <- rowMedians(a.np, na.rm=TRUE)
 ##			  r.expected[index.np, ] <- matrix(meds, length(index.np), ncol(a.np))
 ##		  }
 ##		  lrr <- log2(obs.r/r.expected)
 ##
 ##		  dimnames(bf) <- dimnames(lrr) <- dns
 ##		  res <- list(baf=bf,
 ##			      lrr=lrr)
 ##		  return(res)
 ##	  })
3ad9a188
 
 
 setAs("CNSet", "oligoSetList", function(from, to){
 	constructOligoSetListFrom(from)
 })
 
 
 
 setMethod(OligoSetList, "CNSet", function(object,...){
 	constructOligoSetListFrom(object, ...)
 })
 setMethod(BafLrrSetList, "CNSet", function(object,...){
 	constructBafLrrSetListFrom(object, ...)
 })
 
 
 
 
 
 constructOligoSetListFrom <- function(object, ...){
 	##row.index <- seq_len(nrow(object))
 	##col.index <- seq_len(ncol(object))
 	is.lds <- ifelse(is(calls(object), "ff_matrix") | is(calls(object), "ffdf"), TRUE, FALSE)
 	if(is.lds) stopifnot(isPackageLoaded("ff"))
 	b.r <- calculateRBaf(object, ...)
 	b <- b.r[["baf"]]
 	r <- b.r[["lrr"]]
 	j <- match(colnames(r[[1]]), sampleNames(object))
 	rns <- lapply(r, rownames)
 	fDList <- foreach(featureid=rns) %do%{
 		featureData(object)[match(featureid, featureNames(object)), ]
 	}
 	names(fDList) <- sapply(fDList, function(x) chromosome(x)[1])
 	gtPlist <- gtlist <- vector("list", length(r))
 	for(i in seq_along(r)){
 		gtlist[[i]] <- initializeBigMatrix("call", nr=nrow(r[[i]]), nc=length(j), vmode="integer")
 		gtPlist[[i]] <- initializeBigMatrix("callPr", nr=nrow(r[[i]]), nc=length(j), vmode="integer")
 		featureid <- rownames(r[[i]])
 		ix <- match(featureid, featureNames(object))
 		rownames(gtPlist[[i]]) <- rownames(gtlist[[i]]) <- featureid
 		colnames(gtPlist[[i]]) <- colnames(gtlist[[i]]) <- colnames(r[[i]])
 		for(k in seq_along(j)){
 			gtlist[[i]][, k] <- calls(object)[ix, j[k]]
 			gtPlist[[i]][, k] <- snpCallProbability(object)[ix, j[k]]
 		}
 	}
 	ad <- AssayDataList(baf=b, copyNumber=r, call=gtlist, callProbability=gtPlist)
 	object <- new("oligoSetList",
 		      assayDataList=ad,
 		      featureDataList=fDList,
 		      chromosome=names(fDList),
 		      phenoData=phenoData(object)[j, ],
 		      annotation=annotation(object),
 		      genome=genomeBuild(object))
 	return(object)
 }
 
 
 
 constructBafLrrSetListFrom <- function(object, ...){
 	is.lds <- ifelse(is(calls(object), "ff_matrix") | is(calls(object), "ffdf"), TRUE, FALSE)
 	if(is.lds) stopifnot(isPackageLoaded("ff"))
 	b.r <- calculateRBaf(object, ...)
 	b <- b.r[["baf"]]
 	r <- b.r[["lrr"]]
 	j <- match(colnames(r[[1]]), sampleNames(object))
 	rns <- lapply(r, rownames)
 	featureid <- NULL
 	fDList <- foreach(featureid=rns) %do%{
 		featureData(object)[match(featureid, featureNames(object)), ]
 	}
 	names(fDList) <- sapply(fDList, function(x) chromosome(x)[1])
 	ad <- AssayDataList(baf=b, lrr=r)
 	object <- new("BafLrrSetList",
 		      assayDataList=ad,
 		      featureDataList=fDList,
 		      chromosome=names(fDList),
 		      phenoData=phenoData(object)[j, ],
 		      annotation=annotation(object),
 		      genome=genomeBuild(object))
 	return(object)
 }
 
 computeBR <- constructBafLrrSetListFrom