Browse code

Updated sample.CNSetLM helpfile with coercion to v1.0.1 of the CNSet class definition

git-svn-id: file:///home/git/hedgehog.fhcrc.org/bioconductor/trunk/madman/Rpacks/crlmm@48950 bc3139a8-67e5-0310-9ffc-ced21a209358

Rob Scharp authored on 21/08/2010 02:49:06
Showing10 changed files

... ...
@@ -9,7 +9,7 @@ Description: Faster implementation of CRLMM specific to SNP 5.0 and 6.0 arrays,
9 9
 License: Artistic-2.0
10 10
 Depends: R (>= 2.11.0),
11 11
          methods,
12
-         Biobase (>= 2.7.2),
12
+         Biobase (>= 2.9.0),
13 13
          oligoClasses (>= 1.11.0)
14 14
 Imports: affyio (>= 1.15.2),
15 15
          ellipse,
... ...
@@ -54,10 +54,10 @@ importFrom(mvtnorm, dmvnorm)
54 54
 importFrom(ellipse, ellipse)
55 55
 
56 56
 importFrom(ff, ffdf, physical.ff, physical.ffdf)
57
+importClassesFrom(oligoClasses, ffdf)
57 58
 
58
-exportClasses(ffdf, list)
59 59
 exportMethods(lines)
60
-exportMethods(CA, CB, totalCopyNumber, initialize)
60
+exportMethods(CA, CB, totalCopyNumber)
61 61
 export(crlmm, 
62 62
        crlmmIllumina, 
63 63
        crlmmIllumina2,
... ...
@@ -68,7 +68,6 @@ export(crlmm,
68 68
        snprma2,
69 69
        crlmm2,
70 70
        genotype2, genotypeLD,
71
-       batch,
72 71
        crlmmCopynumber2, crlmmCopynumberLD)
73 72
 export(constructIlluminaCNSet)
74 73
 ##export(linesCNSet)
... ...
@@ -77,7 +76,7 @@ export(computeCN, fit.lm1, fit.lm2, fit.lm3, fit.lm4, construct,
77 76
 export(computeCopynumber, ACN)
78 77
 
79 78
 ## For debugging
80
-exportPattern("^[^\\.]")
79
+## exportPattern("^[^\\.]")
81 80
 
82 81
 
83 82
 
... ...
@@ -3,6 +3,4 @@ setOldClass("ff_matrix")
3 3
 setOldClass("ffdf")
4 4
 setClassUnion("ff_or_matrix", c("ff_matrix", "matrix", "ffdf"))
5 5
 setClassUnion("integerOrMissing", c("integer", "missing", "numeric"))
6
-##setClass("CNSetLM", contains="CNSet", representation(lM="list_or_ffdf"))
7
-##setClass("CNSetTest", contains="SnpSuperSet")
8 6
 
... ...
@@ -16,14 +16,8 @@ setGeneric("totalCopyNumber", function(object, i, j, ...) standardGeneric("total
16 16
 
17 17
 ## The generics below are for internal use with copy number methods
18 18
 ## If we keep them in oligoClasses, we need to export and document
19
-setGeneric("corr", function(object, allele) standardGeneric("corr"))
20
-setGeneric("nu", function(object, allele) standardGeneric("nu"))
21
-setGeneric("phi", function(object, allele) standardGeneric("phi"))
22
-setGeneric("sigma2", function(object, allele) standardGeneric("sigma2"))
23
-setGeneric("tau2", function(object, allele) standardGeneric("tau2"))
24
-
25 19
 setGeneric("nuA", function(object) standardGeneric("nuA"))
26
-setGeneric("nuB", function(object) standardGeneric("nuA"))
20
+setGeneric("nuB", function(object) standardGeneric("nuB"))
27 21
 setGeneric("phiA", function(object) standardGeneric("phiA"))
28 22
 setGeneric("phiB", function(object) standardGeneric("phiB"))
29 23
 setGeneric("sigma2A", function(object) standardGeneric("sigma2A"))
... ...
@@ -33,7 +27,6 @@ setGeneric("tau2B", function(object) standardGeneric("tau2B"))
33 27
 setGeneric("corrAA", function(object) standardGeneric("corrAA"))
34 28
 setGeneric("corrBB", function(object) standardGeneric("corrBB"))
35 29
 setGeneric("corrAB", function(object) standardGeneric("corrAB"))
36
-
37 30
 setGeneric("nuA<-", function(object, value) standardGeneric("nuA<-"))
38 31
 setGeneric("nuB<-", function(object, value) standardGeneric("nuB<-"))
39 32
 setGeneric("phiA<-", function(object, value) standardGeneric("phiA<-"))
... ...
@@ -29,7 +29,7 @@ getFeatureData.Affy <- function(cdfName, copynumber=FALSE){
29 29
 		load(file.path(path, "cnProbes.rda"))
30 30
 		cnProbes <- get("cnProbes")
31 31
 		snpIndex <- seq(along=gns)
32
-		npIndex <- seq(along=rownames(cnProbes)) + max(snpIndex) 
32
+		npIndex <- seq(along=rownames(cnProbes)) + max(snpIndex)
33 33
 		featurenames <- c(gns, rownames(cnProbes))
34 34
 	} else featurenames <- gns
35 35
 	fvarlabels=c("chromosome", "position", "isSnp")
... ...
@@ -76,28 +76,24 @@ construct <- function(filenames, cdfName, copynumber=FALSE,
76 76
 		featureData <- featureData[index, ]
77 77
 	}
78 78
 	nr <- nrow(featureData); nc <- length(sns)
79
-	ffObjects <- list(alleleA=initializeBigMatrix(name="A", nr, nc),
79
+	adObjects <- list(alleleA=initializeBigMatrix(name="A", nr, nc),
80 80
 			  alleleB=initializeBigMatrix(name="B", nr, nc),
81 81
 			  call=initializeBigMatrix(name="call", nr, nc),
82 82
 			  callProbability=initializeBigMatrix(name="callPr", nr,nc))
83
-##			  CA=initializeBigMatrix(name="CA", nr, nc),
84
-##			  CB=initializeBigMatrix(name="CB", nr, nc))
85 83
 	pd <- data.frame(matrix(NA, nc, 3), row.names=sns)
86 84
 	colnames(pd)=c("SKW", "SNR", "gender")
87 85
 	phenoData <- new("AnnotatedDataFrame", data=pd)
88
-	ffObjects <- lapply(ffObjects, function(x,sns) {colnames(x) <- sns; return(x)}, sns=sns)
89
-	callSet <- new("CNSet", 
90
-		       alleleA=ffObjects[["alleleA"]],
91
-		       alleleB=ffObjects[["alleleB"]],
92
-		       call=ffObjects[["call"]],
93
-		       callProbability=ffObjects[["callProbability"]],
94
-##		       CA=ffObjects[["CA"]],
95
-##		       CB=ffObjects[["CB"]],
86
+	adObjects <- lapply(adObjects, function(x,sns) {colnames(x) <- sns; return(x)}, sns=sns)
87
+	callSet <- new("CNSet",
88
+		       alleleA=adObjects[["alleleA"]],
89
+		       alleleB=adObjects[["alleleB"]],
90
+		       call=adObjects[["call"]],
91
+		       callProbability=adObjects[["callProbability"]],
96 92
 		       protocolData=protocolData,
97 93
 		       phenoData=phenoData,
98 94
 		       featureData=featureData,
99 95
 		       annotation=cdfName)
100
-	lM(callSet) <- initializeParamObject(list(featureNames(callSet), unique(protocolData(callSet)$batch)))	
96
+	lM(callSet) <- initializeParamObject(list(featureNames(callSet), unique(protocolData(callSet)$batch)))
101 97
 	return(callSet)
102 98
 }
103 99
 
... ...
@@ -156,7 +152,7 @@ construct <- function(filenames, cdfName, copynumber=FALSE,
156 152
 ##			np.index <- which(isSnp(callSet) == 0)
157 153
 ##			cnrmaRes <- cnrma(filenames=filenames[j],
158 154
 ##					  cdfName=cdfName,
159
-##					  row.names=featureNames(callSet)[np.index],				  
155
+##					  row.names=featureNames(callSet)[np.index],
160 156
 ##					  sns=sns[j],
161 157
 ##					  seed=seed,
162 158
 ##					  verbose=verbose)
... ...
@@ -1608,19 +1604,19 @@ processCEL2 <- function(i, filenames, row.names, A, SKW, seed, cdfName, pkgname)
1608 1604
 ##	return(NP)
1609 1605
 ##}
1610 1606
 
1611
-getFlags <- function(object, PHI.THR){
1612
-	batch <- unique(object$batch)
1613
-	nuA <- getParam(object, "nuA", batch)
1614
-	nuB <- getParam(object, "nuB", batch)
1615
-	phiA <- getParam(object, "phiA", batch)
1616
-	phiB <- getParam(object, "phiB", batch)
1617
-	negativeNus <- nuA < 1 | nuB < 1
1618
-	negativePhis <- phiA < PHI.THR | phiB < PHI.THR
1619
-	negativeCoef <- negativeNus | negativePhis
1620
-	notfinitePhi <- !is.finite(phiA) | !is.finite(phiB)
1621
-	flags <- negativeCoef | notfinitePhi
1622
-	return(flags)
1623
-}
1607
+##getFlags <- function(object, PHI.THR){
1608
+##	batch <- unique(object$batch)
1609
+##	nuA <- getParam(object, "nuA", batch)
1610
+##	nuB <- getParam(object, "nuB", batch)
1611
+##	phiA <- getParam(object, "phiA", batch)
1612
+##	phiB <- getParam(object, "phiB", batch)
1613
+##	negativeNus <- nuA < 1 | nuB < 1
1614
+##	negativePhis <- phiA < PHI.THR | phiB < PHI.THR
1615
+##	negativeCoef <- negativeNus | negativePhis
1616
+##	notfinitePhi <- !is.finite(phiA) | !is.finite(phiB)
1617
+##	flags <- negativeCoef | notfinitePhi
1618
+##	return(flags)
1619
+##}
1624 1620
 
1625 1621
 
1626 1622
 ##instantiateObjects <- function(object, cnOptions){
... ...
@@ -2737,14 +2733,11 @@ constructIlluminaAssayData <- function(np, snp, object, storage.mode="environmen
2737 2733
 	gt <- rbind(gt, emptyMatrix, deparse.level=0)[order.index,]
2738 2734
 	pr <- stripnames(snpCallProbability(object))
2739 2735
 	pr <- rbind(pr, emptyMatrix, deparse.level=0)[order.index, ]
2740
-	emptyMatrix <- matrix(integer(), nrow(A), ncol(A))	
2741 2736
 	aD <- assayDataNew(storage.mode,
2742 2737
 			   alleleA=A,
2743 2738
 			   alleleB=B,
2744 2739
 			   call=gt,
2745
-			   callProbability=pr,
2746
-			   CA=emptyMatrix,
2747
-			   CB=emptyMatrix)
2740
+			   callProbability=pr)
2748 2741
 }
2749 2742
 constructIlluminaCNSet <- function(crlmmResult,
2750 2743
 				   path,
... ...
@@ -2758,14 +2751,20 @@ constructIlluminaCNSet <- function(crlmmResult,
2758 2751
 	new.order <- order(fD$chromosome, fD$position)
2759 2752
 	fD <- fD[new.order, ]
2760 2753
 	aD <- constructIlluminaAssayData(cnAB, res, crlmmResult, order.index=new.order)
2761
-	protocolData(crlmmResult)$batch <- vector("integer", ncol(crlmmResult))
2754
+	##protocolData(crlmmResult)$batch <- vector("integer", ncol(crlmmResult))
2755
+	batch <- vector("integer", ncol(crlmmResult))
2762 2756
 	container <- new("CNSet", 
2763
-			 assayData=aD,
2757
+			 call=aD[["call"]],
2758
+			 callProbability=aD[["callProbability"]],
2759
+			 alleleA=aD[["alleleA"]],
2760
+			 alleleB=aD[["alleleB"]],
2764 2761
 			 phenoData=phenoData(crlmmResult),
2765 2762
 			 protocolData=protocolData(crlmmResult),
2766 2763
 			 featureData=fD,
2764
+			 batch=batch,
2767 2765
 			 annotation="human370v1c")
2768
-	lM(container) <- initializeParamObject(list(featureNames(container), unique(protocolData(container)$batch)))
2766
+##	lM(container) <- initializeParamObject(list(featureNames(container), unique(protocolData(container)$batch)))
2767
+	lM(container) <- initializeLmFrom(container)
2769 2768
 	container
2770 2769
 }
2771 2770
 				   
... ...
@@ -17,11 +17,7 @@ linearParamElementReplace <- function(obj, elt, value) {
17 17
 }
18 18
 
19 19
 
20
-setMethod("nu", c("CNSet", "character"), function(object, allele) nu(lM(object), allele))
21
-setMethod("phi", c("CNSet", "character"), function(object, allele) phi(lM(object), allele))
22
-setMethod("sigma2", c("CNSet", "character"), function(object, allele) phi(lM(object), allele))
23
-setMethod("tau2", c("CNSet", "character"), function(object, allele) phi(lM(object), allele))
24
-setMethod("corr", c("CNSet", "character"), function(object, allele) phi(lM(object), allele))
20
+
25 21
 
26 22
 setMethod("nuA", signature=signature(object="CNSet"), function(object) nu(object, "A"))
27 23
 setMethod("nuB", signature=signature(object="CNSet"), function(object) nu(object, "B"))
... ...
@@ -1,62 +1 @@
1
-setMethod("nu", c("LinearModelParameter", "character"), 
2
-	  function(object, allele){
3
-		  getValue <- function(allele){
4
-			  switch(allele,
5
-				 A="nuA",
6
-				 B="nuB",
7
-				 stop("allele must be 'A' or 'B'"))
8
-		  }	
9
-		  val <- getValue(allele)
10
-		  assayDataElement(object, val)
11
-})
12 1
 
13
-
14
-
15
-setMethod("phi", c("LinearModelParameter", "character"),
16
-	  function(object, allele){
17
-		  getValue <- function(allele){
18
-			  switch(allele,
19
-				 A="phiA",
20
-				 B="phiB",
21
-				 stop("allele must be 'A' or 'B'"))
22
-		  }	
23
-		  val <- getValue(allele)
24
-		  assayDataElement(object, val)
25
-	  })
26
-
27
-setMethod("sigma2", c("LinearModelParameter", "character"),
28
-	  function(object, allele){
29
-		  getValue <- function(allele){
30
-			  switch(allele,
31
-				 A="sigma2A",
32
-				 B="sigma2B",
33
-				 stop("allele must be 'A' or 'B'"))
34
-		  }
35
-		  val <- getValue(allele)
36
-		  assayDataElement(object, val)
37
-	  })
38
-
39
-setMethod("tau2", c("LinearModelParameter", "character"),
40
-	  function(object, allele){
41
-		  getValue <- function(allele){
42
-			  switch(allele,
43
-				 A="tau2A",
44
-				 B="tau2B",
45
-				 stop("allele must be 'A' or 'B'"))
46
-		  }
47
-		  val <- getValue(allele)		  
48
-		  assayDataElement(object, val)
49
-	  })
50
-
51
-setMethod("corr", c("LinearModelParameter", "character"),
52
-	  function(object, allele){
53
-		  getValue <- function(allele){
54
-			  switch(allele,
55
-				 AA="corrAA",
56
-				 AB="corrAB",
57
-				 BB="corrBB",
58
-				 stop("allele must be 'AA', 'AB', or 'BB'"))
59
-		  }
60
-		  val <- getValue(allele)		  
61
-		  assayDataElement(object, val)
62
-	  })
... ...
@@ -34,24 +34,24 @@ setMethod("snpNames", "eSet", function(object){
34 34
 setMethod("chromosome", "eSet", function(object) fData(object)$chromosome)
35 35
 setMethod("position", "eSet", function(object) fData(object)$position)
36 36
 
37
-setMethod("getParam", signature(object="eSet",
38
-				name="character",
39
-				batch="ANY"),
40
-	  function(object, name, batch){
41
-		  if(length(batch) > 1){
42
-			  warning("batch argument to getParam should have length 1.  Only using the first")
43
-			  batch <- batch[1]
44
-		  }
45
-		  getParam.SnpSuperSet(object, name, batch)
46
-})
37
+##setMethod("getParam", signature(object="eSet",
38
+##				name="character",
39
+##				batch="ANY"),
40
+##	  function(object, name, batch){
41
+##		  if(length(batch) > 1){
42
+##			  warning("batch argument to getParam should have length 1.  Only using the first")
43
+##			  batch <- batch[1]
44
+##		  }
45
+##		  getParam.SnpSuperSet(object, name, batch)
46
+##})
47 47
 
48
-setMethod("batch", "eSet", function(object){
49
-	if("batch" %in% varLabels(object)){
50
-		return(object$batch)
51
-	} else {
52
-		return(protocolData(object)$batch)
53
-	}
54
-})
48
+##setMethod("batch", "eSet", function(object){
49
+##	if("batch" %in% varLabels(object)){
50
+##		return(object$batch)
51
+##	} else {
52
+##		return(protocolData(object)$batch)
53
+##	}
54
+##})
55 55
 
56 56
 ##setMethod("combine", signature=signature(x="eSet", y="eSet"),
57 57
 ##	  function(x, y, ...){
... ...
@@ -81,21 +81,21 @@ setMethod("batch", "eSet", function(object){
81 81
 
82 82
 
83 83
 
84
-setMethod("pr", signature(object="eSet",
85
-			  name="character",
86
-			  batch="ANY",
87
-			  value="numeric"), 
88
-	  function(object, name, batch, value){
89
-		  label <- paste(name, batch, sep="_")
90
-		  colindex <- match(label, fvarLabels(object))
91
-		  if(length(colindex) == 1){
92
-			  fData(object)[, colindex] <- value
93
-		  } 
94
-		  if(is.na(colindex)){
95
-			  stop(paste(label, " not found in object"))
96
-		  }
97
-		  if(length(colindex) > 1){
98
-			  stop(paste(label, " not unique"))
99
-		  }
100
-		  object
101
-	  })
84
+##setMethod("pr", signature(object="eSet",
85
+##			  name="character",
86
+##			  batch="ANY",
87
+##			  value="numeric"), 
88
+##	  function(object, name, batch, value){
89
+##		  label <- paste(name, batch, sep="_")
90
+##		  colindex <- match(label, fvarLabels(object))
91
+##		  if(length(colindex) == 1){
92
+##			  fData(object)[, colindex] <- value
93
+##		  } 
94
+##		  if(is.na(colindex)){
95
+##			  stop(paste(label, " not found in object"))
96
+##		  }
97
+##		  if(length(colindex) > 1){
98
+##			  stop(paste(label, " not unique"))
99
+##		  }
100
+##		  object
101
+##	  })
102 102
deleted file mode 100644
... ...
@@ -1,109 +0,0 @@
1
-\name{linearModelParams}
2
-\alias{lM}
3
-\alias{corr}
4
-\alias{nu}
5
-\alias{phi}
6
-\alias{sigma2}
7
-\alias{tau2}
8
-\title{
9
-	Accessors for linear model parameters.
10
-}
11
-\description{
12
-	
13
-	Accessors for linear model parameters.
14
-
15
-}
16
-
17
-\usage{
18
-lM(object)
19
-corr(object, allele)
20
-nu(object, allele)
21
-phi(object, allele)
22
-sigma2(object, allele)
23
-tau2(object, allele)
24
-}
25
-
26
-\arguments{
27
-  \item{object}{  An object of class \code{CNSetLM}}
28
-  \item{allele}{
29
-   Character string. Valid entries for most accessors are 'A' or 'B'.
30
-   For the \code{corr} accessor, valid entries are 'AA', 'AB', or 'BB'.
31
-   }
32
-}
33
-
34
-\details{
35
-
36
-\code{lM}:  Extracts entire list of linear model parameters.
37
-
38
-\code{corr}: The within-genotype correlation of log2(A) and log2(B) intensities.
39
-
40
-\code{nu}: The intercept for the linear model.  The linear model is
41
-fit to the A and B alleles independently.
42
-
43
-\code{phi}: The slope for the linear model.  The linear model is fit
44
-independently to the A and B alleles.
45
-
46
-\code{sigma2}: For allele A, sigma2 is calculated as the squared MAD
47
-of the log2(intensity) for allele 'A' among subjects with genotype AA.
48
-For allele B, sigma2 is calculated as the squared MAD of the
49
-log2(intensity) for allele 'B' among subjects with genotype BB.
50
-sigma2 can be interpreted as a robust estimate of the signal variance.
51
-
52
-\code{tau2}: For allele A, tau2 is calculated as the squared MAD of
53
-the log2(intensity) for allele 'A' among subjects with genotype BB.
54
-For allele B, tau2 is calculated as the squared MAD of the
55
-log2(intensity) for allele 'B' among subjects with genotype AA.  tau2
56
-can be interpeted as a robust estimate of the background variance.
57
-
58
-}
59
-
60
-\value{
61
-
62
-	A matrix or \code{ff} object.
63
-
64
-}
65
-
66
-\author{
67
-R. Scharpf
68
-}
69
-
70
-\seealso{
71
-	\code{\link{CNSet-class}}
72
-}
73
-\examples{
74
-## object with ff class
75
-if(require("ff")){
76
-	data(sample.CNSetLMff)
77
-	invisible(open(sample.CNSetLMff))
78
-	class(lM(sample.CNSetLMff))
79
-	params <- lM(sample.CNSetLMff)
80
-	nuA <- nu(sample.CNSetLMff, "A")
81
-	nuB <- nu(sample.CNSetLMff, "B")
82
-	phA <- phi(sample.CNSetLMff, "A")
83
-	phB <- phi(sample.CNSetLMff, "B")
84
-	sig2A <- sigma2(sample.CNSetLMff, "A") 
85
-	sig2B <- sigma2(sample.CNSetLMff, "B")
86
-	tau2A <- tau2(sample.CNSetLMff, "A")
87
-	tau2B <- tau2(sample.CNSetLMff, "B")
88
-	corrAA <- corr(sample.CNSetLMff, "AA")
89
-	corrBB <- corr(sample.CNSetLMff, "BB")
90
-	corrAB <- corr(sample.CNSetLMff, "AB")
91
-	invisible(close(sample.CNSetLMff))
92
-}
93
-## object with matrix class
94
-data(sample.CNSetLM)
95
-class(lM(sample.CNSetLM))
96
-nuA <- nu(sample.CNSetLM, "A")
97
-nuB <- nu(sample.CNSetLM, "B")
98
-phA <- phi(sample.CNSetLM, "A")
99
-phB <- phi(sample.CNSetLM, "B")
100
-sig2A <- sigma2(sample.CNSetLM, "A") 
101
-sig2B <- sigma2(sample.CNSetLM, "B")
102
-tau2A <- tau2(sample.CNSetLM, "A")
103
-tau2B <- tau2(sample.CNSetLM, "B")
104
-corrAA <- corr(sample.CNSetLM, "AA")
105
-corrBB <- corr(sample.CNSetLM, "BB")
106
-corrAB <- corr(sample.CNSetLM, "AB")
107
-}
108
-\keyword{manip}
109
-
... ...
@@ -13,6 +13,9 @@
13 13
 \usage{data(sample.CNSetLM)}
14 14
 \format{
15 15
 
16
+	This class has been deprecated.  See example below for how to
17
+	update an existing 'CNSetLM' object to class 'CNSet'.
18
+ 
16 19
 	The data illustrates the \code{CNSetLM-class}, with
17 20
 	\code{assayData} containing the quantile-normalized
18 21
 	intensities for the A and B alleles, genotype calls and
... ...
@@ -22,6 +25,10 @@
22 25
 	
23 26
 }
24 27
 \examples{
28
+## class CNSetLM has been deprecated
25 29
 data(sample.CNSetLM)
30
+## update to class CNSet
31
+cnSet <- as(sample.CNSetLM, "CNSet")
32
+all(isCurrent(cnSet))
26 33
 }
27 34
 \keyword{datasets}