@@ -1,4 +1,4 @@

-2009-03-01  Benilton Carvalho <bcarvalh@jhsph.edu>

+2009-03-01  Benilton Carvalho <bcarvalh@jhsph.edu> - committed version 1.0.40

 * Added TODO file

@@ -9,3 +9,11 @@

 * Files with suffix XYZ-functions.R should have *functions* for methodology XYZ

 * Files with suffix XYZ-methods.R should have *methods* for methodology XYZ

+

+* Exported assayDataNew from Biobase

+

+* Added inst/scripts to store vignettes that can't be built by BioC

+due to the fact that they depend on external data not available as

+BioC data packages

+

+* Added crlmmSet class and calls/confs methods

@@ -1,14 +1,14 @@

 Package: crlmm

 Type: Package

 Title: Genotype Calling via CRLMM Algorithm

-Version: 1.0.38

+Version: 1.0.40

 Date: 2008-12-28

 Author: Rafael A Irizarry

 Maintainer: Benilton S Carvalho <bcarvalh@jhsph.edu>

 Description: Faster implementation of CRLMM specific to SNP 5.0 and 6.0 arrays.

-Imports: affyio, preprocessCore, utils, stats, genefilter, splines, Biobase

 License: Artistic-2.0

-Suggests: hapmapsnp5, genomewidesnp5Crlmm, genomewidesnp6Crlmm, methods

+Imports: affyio, preprocessCore, utils, stats, genefilter, splines, Biobase

+Suggests: hapmapsnp5, genomewidesnp5Crlmm, genomewidesnp6Crlmm, methods, GGdata, snpMatrix

 Collate: AllClasses.R

          crlmm-methods.R

 	 cnrma-functions.R

@@ -2,10 +2,10 @@ useDynLib("crlmm", .registration=TRUE)

 export("crlmm", "list.celfiles", "computeCnBatch")

 exportMethods("list2crlmmSet", "calls", "confs")

 exportClasses("crlmmSet")

+import(Biobase)

 importFrom(affyio, read.celfile.header, read.celfile)

 importFrom(preprocessCore, normalize.quantiles.use.target)

 importFrom(utils, data, packageDescription, setTxtProgressBar, txtProgressBar)

 importFrom(stats, coef, cov, dnorm, kmeans, lm, mad, median, quantile, sd)

-importFrom(Biobase, rowMedians)

 importFrom(genefilter, rowSds)

@@ -3,11 +3,27 @@

 #####################################

 - Decide on output format (BC vote: eSet-like)

 - Helper to convert to snpMatrix

-- Vignette with non-trivial use of crlmm (use Vince's)

-- Add RS ids

+- Add RS ids to annotation packages

 - Allele plots

 - M v S plots

+

 #####################################

 ### FOR CNRMA

 #####################################

+- fix the following 2 items

+* checking R code for possible problems ... NOTE

+cnrma: no visible global function definition for ‘getCnvFid’

+cnrma: no visible binding for global variable ‘reference’

+coefs: no visible binding for global variable ‘nu’

+coefs: no visible binding for global variable ‘phi’

+coefs: no visible binding for global variable ‘nu.se’

+coefs: no visible binding for global variable ‘phi.se’

+

+* checking for unstated dependencies in R code ... WARNING

+'library' or 'require' calls not declared from:

+  affyio splines Biobase genefilter

+

+* checking for missing documentation entries ... WARNING

+Undocumented code objects:

+  computeCnBatch

new file mode 100644

@@ -0,0 +1,10 @@

+all: copynumber crlmmDownstream clean

+

+copynumber: copynumber.tex

+	cp -p ../scripts/copynumber.pdf .

+

+crlmmDownstream: crlmmDownstream.tex

+	cp -p ../scripts/crlmmDownstream.pdf .

+

+clean:

+	-$(RM) -f *.out *.bbl *.log *.aux *.blg *.brf *.toc *.tex

@@ -1,378 +1,11 @@

-%\VignetteIndexEntry{crlmm copy number Vignette}

+%\VignetteIndexEntry{crlmm Vignette - Copy Number}

 %\VignetteDepends{crlmm, genomewidesnp6Crlmm}

 %\VignetteKeywords{crlmm, SNP 6}

 %\VignettePackage{crlmm}

-\documentclass{article}

-\usepackage{graphicx}

-\newcommand{\Rfunction}[1]{{\texttt{#1}}}

-\newcommand{\Rmethod}[1]{{\texttt{#1}}}

-\newcommand{\Rcode}[1]{{\texttt{#1}}}

-\newcommand{\Robject}[1]{{\texttt{#1}}}

-\newcommand{\Rpackage}[1]{{\textsf{#1}}}

-\newcommand{\Rclass}[1]{{\textit{#1}}}

-\newcommand{\oligo}{\Rpackage{oligo }}

-\newcommand{\R}{\textsf{R}}

+%%% The real .Rnw is under inst/scripts.

+%%% It cannot be added here because it depends on external

+%%% data not yet available as a BioC package

+\documentclass{article}

 \begin{document}

-\title{Estimating copy number for Affymetrix 6.0  with the crlmm Package}

-\date{February, 2009}

-\author{Rob Scharpf}

-\maketitle

-

-<<setup, echo=FALSE, results=hide>>=

-options(width=60)

-options(continue=" ")

-options(prompt="R> ")

-@ 

-

-\section{Estimating copy number}

-

-General requirements: in addition to the \R{} packages indicated

-below, processing a large number of samples requires a high-end

-computer with a large amount of RAM. Currently, the maximum number of

-samples that can be genotyped at once is approximately 2000 and this

-requires roughly 32G of RAM.

-

-\paragraph{Suggested work flow.} I typically submit code for

-preprocessing and genotyping as a batch job to a computing cluster.

-When the preprocessing is complete, I pick a chromosome and work

-interactively from the cluster to calculate copy number and make some

-of the suggested diagnostic plots.

-

-\subsection{Preprocess and genotype the samples}

-

-First, load the required libraries.

-

-<<requiredPackages>>=

-library(crlmm)

-library(genomewidesnp6Crlmm)

-library(ellipse)

-@ 

-

-Specify an output directory and provide the complete path for the CEL

-file names.

-

-<<>>=

-outdir <- "/thumper/ctsa/snpmicroarray/rs/data/gain/bipolar/GRU-EA"

-datadir <- "/thumper/ctsa/snpmicroarray/GAIN/Bipolar/GRU-EA"

-fns <- list.files(datadir, pattern=".CEL", full.names=TRUE)

-@ 

-

-Genotype the samples with CRLMM. Submit this job to the cluster and save the

-output to \Robject{outdir}.

-

-<<genotype, eval=FALSE>>=

-res2 <- crlmm(filenames=fns,  save.it=TRUE, intensityFile=file.path(outdir, "intensities.rda"))

-save(res2, file=file.path(outdir, "res2.rda"))

-@ 

-

-Quantile normalize the nonpolymorphic probes and save the output:

-

-<<quantileNormalizeCnProbes, eval=FALSE>>=

-res3 <- cnrma(fns)

-save(res3, file=file.path(outdir, "res3.rda"))

-@ 

-

-\subsection{Copy number}

-

-Copy number can be assessed one chromosome at a time.  Here we specify

-chromosome 15 and and load a list of indices used to subset the files

-saved in the previous section. The first element in the list

-correspond to indices of polymorphic probes on chromosome 15; the

-second element corresponds to indices of nonpolymorphic probes on

-chromosome 15.

-

-<<chromosomeIndex>>=

-CHR <- 15

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

-data(list=CHR_INDEX, package="genomewidesnp6Crlmm")

-str(index)

-@ 

-

-Next we load 3 files that were saved from the preprocessing step and

-then subset these lists using the above indices to extract the

-preprocessed intensities and genotypes needed for estimating copy

-number.  Specifically, we require 6 items:

-

-\begin{itemize}

-  \item quantile-normalized A intensities (I1 x J)

-  \item quantile-normalized B intensities (I1 x J)

-  \item quantile-normalized intensities from nonpolymorphic (NP) probes (I2 x J)

-  \item genotype calls (I1 x J)

-  \item confidence scores of the genotype calls  (I1 x J)

-  \item signal to noise ratio (SNR) of the samples (J)

-  \end{itemize}

-  

-  These items are extracted as follows:

-

-<<eval=FALSE>>=

-load(file.path(outdir, "intensities.rda"))

-load(file.path(outdir, "res2.rda"))

-load(file.path(outdir, "res3.rda"))

-@ 

-

-Depending on the number of samples, the above objects can be large and

-take a while to load.  For the figures in this vignette, I am loading

-only 5000 SNP-level summaries from chromosome 22 (28MB) and the

-quantile-normalized intensities for 5000 nonpolymorphic probes on

-chromosome 22 (8 MB total).

-

-<<loadTmpFiles, eval=FALSE, echo=FALSE>>=

-load("/thumper/ctsa/snpmicroarray/rs/data/gain/bipolar/tmp.rda")

-assign("res", tmp, envir=.GlobalEnv)

-rm(tmp); gc()

-load("/thumper/ctsa/snpmicroarray/rs/data/gain/bipolar/tmpcn.rda")

-A <- res$A

-B <- res$B

-calls <- res$calls

-conf <- res$conf

-SNR <- res$SNR

-NP <- res3$NP

-@ 

-

-<<snpAndCnSummaries, eval=FALSE>>=

-A <- res$A

-B <- res$B

-calls <- res2$calls

-conf <- res2$conf

-SNR <- res2$SNR

-NP <- res3$NP

-rm(res, res2, res3)

-gc()

-@ 

-

-Make a histogram of the signal to noise ratio for these samples:

-

-<<plotSnr, fig=TRUE>>=

-hist(SNR, xlab="SNR", main="")

-@ 

-

-We suggest excluding samples with a signal to noise ratio less than 5.

-We then extract the probes (rows) that are on chromosome 15 and the

-samples (columns) that have a suitable signal to noise ratio.

-

-<<subsetCrlmmOutput, eval=FALSE>>=

-A <- A[index[[1]], SNR > 5]

-B <- B[index[[1]], SNR > 5]

-calls <- calls[index[[1]], SNR > 5]

-conf <- conf[index[[1]], SNR > 5]

-NP <- NP[index[[2]], SNR > 5]

-rm(res, res2, res3); gc()

-@ 

-

-We want to adjust for batch effects.  In our experience, the chemistry

-plate is a good surrogate for batch, although this should be assessed

-on a study by study basis.  For samples processed by Broad, the plate

-is indicated by a capitalized five-letter word and is part of the

-sample name.  It is important that the plate (batch) variable is

-ordered the same as filenames in the preprocessed data.  For instance,

-to indicate plate in this example one can use the command

-

-<<specifyBatch>>=

-sns <- colnames(calls)

-sns[1]

-plates <- sapply(sns, function(x) strsplit(x, "_")[[1]][2])

-table(plates)

-@ 

-

-We are now ready to estimate the copy number for each batch.  In the

-current version of this package, one specifies an environment to which

-intermediate \R{} objects for copy number estimation are

-stored. Allele-specific estimates of copy number are also stored in

-this environment.

-

-<<copyNumberByBatch, eval=FALSE>>=

-e1 <- new.env()

-computeCnBatch(A=A,

-	       B=B,

-	       calls=calls,

-	       conf=conf,

-	       NP=NP,

-	       plate=plates,

-	       envir=e1, chrom=CHR, DF.PRIOR=75)

-@ 

-

-The DF.PRIOR indicates how much we will shrink SNP-specific estimates

-of the variance and correlation.

-

-<<accessingEstimates>>=

-copyA <- get("CA", e1)

-copyB <- get("CB", e1)

-copyT <- (copyA+copyB)/100

-copyT[copyT < 0] <- 0

-copyT[copyT > 6] <- 6

-@ 

-

-\section{Suggested plots}

-

-\paragraph{One sample at a time}

-

-<<annotation>>=

-data(snpProbes, package="genomewidesnp6Crlmm")

-data(cnProbes, package="genomewidesnp6Crlmm")

-position <- snpProbes[match(rownames(calls), rownames(snpProbes)), "position"]

-@ 

-

-To plot physical position versus copy number for the first sample:

-

-<<oneSample, fig=TRUE, width=8, height=4, include=FALSE>>=

-require(SNPchip)

-par(las=1)

-plotCytoband(as.character(CHR), ylim=c(0,7), cytoband.ycoords=c(6.5,7), label.cytoband=FALSE, main="Chr 15")

-points(position, copyT[, 1], pch=".", cex=2, , xaxt="n", col="grey70")

-axis(1, at=pretty(range(position)), labels=pretty(range(position))/1e6)

-axis(2, at=0:5, labels=0:5)

-mtext("position (Mb)", 1, line=2)

-mtext(expression(C[A] + C[B]), 2, line=2, las=3)

-@ 

-

-\begin{figure}

-  \includegraphics[width=0.9\textwidth]{copynumber-oneSample}

-  \caption{Total copy number (y-axis) for chromosome 22 plotted

-    against physical position (x-axis) for one sample.  }

-\end{figure}

-

-\paragraph{One SNP at a time}

-

-<<intermediateFiles>>=

-tau2A <- get("tau2A", e1)

-tau2B <- get("tau2B", e1)

-sig2A <- get("sig2A", e1)

-sig2B <- get("sig2B", e1)

-nuA <- get("nuA", e1)

-phiA <- get("phiA", e1)

-nuB <- get("nuB", e1)

-phiB <- get("phiB", e1)

-corr <- get("corr", e1)

-corrA.BB <- get("corrA.BB", e1)

-corrB.AA <- get("corrA.BB", e1)

-A <- get("A", e1)

-B <- get("B", e1)

-@ 

-

-Here, we plot the prediction regions for total copy number 2 and 3 for

-the first plate. Black plotting symbols are estimates from the first

-plate; light grey are points from other plates. (You could also draw

-prediction regions for 0-4 copies, but it gets crowded).  Notice that

-there is little evidence of a plate effect for this SNP.

-

-<<predictionRegion, fig=TRUE, width=8, height=8, include=FALSE>>=

-par(las=1)

-p <- 1 ##indicates plate

-J <- grep(unique(plates)[p], sns) ##sample indices for this plate

-ylim <- c(6.5,13)

-I <- which(phiA > 10 & phiB > 10)

-i <- I[1]

-plot(log2(A[i, ]), log2(B[i, ]), pch=as.character(calls[i, ]), col="grey60", cex=0.9, ylim=ylim, xlim=ylim, xlab="A", ylab="B")

-points(log2(A[i, J]), log2(B[i, J]), col="black", pch=as.character(calls[i, J]))

-for(CT in 2:3){

-	if(CT == 2) ellipse.col <- "brown" else ellipse.col <- "purple"

-	for(CA in 0:CT){

-		CB <- CT-CA

-		A.scale <- sqrt(tau2A[i, p]*(CA==0) + sig2A[i, p]*(CA > 0))

-		B.scale <- sqrt(tau2B[i, p]*(CB==0) + sig2B[i, p]*(CB > 0))

-		scale <- c(A.scale, B.scale)

-		if(CA == 0 & CB > 0) rho <- corrA.BB[i, p]

-		if(CA > 0 & CB == 0) rho <- corrB.AA[i, p]

-		if(CA > 0 & CB > 0) rho <- corr[i, p]		

-		lines(ellipse(x=rho, centre=c(log2(nuA[i, p]+CA*phiA[i, p]), log2(nuB[i, p]+CB*phiB[i, p])),

-			      scale=scale), col=ellipse.col, lwd=2)

-	}

-}

-legend("topright", lwd=3, col=c("black", "purple"), legend=c("2 copies", "3 copies"), bty="n")

-@ 

-

-\begin{figure}

-  \includegraphics[width=0.9\textwidth]{copynumber-predictionRegion}

-  \caption{Prediction regions for copy number 2 and 3 for one SNP.

-    The plate effects are negligible for this SNP and we only plot the

-    prediction region for the first plate.}

-\end{figure}

-

-Here is one way to identify a SNP with a large plate effect -- look

-for shifts in the prediction regions (here I'm lucking for large

-shifts in the A direction). 

-

-<<shifts>>=

-shiftA <- nuA+2*phiA 

-maxA <- apply(shiftA, 1, max, na.rm=T)

-minA <- apply(shiftA, 1, min, na.rm=T)

-d <- maxA - minA

-hist(d)

-index <- which(d > 3000)[2]

-plate1 <- which(shiftA[index, ] == maxA[index])

-plate2 <- which(shiftA[index, ] == minA[index])

-@ 

-

-Now plot the predictions for the two plates.  All the ellipses are

-prediction regions for copy number 2.  Note that while I looked for

-shifts in the A direction, the shift often occurs in both the B and A

-dimensions.

-

-<<plateEffect, fig=TRUE, width=8, height=8, include=FALSE>>=

-par(las=1)

-i <- index

-J1 <- grep(unique(plates)[plate1], sns) ##sample indices for this plate

-J2 <- grep(unique(plates)[plate2], sns) ##sample indices for this plate

-ylim <- c(6.5,13)

-plot(log2(A[i, ]), log2(B[i, ]), pch=as.character(calls[i, ]), col="grey60", cex=0.9, ylim=ylim, xlim=ylim, xlab="A", ylab="B")

-points(log2(A[i, J1]), log2(B[i, J1]), col="brown", pch=as.character(calls[i, J1]))

-points(log2(A[i, J2]), log2(B[i, J2]), col="blue", pch=as.character(calls[i, J2]))

-p <- plate1

-CT <- 2

-ellipse.col <- "brown" 

-for(CA in 0:CT){

-	CB <- CT-CA

-	A.scale <- sqrt(tau2A[i, p]*(CA==0) + sig2A[i, p]*(CA > 0))

-	B.scale <- sqrt(tau2B[i, p]*(CB==0) + sig2B[i, p]*(CB > 0))

-	scale <- c(A.scale, B.scale)

-	lines(ellipse(x=corr[i, p]*(CA>0)*(CB>0), centre=c(log2(nuA[i, p]+CA*phiA[i, p]), log2(nuB[i, p]+CB*phiB[i, p])),

-		      scale=scale), col=ellipse.col, lwd=2)

-}

-ellipse.col <- "blue" 

-p <- plate2

-for(CA in 0:CT){

-	CB <- CT-CA

-	A.scale <- sqrt(tau2A[i, p]*(CA==0) + sig2A[i, p]*(CA > 0))

-	B.scale <- sqrt(tau2B[i, p]*(CB==0) + sig2B[i, p]*(CB > 0))

-	scale <- c(A.scale, B.scale)

-	lines(ellipse(x=corr[i, p]*(CA>0)*(CB>0), centre=c(log2(nuA[i, p]+CA*phiA[i, p]), log2(nuB[i, p]+CB*phiB[i, p])),

-		      scale=scale), col=ellipse.col, lwd=2)

-}

-legend("topright", lwd=3, col=c("brown", "blue"), legend=c("2 copies, plate 1", "2 copies, plate 2 "), bty="n")

-@ 

-

-\begin{figure}

-  \includegraphics[width=0.9\textwidth]{copynumber-plateEffect}

-  \caption{The prediction regions for copy number 2 shift when there

-    is a large plate effect.}

-\end{figure}

-

-Alternatively, loop through a large number of SNPs to see how these

-regions move

-<<codeForLoop, eval=FALSE>>=

-par(las=1, pty="s", ask=TRUE)

-for(i in I[1:50]){

-	plot(log2(A[i, ]), log2(B[i, ]), pch=as.character(calls[i, ]), col=col[DS], cex=0.9, ylim=ylim, xlim=ylim)

-	for(CT in 2:3){

-		if(CT == 2) ellipse.col <- "brown" else ellipse.col <- "purple"

-		for(CA in 0:CT){

-			CB <- CT-CA

-			A.scale <- sqrt(tau2A[i, p]*(CA==0) + sig2A[i, p]*(CA > 0))

-			B.scale <- sqrt(tau2B[i, p]*(CB==0) + sig2B[i, p]*(CB > 0))

-			scale <- c(A.scale, B.scale)

-			lines(ellipse(x=corr[i, p]*(CA>0)*(CB>0), centre=c(log2(nuA[i, p]+CA*phiA[i, p]), log2(nuB[i, p]+CB*phiB[i, p])),

-				      scale=scale), col=ellipse.col, lwd=2)

-		}

-	}

-	legend("topright", lwd=3, col=col, legend=c("3 copies", "2 copies"), bty="n")

-}

-@ 

-

-\section{Session information}

-<<>>=

-sessionInfo()

-@ 

-

-

 \end{document}

new file mode 100644

@@ -0,0 +1,10 @@

+%\VignetteIndexEntry{crlmm Vignette - Downstream Analysis}

+%\VignetteKeywords{genotype, crlmm, SNP 5, SNP 6}

+%\VignettePackage{crlmm}

+

+%%% The real .Rnw is under inst/scripts.

+%%% It cannot be added here because it depends on external

+%%% data not yet available as a BioC package

+\documentclass[12pt]{article}

+\begin{document}

+\end{document}

new file mode 100644

Binary files /dev/null and b/inst/doc/crlmmDownstream.pdf differ

similarity index 98%

rename from inst/doc/crlmm.Rnw

rename to inst/doc/genotyping.Rnw

@@ -1,4 +1,4 @@

-%\VignetteIndexEntry{crlmm Vignette}

+%\VignetteIndexEntry{crlmm Vignette - Genotyping}

 %\VignetteDepends{crlmm, hapmapsnp5, genomewidesnp5Crlmm}

 %\VignetteKeywords{genotype, crlmm, SNP 5, SNP 6}

 %\VignettePackage{crlmm}

similarity index 100%

rename from inst/doc/crlmm.pdf

rename to inst/doc/genotyping.pdf

new file mode 100644

@@ -0,0 +1,378 @@

+%\VignetteIndexEntry{crlmm Vignette - Copy Number}

+%\VignetteDepends{crlmm, genomewidesnp6Crlmm}

+%\VignetteKeywords{crlmm, SNP 6}

+%\VignettePackage{crlmm}

+\documentclass{article}

+\usepackage{graphicx}

+\newcommand{\Rfunction}[1]{{\texttt{#1}}}

+\newcommand{\Rmethod}[1]{{\texttt{#1}}}

+\newcommand{\Rcode}[1]{{\texttt{#1}}}

+\newcommand{\Robject}[1]{{\texttt{#1}}}

+\newcommand{\Rpackage}[1]{{\textsf{#1}}}

+\newcommand{\Rclass}[1]{{\textit{#1}}}

+\newcommand{\oligo}{\Rpackage{oligo }}

+\newcommand{\R}{\textsf{R}}

+

+\begin{document}

+\title{Estimating copy number for Affymetrix 6.0  with the crlmm Package}

+\date{February, 2009}

+\author{Rob Scharpf}

+\maketitle

+

+<<setup, echo=FALSE, results=hide>>=

+options(width=60)

+options(continue=" ")

+options(prompt="R> ")

+@ 

+

+\section{Estimating copy number}

+

+General requirements: in addition to the \R{} packages indicated

+below, processing a large number of samples requires a high-end

+computer with a large amount of RAM. Currently, the maximum number of

+samples that can be genotyped at once is approximately 2000 and this

+requires roughly 32G of RAM.

+

+\paragraph{Suggested work flow.} I typically submit code for

+preprocessing and genotyping as a batch job to a computing cluster.

+When the preprocessing is complete, I pick a chromosome and work

+interactively from the cluster to calculate copy number and make some

+of the suggested diagnostic plots.

+

+\subsection{Preprocess and genotype the samples}

+

+First, load the required libraries.

+

+<<requiredPackages>>=

+library(crlmm)

+library(genomewidesnp6Crlmm)

+library(ellipse)

+@ 

+

+Specify an output directory and provide the complete path for the CEL

+file names.

+

+<<>>=

+outdir <- "/thumper/ctsa/snpmicroarray/rs/data/gain/bipolar/GRU-EA"

+datadir <- "/thumper/ctsa/snpmicroarray/GAIN/Bipolar/GRU-EA"

+fns <- list.files(datadir, pattern=".CEL", full.names=TRUE)

+@ 

+

+Genotype the samples with CRLMM. Submit this job to the cluster and save the

+output to \Robject{outdir}.

+

+<<genotype, eval=FALSE>>=

+res2 <- crlmm(filenames=fns,  save.it=TRUE, intensityFile=file.path(outdir, "intensities.rda"))

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

+@ 

+

+Quantile normalize the nonpolymorphic probes and save the output:

+

+<<quantileNormalizeCnProbes, eval=FALSE>>=

+res3 <- cnrma(fns)

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

+@ 

+

+\subsection{Copy number}

+

+Copy number can be assessed one chromosome at a time.  Here we specify

+chromosome 15 and and load a list of indices used to subset the files

+saved in the previous section. The first element in the list

+correspond to indices of polymorphic probes on chromosome 15; the

+second element corresponds to indices of nonpolymorphic probes on

+chromosome 15.

+

+<<chromosomeIndex>>=

+CHR <- 15

+CHR_INDEX <- paste(CHR, "index", sep="")

+data(list=CHR_INDEX, package="genomewidesnp6Crlmm")

+str(index)

+@ 

+

+Next we load 3 files that were saved from the preprocessing step and

+then subset these lists using the above indices to extract the

+preprocessed intensities and genotypes needed for estimating copy

+number.  Specifically, we require 6 items:

+

+\begin{itemize}

+  \item quantile-normalized A intensities (I1 x J)

+  \item quantile-normalized B intensities (I1 x J)

+  \item quantile-normalized intensities from nonpolymorphic (NP) probes (I2 x J)

+  \item genotype calls (I1 x J)

+  \item confidence scores of the genotype calls  (I1 x J)

+  \item signal to noise ratio (SNR) of the samples (J)

+  \end{itemize}

+  

+  These items are extracted as follows:

+

+<<eval=FALSE>>=

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

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

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

+@ 

+

+Depending on the number of samples, the above objects can be large and

+take a while to load.  For the figures in this vignette, I am loading

+only 5000 SNP-level summaries from chromosome 22 (28MB) and the

+quantile-normalized intensities for 5000 nonpolymorphic probes on

+chromosome 22 (8 MB total).

+

+<<loadTmpFiles, eval=FALSE, echo=FALSE>>=

+load("/thumper/ctsa/snpmicroarray/rs/data/gain/bipolar/tmp.rda")

+assign("res", tmp, envir=.GlobalEnv)

+rm(tmp); gc()

+load("/thumper/ctsa/snpmicroarray/rs/data/gain/bipolar/tmpcn.rda")

+A <- res$A

+B <- res$B

+calls <- res$calls

+conf <- res$conf

+SNR <- res$SNR

+NP <- res3$NP

+@ 

+

+<<snpAndCnSummaries, eval=FALSE>>=

+A <- res$A

+B <- res$B

+calls <- res2$calls

+conf <- res2$conf

+SNR <- res2$SNR

+NP <- res3$NP

+rm(res, res2, res3)

+gc()

+@ 

+

+Make a histogram of the signal to noise ratio for these samples:

+

+<<plotSnr, fig=TRUE>>=

+hist(SNR, xlab="SNR", main="")

+@ 

+

+We suggest excluding samples with a signal to noise ratio less than 5.

+We then extract the probes (rows) that are on chromosome 15 and the

+samples (columns) that have a suitable signal to noise ratio.

+

+<<subsetCrlmmOutput, eval=FALSE>>=

+A <- A[index[[1]], SNR > 5]

+B <- B[index[[1]], SNR > 5]

+calls <- calls[index[[1]], SNR > 5]

+conf <- conf[index[[1]], SNR > 5]

+NP <- NP[index[[2]], SNR > 5]

+rm(res, res2, res3); gc()

+@ 

+

+We want to adjust for batch effects.  In our experience, the chemistry

+plate is a good surrogate for batch, although this should be assessed

+on a study by study basis.  For samples processed by Broad, the plate

+is indicated by a capitalized five-letter word and is part of the

+sample name.  It is important that the plate (batch) variable is

+ordered the same as filenames in the preprocessed data.  For instance,

+to indicate plate in this example one can use the command

+

+<<specifyBatch>>=

+sns <- colnames(calls)

+sns[1]

+plates <- sapply(sns, function(x) strsplit(x, "_")[[1]][2])

+table(plates)

+@ 

+

+We are now ready to estimate the copy number for each batch.  In the

+current version of this package, one specifies an environment to which

+intermediate \R{} objects for copy number estimation are

+stored. Allele-specific estimates of copy number are also stored in

+this environment.

+

+<<copyNumberByBatch, eval=FALSE>>=

+e1 <- new.env()

+computeCnBatch(A=A,

+	       B=B,

+	       calls=calls,

+	       conf=conf,

+	       NP=NP,

+	       plate=plates,

+	       envir=e1, chrom=CHR, DF.PRIOR=75)

+@ 

+

+The DF.PRIOR indicates how much we will shrink SNP-specific estimates

+of the variance and correlation.

+

+<<accessingEstimates>>=

+copyA <- get("CA", e1)

+copyB <- get("CB", e1)

+copyT <- (copyA+copyB)/100

+copyT[copyT < 0] <- 0

+copyT[copyT > 6] <- 6

+@ 

+

+\section{Suggested plots}

+

+\paragraph{One sample at a time}

+

+<<annotation>>=

+data(snpProbes, package="genomewidesnp6Crlmm")

+data(cnProbes, package="genomewidesnp6Crlmm")

+position <- snpProbes[match(rownames(calls), rownames(snpProbes)), "position"]

+@ 

+

+To plot physical position versus copy number for the first sample:

+

+<<oneSample, fig=TRUE, width=8, height=4, include=FALSE>>=

+require(SNPchip)

+par(las=1)

+plotCytoband(as.character(CHR), ylim=c(0,7), cytoband.ycoords=c(6.5,7), label.cytoband=FALSE, main="Chr 15")

+points(position, copyT[, 1], pch=".", cex=2, , xaxt="n", col="grey70")

+axis(1, at=pretty(range(position)), labels=pretty(range(position))/1e6)

+axis(2, at=0:5, labels=0:5)

+mtext("position (Mb)", 1, line=2)

+mtext(expression(C[A] + C[B]), 2, line=2, las=3)

+@ 

+

+\begin{figure}

+  \includegraphics[width=0.9\textwidth]{copynumber-oneSample}

+  \caption{Total copy number (y-axis) for chromosome 22 plotted

+    against physical position (x-axis) for one sample.  }

+\end{figure}

+

+\paragraph{One SNP at a time}

+

+<<intermediateFiles>>=

+tau2A <- get("tau2A", e1)

+tau2B <- get("tau2B", e1)

+sig2A <- get("sig2A", e1)

+sig2B <- get("sig2B", e1)

+nuA <- get("nuA", e1)

+phiA <- get("phiA", e1)

+nuB <- get("nuB", e1)

+phiB <- get("phiB", e1)

+corr <- get("corr", e1)

+corrA.BB <- get("corrA.BB", e1)

+corrB.AA <- get("corrA.BB", e1)

+A <- get("A", e1)

+B <- get("B", e1)

+@ 

+

+Here, we plot the prediction regions for total copy number 2 and 3 for

+the first plate. Black plotting symbols are estimates from the first

+plate; light grey are points from other plates. (You could also draw

+prediction regions for 0-4 copies, but it gets crowded).  Notice that

+there is little evidence of a plate effect for this SNP.

+

+<<predictionRegion, fig=TRUE, width=8, height=8, include=FALSE>>=

+par(las=1)

+p <- 1 ##indicates plate

+J <- grep(unique(plates)[p], sns) ##sample indices for this plate

+ylim <- c(6.5,13)

+I <- which(phiA > 10 & phiB > 10)

+i <- I[1]

+plot(log2(A[i, ]), log2(B[i, ]), pch=as.character(calls[i, ]), col="grey60", cex=0.9, ylim=ylim, xlim=ylim, xlab="A", ylab="B")

+points(log2(A[i, J]), log2(B[i, J]), col="black", pch=as.character(calls[i, J]))

+for(CT in 2:3){

+	if(CT == 2) ellipse.col <- "brown" else ellipse.col <- "purple"

+	for(CA in 0:CT){

+		CB <- CT-CA

+		A.scale <- sqrt(tau2A[i, p]*(CA==0) + sig2A[i, p]*(CA > 0))

+		B.scale <- sqrt(tau2B[i, p]*(CB==0) + sig2B[i, p]*(CB > 0))

+		scale <- c(A.scale, B.scale)

+		if(CA == 0 & CB > 0) rho <- corrA.BB[i, p]

+		if(CA > 0 & CB == 0) rho <- corrB.AA[i, p]

+		if(CA > 0 & CB > 0) rho <- corr[i, p]		

+		lines(ellipse(x=rho, centre=c(log2(nuA[i, p]+CA*phiA[i, p]), log2(nuB[i, p]+CB*phiB[i, p])),

+			      scale=scale), col=ellipse.col, lwd=2)

+	}

+}

+legend("topright", lwd=3, col=c("black", "purple"), legend=c("2 copies", "3 copies"), bty="n")

+@ 

+

+\begin{figure}

+  \includegraphics[width=0.9\textwidth]{copynumber-predictionRegion}

+  \caption{Prediction regions for copy number 2 and 3 for one SNP.

+    The plate effects are negligible for this SNP and we only plot the

+    prediction region for the first plate.}

+\end{figure}

+

+Here is one way to identify a SNP with a large plate effect -- look

+for shifts in the prediction regions (here I'm lucking for large

+shifts in the A direction). 

+

+<<shifts>>=

+shiftA <- nuA+2*phiA 

+maxA <- apply(shiftA, 1, max, na.rm=T)

+minA <- apply(shiftA, 1, min, na.rm=T)

+d <- maxA - minA

+hist(d)

+index <- which(d > 3000)[2]

+plate1 <- which(shiftA[index, ] == maxA[index])

+plate2 <- which(shiftA[index, ] == minA[index])

+@ 

+

+Now plot the predictions for the two plates.  All the ellipses are

+prediction regions for copy number 2.  Note that while I looked for

+shifts in the A direction, the shift often occurs in both the B and A

+dimensions.

+

+<<plateEffect, fig=TRUE, width=8, height=8, include=FALSE>>=

+par(las=1)

+i <- index

+J1 <- grep(unique(plates)[plate1], sns) ##sample indices for this plate

+J2 <- grep(unique(plates)[plate2], sns) ##sample indices for this plate

+ylim <- c(6.5,13)

+plot(log2(A[i, ]), log2(B[i, ]), pch=as.character(calls[i, ]), col="grey60", cex=0.9, ylim=ylim, xlim=ylim, xlab="A", ylab="B")

+points(log2(A[i, J1]), log2(B[i, J1]), col="brown", pch=as.character(calls[i, J1]))

+points(log2(A[i, J2]), log2(B[i, J2]), col="blue", pch=as.character(calls[i, J2]))

+p <- plate1

+CT <- 2

+ellipse.col <- "brown" 

+for(CA in 0:CT){

+	CB <- CT-CA

+	A.scale <- sqrt(tau2A[i, p]*(CA==0) + sig2A[i, p]*(CA > 0))

+	B.scale <- sqrt(tau2B[i, p]*(CB==0) + sig2B[i, p]*(CB > 0))

+	scale <- c(A.scale, B.scale)

+	lines(ellipse(x=corr[i, p]*(CA>0)*(CB>0), centre=c(log2(nuA[i, p]+CA*phiA[i, p]), log2(nuB[i, p]+CB*phiB[i, p])),

+		      scale=scale), col=ellipse.col, lwd=2)

+}

+ellipse.col <- "blue" 

+p <- plate2

+for(CA in 0:CT){

+	CB <- CT-CA

+	A.scale <- sqrt(tau2A[i, p]*(CA==0) + sig2A[i, p]*(CA > 0))

+	B.scale <- sqrt(tau2B[i, p]*(CB==0) + sig2B[i, p]*(CB > 0))

+	scale <- c(A.scale, B.scale)

+	lines(ellipse(x=corr[i, p]*(CA>0)*(CB>0), centre=c(log2(nuA[i, p]+CA*phiA[i, p]), log2(nuB[i, p]+CB*phiB[i, p])),

+		      scale=scale), col=ellipse.col, lwd=2)

+}

+legend("topright", lwd=3, col=c("brown", "blue"), legend=c("2 copies, plate 1", "2 copies, plate 2 "), bty="n")

+@ 

+

+\begin{figure}

+  \includegraphics[width=0.9\textwidth]{copynumber-plateEffect}

+  \caption{The prediction regions for copy number 2 shift when there

+    is a large plate effect.}

+\end{figure}

+

+Alternatively, loop through a large number of SNPs to see how these

+regions move

+<<codeForLoop, eval=FALSE>>=

+par(las=1, pty="s", ask=TRUE)

+for(i in I[1:50]){

+	plot(log2(A[i, ]), log2(B[i, ]), pch=as.character(calls[i, ]), col=col[DS], cex=0.9, ylim=ylim, xlim=ylim)

+	for(CT in 2:3){

+		if(CT == 2) ellipse.col <- "brown" else ellipse.col <- "purple"

+		for(CA in 0:CT){

+			CB <- CT-CA

+			A.scale <- sqrt(tau2A[i, p]*(CA==0) + sig2A[i, p]*(CA > 0))

+			B.scale <- sqrt(tau2B[i, p]*(CB==0) + sig2B[i, p]*(CB > 0))

+			scale <- c(A.scale, B.scale)

+			lines(ellipse(x=corr[i, p]*(CA>0)*(CB>0), centre=c(log2(nuA[i, p]+CA*phiA[i, p]), log2(nuB[i, p]+CB*phiB[i, p])),

+				      scale=scale), col=ellipse.col, lwd=2)

+		}