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

 * Fixed downstream vignette to account for the SnpSet object being returned by crlmm

 * Fixed some minor status messages

+

+2009-04-09 R Scharpf - committed version 1.0.81

+

+* added a skeleton of a copy number vignette for the illumina platform (illumina_copynumber.Rnw) to inst/scripts

@@ -1,7 +1,7 @@

 Package: crlmm

 Type: Package

 Title: Genotype Calling (CRLMM) and Copy Number Analysis tool for Affymetrix SNP 5.0 and 6.0 and Illumina arrays.

-Version: 1.0.80

+Version: 1.0.81

 Date: 2008-12-30

 Author: Rafael A Irizarry, Benilton S Carvalho <bcarvalh@jhsph.edu>, Robert Scharpf <rscharpf@jhsph.edu>, Matt Ritchie <mritchie@wehi.EDU.AU>

 Maintainer: Benilton S Carvalho <bcarvalh@jhsph.edu>, Robert Scharpf <rscharpf@jhsph.edu>, Matt Ritchie <mritchie@wehi.EDU.AU>

@@ -188,7 +188,8 @@ goodSnps <- function(phi.thr, envir, fewAA=20, fewBB=20){

 instantiateObjects <- function(calls, conf, NP, plate, envir, chrom, A, B,

 			       gender, SNRmin=5, SNR,

-                               pkgname){

+                               pkgname,

+			       locusSet){

 	pkgname <- paste(pkgname, "Crlmm", sep="")

 	envir[["chrom"]] <- chrom

 	CHR_INDEX <- paste(chrom, "index", sep="")

@@ -615,7 +616,7 @@ oneBatch <- function(plateIndex,

 	CHR <- envir[["chrom"]]

 	if(bias.adj){

 		nuA <- envir[["nuA"]]

-		if(all(is.na(nuA))) {

+		if(all(is.na(nuA))){

 			message("Background and signal coefficients have not yet been estimated -- can not do bias correction yet")

 			stop("Must run computeCopynumber a second time with bias.adj=TRUE to do the adjustment")

 		}

@@ -1021,7 +1022,7 @@ biasAdj <- function(plateIndex, envir, priorProb, PROP=0.75){

 	hemDel <- hemDel/total

 	norm <- norm/total

 	amp <- amp/total

-	envir[["posteriorProb"]] <- list(hemDel=hemDel, norm=norm, amp=amp) 

+	envir[["posteriorProb"]] <- list(hemDel=hemDel, norm=norm, amp=amp)

 	tmp <- array(NA, dim=c(nrow(A), ncol(A), 4))

 	tmp[, , 1] <- homDel

 	tmp[, , 2] <- hemDel

@@ -1047,18 +1048,21 @@ biasAdj <- function(plateIndex, envir, priorProb, PROP=0.75){

 	##whether more are up or down

 	if(CHR != 23){

 		moreup <- rowSums(tmp2 > 3) > rowSums(tmp2 < 3) ##3 is normal

-		##notUp <-  tmp2[ii & moreup, ] <= 3

-		##notDown <- tmp2[ii & !moreup, ] >= 3

-		

+		down <-  tmp2[ii & moreup, ] <= 3

+		up <- tmp2[ii & !moreup, ] >= 3

 		##What if we just keep X% of the ones with the highest

 		##posterior probability for normal

-		prNorm <- tmp[, , 3]

-		rankNorm <- t(apply(prNorm, 1, order, decreasing=TRUE))

-		percentage <- round(0.75*ncol(prNorm),0)

-		NORM <- rankNorm <= percentage

-##		NORM <- matrix(NA, nrow(A), ncol(A))

-##		NORM[ii & moreup, ] <- notUp

-##		NORM[ii & !moreup, ] <- notDown

+##		prNorm <- tmp[, , 3]

+##		rankNorm <- t(apply(prNorm, 1, order, decreasing=TRUE))

+

+		rankUP <- t(apply(tmp[moreup & ii, , 4]/tmp[moreup & ii, , 3], 1, order, decreasing=TRUE))

+		rankDown <- t(apply(tmp[!moreup & ii, , 2]/tmp[!moreup & ii, , 3], 1, order, decreasing=TRUE))		

+		maxNumberToDrop <- round(0.25*ncol(tmp2),0)

+##		NORM <- rankNorm <= percentage

+		NORM <- matrix(TRUE, nrow(A), ncol(A))

+		NORM[ii & moreup, ] <- rankUP >= maxNumberToDrop

+		NORM[ii & !moreup, ] <- rankDown >= maxNumberToDrop

+##		NORM[ii & !moreup, ] <- up

 		##Define NORM so that we can iterate this step

 		##NA's in the previous iteration (normal) will be propogated

 		normal <- NORM*normal

new file mode 100755

@@ -0,0 +1,80 @@

+<<>>=

+# start R from /home/bst/other/mritchie/R/R-300309/bin/R

+# as this has current crlmm and region package installed

+# if you want to ue your own R, the current region package is at

+# /thumper/ctsa/snpmicroarray/illumina/crlmm/370k/human370v1cCrlmm/

+library(Biobase)

+library(crlmm)

+setwd("/thumper/ctsa/snpmicroarray/illumina/IDATS/370k")

+@ 

+

+<<readIdat>>=

+samplesheet5 = read.csv("HumanHap370Duo_Sample_Map.csv", header=TRUE, as.is=TRUE)[-c(28:46,61:75,78:79),]

+if(!exists("RG")){

+	## remove samples which I don't have .idats for

+	## subset further to allow for quicker testing

+	##samplesheet5 = samplesheet5[1:5,]

+	## read in the .idats

+	RG <- readIdatFiles(samplesheet5, arrayInfoColNames=list(barcode=NULL, position="SentrixPosition"), saveDate=TRUE)

+}

+@ 

+

+Get dates for each array - scan date is likely to be most meaningful

+

+<<dates>>=

+pd <- pData(RG)

+scandatetime = strptime(as.character(pd$ScanDate), "%m/%d/%Y %H:%M:%S %p")

+decodedatetime = strptime(as.character(pd$DecodeDate), "%m/%d/%Y %H:%M:%S %p")

+table(format(scandatetime, "%d %b %Y"))

+table(format(scandatetime, "%b"))

+@ 

+

+Note: the code below is run for testing purposes only. None of these

+functions are exported, so would not routinely be run directly by the

+user.  A typical analysis would involve runnning readIdatFiles()

+followed by crlmmIllumina()

+

+<<crlmm>>=

+if(!exists("res.rda")){

+	outdir <- "/thumper/ctsa/snpmicroarray/rs/data/hapmap/illumina/HumanCNV370-Duo"	

+	if(!file.exists(file.path(outdir, "res.rda"))){

+		crlmmOut <- crlmmIllumina(RG=RG, cdfName="human370v1c", sns=pData(RG)$ID, returnParams=TRUE, save.it=TRUE, intensityFile=file.path(outdir, "res.rda"))

+		save(crlmmOut, file=file.path(outdir, "crlmmOut.rda"))				

+	} else{

+		message("Loading...")		

+		load(file.path(outdir, "res.rda"))

+		load(file.path(outdir, "crlmmOut.rda"))		

+	}

+}

+@ 

+

+<<SNR>>=

+hist(crlmmOut$SNR) ##approx. 5-fold higher than what we see in Affy!

+@ 

+

+<<>>=

+##cnAB was assigned to the global environment

+NP = (cnAB$A+cnAB$B)/2 # average normalized A and B intensities

+colnames(NP) <- colnames(crlmm:::calls(crlmmOut))

+cnenv = new.env()

+chr22 <- computeCopynumber(chrom=22, A=res[["A"]], B=res[["B"]], 

+			   calls=crlmm:::calls(crlmmOut),

+			   conf=confs(crlmmOut), 

+			   plate=rep(1, ncol(NP)),

+			   NP=NP, 

+			   cdfName="human370v1c", 

+			   SNR=res[["SNR"]], 

+			   envir=cnenv)

+@ 

+

+<<>>=

+copyA=cnenv[["CA"]]

+copyB=cnenv[["CB"]]

+copyT=(copyA + copyB)/100

+copyT[copyT > 6] <- 6

+copyT[copyT < 0] <- 0

+plot(copyT[, 1], pch=".")

+@ 

+

+

+