@@ -681,7 +681,7 @@ setMethod(OligoSetList, "CNSet", function(object,...){

 	constructOligoSetListFrom(object, ...)

 })

 setMethod(BafLrrSetList, "CNSet", function(object,...){

-	constructBafLrrSetListFrom(object, ...)

+  constructBafLrrSetListFrom(object, ...)

 })

@@ -689,41 +689,41 @@ setMethod(BafLrrSetList, "CNSet", function(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)

+  ##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)

 }

@@ -571,7 +571,7 @@ calculateRBafCNSet <- function(object, batch.name, chrom){

 		for(j in seq_along(sampleindex)){

 			bname <- batch.name[j]

 			J <- sampleindex[[j]]

-			res <- crlmm:::calculateRTheta(object=object,

+			res <- calculateRTheta(object=object, # crlmm:::

 						       batch.name=bname,

 						       feature.index=i)

 			k <- match(sampleNames(object)[J], sns)

@@ -1,5 +1,6 @@

 setMethod("posteriorMean", signature(object="CNSet"), function(object) assayDataElement(object, "posteriorMean"))

 setReplaceMethod("posteriorMean", signature(object="CNSet", value="matrix"), function(object, value) assayDataElementReplace(object, "posteriorMean", value))

+##setReplaceMethod("mixtureParams", signature(object="CNSet", value="ANY"), function(object, value) object@mixtureParams <- mixtureParams)

 linearParamElementReplace <- function(obj, elt, value) {

     storage.mode <- storageMode(batchStatistics(obj))

@@ -668,3 +669,87 @@ calculateRBafCNSet <- function(object, batch.name, chrom){

 ##			      lrr=lrr)

 ##		  return(res)

 ##	  })

+

+

+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

@@ -583,6 +583,7 @@ calculateRBafCNSet <- function(object, batch.name, chrom){

 	if(isPackageLoaded("ff")){

 		pkgs <- c("oligoClasses", "ff", "Biobase", "crlmm")

 	} else pkgs <- c("oligoClasses", "Biobase", "crlmm")

+	i <- NULL

 	res <- foreach(i=indexlist, chr=names(indexlist), .packages=pkgs) %dopar% {

 		processByChromosome(i, chr, path)

 	}

@@ -548,7 +548,7 @@ calculateRBafCNSet <- function(object, batch.name, chrom){

 	indexlist <- split(index, chr[index])

 	J <- which(batch(object) %in% batch.name)

 	sns <- sampleNames(object)[J]

-	sampleindex <- split(J, batch(object)[J])

+	sampleindex <- split(J, factor(batch(object)[J], levels=unique(batch(object)[J])))

 	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)

@@ -535,7 +535,11 @@ setMethod("calculateRBaf", signature(object="CNSet"),

 	  })

 calculateRBafCNSet <- function(object, batch.name, chrom){

-	if(missing(batch.name)) batch.name <- batchNames(object)

+	if(missing(batch.name)) {

+		batch.name <- batchNames(object)

+		if("grandMean" %in% batch.name)

+			batch.name <- batch.name[-length(batch.name)]

+	}

 	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)

@@ -549,35 +553,42 @@ calculateRBafCNSet <- function(object, batch.name, chrom){

 	## if ff package is loaded, these will be ff objects

 	chr <- names(indexlist)

 	rlist <- blist <- vector("list", length(indexlist))

-	for(i in seq_along(indexlist)){

-		I <- indexlist[[i]]

-		nr <- length(I)

-		CHR <- names(indexlist)[i]

-		bafname <- paste("baf_chr", CHR, sep="")

-		rname <- paste("lrr_chr", CHR, sep="")

+	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="")

 		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

-		ix <- order(position(object)[I])

-		I <- I[ix]

-		rownames(rmatrix) <- rownames(bmatrix) <- featureNames(object)[I]

+		ix <- order(position(object)[i])

+		i <- i[ix]

+		rownames(rmatrix) <- rownames(bmatrix) <- featureNames(object)[i]

 		for(j in seq_along(sampleindex)){

 			bname <- batch.name[j]

 			J <- sampleindex[[j]]

-			res <- calculateRTheta(object=object,

-					       batch.name=bname,

-					       feature.index=I)

+			res <- crlmm:::calculateRTheta(object=object,

+						       batch.name=bname,

+						       feature.index=i)

 			k <- match(sampleNames(object)[J], sns)

 			bmatrix[, k] <- res[["baf"]]

 			rmatrix[, k] <- res[["lrr"]]

-			##colnames(bmatrix)[J] <- colnames(rmatrix)[J] <- sampleNames(object)[J]

 		}

-		blist[[i]] <- bmatrix

-		rlist[[i]] <- rmatrix

+		list(bmatrix, rmatrix)

+	}

+	## calcualte R BAF by chromosome

+	if(isPackageLoaded("ff")){

+		pkgs <- c("oligoClasses", "ff", "Biobase", "crlmm")

+	} else pkgs <- c("oligoClasses", "Biobase", "crlmm")

+	res <- foreach(i=indexlist, chr=names(indexlist), .packages=pkgs) %dopar% {

+		processByChromosome(i, chr, path)

 	}

-	res <- list(baf=blist,

-		    lrr=rlist)

+	blist <- lapply(res, "[[", 1)

+	rlist <- lapply(res, "[[", 2)

+	res <- list(baf=blist, lrr=rlist)

 	return(res)

 }

@@ -1,63 +1,6 @@

 setMethod("posteriorMean", signature(object="CNSet"), function(object) assayDataElement(object, "posteriorMean"))

 setReplaceMethod("posteriorMean", signature(object="CNSet", value="matrix"), function(object, value) assayDataElementReplace(object, "posteriorMean", value))

-setAs("CNSet", "oligoSnpSet", function(from, to){

-	cnSet2oligoSnpSet(from)

-})

-

-cnSet2oligoSnpSet <- 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 <- integerMatrix(b.r[[1]], 1000)

-	r <- integerMatrix(b.r[[2]], 100)

-##	if(is.lds){

-##		## initialize a big matrix for raw copy number

-##		message("creating an ff object for storing total copy number")

-##		tcn <- initializeBigMatrix(name="total_cn", nrow(object), ncol(object), vmode="double")

-##		for(j in 1:ncol(object)){

-##			tcn[, j] <- totalCopynumber(object, i=row.index, j=j)

-##		}

-##	} else {

-##		if(ncol(object) > 5){

-##			##this can be memory intensive, so we try to be careful

-##			col.index <- splitIndicesByLength(seq(length=ncol(object)), 5)

-##			tcn <- matrix(NA, nrow(object), ncol(object))

-##			dimnames(tcn) <- list(featureNames(object), sampleNames(object))

-##			rows <- 1:nrow(object)

-##			for(i in seq_along(col.index)){

-##				cat(".")

-##				j <- col.index[[i]]

-##				cnSet <- object[, j]

-##				tcn[, j] <- totalCopynumber(cnSet, i=row.index, j=1:ncol(cnSet))

-##				rm(cnSet); gc()

-##			}

-##			cat("\n")

-##		} else {

-##			tcn <- totalCopynumber(object, i=row.index, j=col.index)

-##		}

-##	}

-##	message("Transforming copy number to log2 scale")

-##	tcn[tcn < 0.1] <- 0.1

-##	tcn[tcn > 8] <- 8

-##	log.tcn <- log2(tcn)

-	tmp <- new("oligoSnpSet",

-		   copyNumber=integerMatrix(b.r[[2]], scale=100),

-		   call=calls(object),

-		   callProbability=snpCallProbability(object),

-		   annotation=annotation(object),

-		   featureData=featureData(object),

-		   phenoData=phenoData(object),

-		   experimentData=experimentData(object),

-		   protocolData=protocolData(object),

-		   is.integer=FALSE)

-	tmp <- assayDataElementReplace(tmp, "baf", integerMatrix(b.r[[1]], scale=1000))

-	return(tmp)

-}

-

-

 linearParamElementReplace <- function(obj, elt, value) {

     storage.mode <- storageMode(batchStatistics(obj))

     switch(storage.mode,

@@ -185,7 +128,7 @@ ACN <- function(object, allele, i , j){

 	}

 	## Define batch.index and sample.index

 	if(!missing.j) {

-		batches <- unique(as.character(batch(object))[j])

+		batches <- unique(as.character(batch(object)[j]))

 		##batches <- as.character(batch(object)[j])

 		batch.index <- match(batches, batchNames(object))

 	} else {

@@ -587,78 +530,129 @@ setMethod("xyplot", signature(x="formula", data="CNSet"),

 })

 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]

+	  function(object, batch.name, chrom){

+		  calculateRBafCNSet(object, batch.name, chrom)

+	  })

-		  index.np <- which(is.np)

-		  if(length(index.np) > 0){

-			  a.np <- A(object)[index.np, J]

-			  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)

+calculateRBafCNSet <- function(object, batch.name, chrom){

+	if(missing(batch.name)) batch.name <- batchNames(object)

+	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]

+	sampleindex <- split(J, batch(object)[J])

+	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))

+	for(i in seq_along(indexlist)){

+		I <- indexlist[[i]]

+		nr <- length(I)

+		CHR <- names(indexlist)[i]

+		bafname <- paste("baf_chr", CHR, sep="")

+		rname <- paste("lrr_chr", CHR, sep="")

+		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

+		ix <- order(position(object)[I])

+		I <- I[ix]

+		rownames(rmatrix) <- rownames(bmatrix) <- featureNames(object)[I]

+		for(j in seq_along(sampleindex)){

+			bname <- batch.name[j]

+			J <- sampleindex[[j]]

+			res <- calculateRTheta(object=object,

+					       batch.name=bname,

+					       feature.index=I)

+			k <- match(sampleNames(object)[J], sns)

+			bmatrix[, k] <- res[["baf"]]

+			rmatrix[, k] <- res[["lrr"]]

+			##colnames(bmatrix)[J] <- colnames(rmatrix)[J] <- sampleNames(object)[J]

+		}

+		blist[[i]] <- bmatrix

+		rlist[[i]] <- rmatrix

+	}

+	res <- list(baf=blist,

+		    lrr=rlist)

+	return(res)

+}

-		  dimnames(bf) <- dimnames(lrr) <- dns

-		  res <- list(baf=bf,

-			      lrr=lrr)

-		  return(res)

-	  })

+##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)

+##	  })

@@ -11,6 +11,8 @@ cnSet2oligoSnpSet <- 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 <- integerMatrix(b.r[[1]], 1000)

+	r <- integerMatrix(b.r[[2]], 100)

 ##	if(is.lds){

 ##		## initialize a big matrix for raw copy number

 ##		message("creating an ff object for storing total copy number")

@@ -438,36 +440,36 @@ setMethod("calculatePosteriorMean", signature(object="CNSet"),

 		  return(pm)

 	  })

-		  .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

-		  }

+##.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

+##		  }

 ## for a given copy number, return a named list of bivariate normal prediction regions

 ##   - elements of list are named by genotype

@@ -42,17 +42,16 @@ cnSet2oligoSnpSet <- function(object){

 ##	tcn[tcn > 8] <- 8

 ##	log.tcn <- log2(tcn)

 	tmp <- new("oligoSnpSet",

-		   ##copyNumber=log.tcn,

-		   copyNumber=b.r[[2]],

-		   ##baf=b.r[[1]],

+		   copyNumber=integerMatrix(b.r[[2]], scale=100),

 		   call=calls(object),

 		   callProbability=snpCallProbability(object),

 		   annotation=annotation(object),

 		   featureData=featureData(object),

 		   phenoData=phenoData(object),

 		   experimentData=experimentData(object),

-		   protocolData=protocolData(object))

-	tmp <- assayDataElementReplace(tmp, "baf", b.r[[1]])

+		   protocolData=protocolData(object),

+		   is.integer=FALSE)

+	tmp <- assayDataElementReplace(tmp, "baf", integerMatrix(b.r[[1]], scale=1000))

 	return(tmp)

 }

@@ -1,6 +1,62 @@

 setMethod("posteriorMean", signature(object="CNSet"), function(object) assayDataElement(object, "posteriorMean"))

 setReplaceMethod("posteriorMean", signature(object="CNSet", value="matrix"), function(object, value) assayDataElementReplace(object, "posteriorMean", value))

+setAs("CNSet", "oligoSnpSet", function(from, to){

+	cnSet2oligoSnpSet(from)

+})

+

+cnSet2oligoSnpSet <- 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)

+##	if(is.lds){

+##		## initialize a big matrix for raw copy number

+##		message("creating an ff object for storing total copy number")

+##		tcn <- initializeBigMatrix(name="total_cn", nrow(object), ncol(object), vmode="double")

+##		for(j in 1:ncol(object)){

+##			tcn[, j] <- totalCopynumber(object, i=row.index, j=j)

+##		}

+##	} else {

+##		if(ncol(object) > 5){

+##			##this can be memory intensive, so we try to be careful

+##			col.index <- splitIndicesByLength(seq(length=ncol(object)), 5)

+##			tcn <- matrix(NA, nrow(object), ncol(object))

+##			dimnames(tcn) <- list(featureNames(object), sampleNames(object))

+##			rows <- 1:nrow(object)

+##			for(i in seq_along(col.index)){

+##				cat(".")

+##				j <- col.index[[i]]

+##				cnSet <- object[, j]

+##				tcn[, j] <- totalCopynumber(cnSet, i=row.index, j=1:ncol(cnSet))

+##				rm(cnSet); gc()

+##			}

+##			cat("\n")

+##		} else {

+##			tcn <- totalCopynumber(object, i=row.index, j=col.index)

+##		}

+##	}

+##	message("Transforming copy number to log2 scale")

+##	tcn[tcn < 0.1] <- 0.1

+##	tcn[tcn > 8] <- 8

+##	log.tcn <- log2(tcn)

+	tmp <- new("oligoSnpSet",

+		   ##copyNumber=log.tcn,

+		   copyNumber=b.r[[2]],

+		   ##baf=b.r[[1]],

+		   call=calls(object),

+		   callProbability=snpCallProbability(object),

+		   annotation=annotation(object),

+		   featureData=featureData(object),

+		   phenoData=phenoData(object),

+		   experimentData=experimentData(object),

+		   protocolData=protocolData(object))

+	tmp <- assayDataElementReplace(tmp, "baf", b.r[[1]])

+	return(tmp)

+}

+

+

 linearParamElementReplace <- function(obj, elt, value) {

     storage.mode <- storageMode(batchStatistics(obj))

     switch(storage.mode,

@@ -554,7 +554,8 @@ setMethod("calculateRBaf", signature(object="CNSet"),

 		  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, ] <- a[is.np, ]

+		  b[is.np, ] <- 0L

 		  dns <- dimnames(a)

 		  dimnames(a) <- dimnames(b) <- NULL

 		  obs.theta <- atan2(b, a)*2/pi

@@ -592,9 +592,11 @@ setMethod("calculateRBaf", signature(object="CNSet"),

 		  r.expected[grBB] <- r.bb[grBB]

 		  index.np <- which(is.np)

-		  a.np <- A(object)[index.np, ]

-		  meds <- rowMedians(a.np, na.rm=TRUE)

-		  r.expected[index.np, ] <- matrix(meds, length(index.np), ncol(a.np))

+		  if(length(index.np) > 0){

+			  a.np <- A(object)[index.np, J]

+			  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

@@ -576,9 +576,9 @@ setMethod("calculateRBaf", signature(object="CNSet"),

 		  bf[grBB] <- 1

 		  r.expected <- matrix(NA, nrow(object), ncol(a))

-		  r.aa <- matrix(RTheta.aa[, "R"], nrow(object), ncol(object), byrow=FALSE)

-		  r.ab <- matrix(RTheta.ab[, "R"], nrow(object), ncol(object), byrow=FALSE)

-		  r.bb <- matrix(RTheta.bb[, "R"], nrow(object), ncol(object), byrow=FALSE)

+		  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

@@ -519,6 +519,7 @@ setMethod("xyplotcrlmm", signature(x="formula", data="CNSet", predictRegion="lis

 		  ##predictRegion is an argument of ABpanel

 		  xyplot(x, df, predictRegion=predictRegion, ...)

 	  })

+

 setMethod("xyplot", signature(x="formula", data="CNSet"),

 	  function(x, data, ...){

 		  if("predictRegion" %in% names(list(...))){

@@ -552,6 +553,8 @@ setMethod("calculateRBaf", signature(object="CNSet"),

 		  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, ]

 		  dns <- dimnames(a)

 		  dimnames(a) <- dimnames(b) <- NULL

 		  obs.theta <- atan2(b, a)*2/pi

@@ -588,7 +591,7 @@ setMethod("calculateRBaf", signature(object="CNSet"),

 		  r.expected[lessAA] <- r.aa[lessAA]

 		  r.expected[grBB] <- r.bb[grBB]

-		  index.np <- which(!isSnp(object))

+		  index.np <- which(is.np)

 		  a.np <- A(object)[index.np, ]

 		  meds <- rowMedians(a.np, na.rm=TRUE)

 		  r.expected[index.np, ] <- matrix(meds, length(index.np), ncol(a.np))

@@ -530,68 +530,72 @@ setMethod("xyplot", signature(x="formula", data="CNSet"),

 setMethod("calculateRBaf", signature(object="CNSet"),

 	  function(object, batch.name){

-	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)

+		  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)

+		  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)

+		  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

-	dns <- dimnames(a)

-	dimnames(a) <- dimnames(b) <- NULL

-	obs.theta <- atan2(b, a)*2/pi

+		  a <- A(object)[, J, drop=FALSE]

+		  b <- B(object)[, J, drop=FALSE] ## NA's for b where nonpolymorphic

+		  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

+		  lessAA <- obs.theta < theta.aa

+		  lessAB <- obs.theta < theta.ab

+		  lessBB <- obs.theta < theta.bb

+		  grAA <- !lessAA

+		  grAB <- !lessAB