Browse code

added plot-methods.R. Added back a versionNumber check (from r45050) in crlmm-illumina.R

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

Rob Scharp authored on 25/05/2010 17:47:57
Showing3 changed files

... ...
@@ -1,7 +1,7 @@
1 1
 Package: crlmm
2 2
 Type: Package
3 3
 Title: Genotype Calling (CRLMM) and Copy Number Analysis tool for Affymetrix SNP 5.0 and 6.0 and Illumina arrays.
4
-Version: 1.7.2
4
+Version: 1.7.3
5 5
 Date: 2010-05-25
6 6
 Author: Rafael A Irizarry, Benilton S Carvalho <carvalho@bclab.org>, Robert Scharpf <rscharpf@jhsph.edu>, Matt Ritchie <mritchie@wehi.edu.au>
7 7
 Maintainer: Benilton S Carvalho <carvalho@bclab.org>, Robert Scharpf <rscharpf@jhsph.edu>, Matt Ritchie <mritchie@wehi.EDU.AU>
... ...
@@ -369,6 +369,9 @@ readIDAT <- function(idatFile){
369 369
 
370 370
   versionNumber <- readBin(tempCon, "integer", n=1, size=8, 
371 371
                            endian="little", signed=FALSE)
372
+
373
+  if(versionNumber<3)
374
+	  stop("Older style IDAT files not supported:  consider updating your scanner settings")
372 375
   
373 376
   nFields <- readBin(tempCon, "integer", n=1, size=4, 
374 377
                      endian="little", signed=FALSE)
375 378
new file mode 100644
... ...
@@ -0,0 +1,79 @@
1
+setOldClass("ellipse")
2
+setMethod("plot", c("numeric", "CNSetLM"), function(x, y, copynumber=2, batch, plot.it=TRUE, ...){
3
+	plotCNSetLM(x, y, copynumber, batch, plot.it, ...)
4
+	})
5
+plotCNSetLM <- function(x, y, copynumber, batch, plot.it, ...){
6
+	require(ellipse)
7
+	object <- y
8
+	I <- x
9
+	if(missing(batch)) batch <- seq(along=unique(protocolData(object)$batch))
10
+	logA <- log2(as.matrix(A(object)[I, ]))
11
+	logB <- log2(as.matrix(B(object)[I, ]))
12
+	if(!"ylim" %in% names(list(...))){
13
+		ylim <- range(c(as.numeric(logA), as.numeric(logB)), na.rm=TRUE)
14
+	}
15
+	hasxlab <- function(...) !all(is.na(pmatch(names(list(...)), "lab")))
16
+	hascol <- function(...) !all(is.na(pmatch(names(list(...)), "col")))
17
+	hasbg <- function(...) !all(is.na(pmatch(names(list(...)), "bg")))	
18
+##	if(color.by.genotype){
19
+##		gt <- as.matrix(calls(object)[I, ])
20
+##		if(!hascol(...)) col <- gt
21
+##	}
22
+	ffIsLoaded <- class(calls(object))[[1]] == "ff"
23
+	if(ffIsLoaded){
24
+		nuA <- (physical(lM(object))[["nuA"]])[I, , drop=FALSE]
25
+		nuB <- (physical(lM(object))[["nuB"]])[I, , drop=FALSE]
26
+		phiA <- (physical(lM(object))[["phiA"]])[I, ,drop=FALSE]
27
+		phiB <- (physical(lM(object))[["phiB"]])[I, ,drop=FALSE]
28
+		tau2A <- (physical(lM(object))[["tau2A"]])[I, ,drop=FALSE]
29
+		tau2B <- (physical(lM(object))[["tau2B"]])[I, ,drop=FALSE]
30
+		sigma2A <- (physical(lM(object))[["sig2A"]])[I, ,drop=FALSE]
31
+		sigma2B <- (physical(lM(object))[["sig2B"]])[I, ,drop=FALSE]
32
+		corrAB <- (physical(lM(object))[["corrAB"]])[I, ,drop=FALSE]
33
+		corrAA <- (physical(lM(object))[["corrAA"]])[I, ,drop=FALSE]
34
+		corrBB <- (physical(lM(object))[["corrBB"]])[I, ,drop=FALSE]
35
+	} else {
36
+		nuA <- lM(object)[["nuA"]][I, , drop=FALSE]
37
+		nuB <- lM(object)[["nuB"]][I, , drop=FALSE]
38
+		phiA <- lM(object)[["phiA"]][I, ,drop=FALSE]
39
+		phiB <- lM(object)[["phiB"]][I, ,drop=FALSE]
40
+		tau2A <- lM(object)[["tau2A"]][I, ,drop=FALSE]
41
+		tau2B <- lM(object)[["tau2B"]][I, ,drop=FALSE]
42
+		sigma2A <- lM(object)[["sig2A"]][I, ,drop=FALSE]
43
+		sigma2B <- lM(object)[["sig2B"]][I, ,drop=FALSE]
44
+		corrAB <- lM(object)[["corrAB"]][I, ,drop=FALSE]
45
+		corrAA <- lM(object)[["corrAA"]][I, ,drop=FALSE]
46
+		corrBB <- lM(object)[["corrBB"]][I, ,drop=FALSE]
47
+	}
48
+	for(i in seq(along=I)){
49
+		for(b in batch){
50
+			jj <- which(batch(object) == batch)
51
+			if(plot.it){
52
+				if(hasxlab(...)){
53
+					plot(logA[i, jj], logB[i, jj], ...)
54
+				} else{
55
+					plot(logA[i, jj], logB[i, jj], xlab="log(A)", ylab="log(B)", ...)
56
+				}
57
+			}
58
+			if(all(is.na(nuA[, b]))) {
59
+				message("Parameter estimates for batch ", b, " not available")
60
+				next()
61
+			}
62
+			for(CN in copynumber){
63
+				for(CA in 0:CN){
64
+					CB <- CN-CA
65
+					A.scale <- sqrt(tau2A[i, b]*(CA==0) + sigma2A[i, b]*(CA > 0))
66
+					B.scale <- sqrt(tau2B[i, b]*(CB==0) + sigma2B[i, b]*(CB > 0))
67
+					scale <- c(A.scale, B.scale)
68
+					if(CA == 0 & CB > 0) rho <- corrBB[i, b]
69
+					if(CA > 0 & CB == 0) rho <- corrAA[i, b]
70
+					if(CA > 0 & CB > 0) rho <- corrAB[i, b]
71
+					if(CA == 0 & CB == 0) rho <- 0
72
+					lines(ellipse(x=rho, centre=c(log2(nuA[i, b]+CA*phiA[i, b]),
73
+							     log2(nuB[i, b]+CB*phiB[i,b])),
74
+						      scale=scale), ...)
75
+				}
76
+			}	
77
+		}
78
+	}
79
+}