git-svn-id: file:///home/git/hedgehog.fhcrc.org/bioconductor/trunk/madman/Rpacks/crlmm@48617 bc3139a8-67e5-0310-9ffc-ced21a209358
... | ... |
@@ -35,7 +35,7 @@ importFrom(oligoClasses, chromosome2integer, celfileDate, list.celfiles, |
35 | 35 |
copyNumber, initializeBigMatrix, initializeBigVector) |
36 | 36 |
|
37 | 37 |
importFrom(graphics, abline, axis, layout, legend, mtext, par, plot, |
38 |
- polygon, rect, segments, text, points, boxplot) |
|
38 |
+ polygon, rect, segments, text, points, boxplot, lines) |
|
39 | 39 |
|
40 | 40 |
importFrom(grDevices, grey) |
41 | 41 |
|
... | ... |
@@ -56,7 +56,7 @@ importFrom(ellipse, ellipse) |
56 | 56 |
importFrom(ff, ffdf) |
57 | 57 |
|
58 | 58 |
exportClasses(CNSetLM, ffdf, list) |
59 |
-exportMethods(open, "[", show, lM) |
|
59 |
+exportMethods(open, "[", show, lM, lines) |
|
60 | 60 |
export(crlmm, |
61 | 61 |
crlmmCopynumber, |
62 | 62 |
crlmmIllumina, |
... | ... |
@@ -72,7 +72,7 @@ export(crlmm, |
72 | 72 |
crlmmCopynumber2, crlmmCopynumberLD) |
73 | 73 |
##export(initializeParamObject, biasAdjNP2) |
74 | 74 |
export(construct) |
75 |
-export(plotCNSetLM) |
|
75 |
+ |
|
76 | 76 |
|
77 | 77 |
|
78 | 78 |
|
... | ... |
@@ -1,84 +1,64 @@ |
1 | 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, ...){ |
|
2 |
+setMethod("lines", c("CNSetLM"), function(x, y, batch, copynumber, ...){ |
|
3 |
+ linesCNSetLM(x, y, batch, copynumber, ...) |
|
4 |
+}) |
|
5 |
+linesCNSetLM <- function(x, y, batch, copynumber, x.axis="A", ...){ |
|
6 | 6 |
require(ellipse) |
7 |
- object <- y |
|
8 |
- stopifnot(length(batch) == 1) |
|
9 |
- stopifnot(batch %in% batch(object)) |
|
10 |
- sample.index <- which(batch(object) %in% batch) |
|
11 |
- I <- x |
|
12 |
- if(missing(batch)) batch <- seq(along=unique(protocolData(object)$batch)) |
|
13 |
- logA <- log2(as.matrix(A(object)[I, sample.index])) |
|
14 |
- dimnames(logA) <- NULL |
|
15 |
- logB <- log2(as.matrix(B(object)[I, sample.index])) |
|
16 |
- dimnames(logB) <- NULL |
|
17 |
- if(!"ylim" %in% names(list(...))){ |
|
18 |
- ylim <- range(c(as.numeric(logA), as.numeric(logB)), na.rm=TRUE) |
|
19 |
- } |
|
20 |
- hasxlab <- function(...) !all(is.na(pmatch(names(list(...)), "lab"))) |
|
21 |
- hascol <- function(...) !all(is.na(pmatch(names(list(...)), "col"))) |
|
22 |
- hasbg <- function(...) !all(is.na(pmatch(names(list(...)), "bg"))) |
|
23 |
-## if(color.by.genotype){ |
|
24 |
-## gt <- as.matrix(calls(object)[I, ]) |
|
25 |
-## if(!hascol(...)) col <- gt |
|
26 |
-## } |
|
7 |
+ object <- x |
|
8 |
+ I <- y |
|
9 |
+ stopifnot(length(I) == 1) |
|
27 | 10 |
calls.class <- class(calls(object)[[1]]) |
28 | 11 |
ffIsLoaded <- calls.class[1] == "ff_matrix" | calls.class[1] == "ffdf" | calls.class[1]=="ff" |
29 | 12 |
column <- grep(batch, unique(batch(object))) |
30 | 13 |
stopifnot(length(column) == 1) |
31 | 14 |
if(ffIsLoaded){ |
32 |
- nuA <- (physical(lM(object))[["nuA"]])[I, column , drop=FALSE] |
|
33 |
- nuB <- (physical(lM(object))[["nuB"]])[I, column , drop=FALSE] |
|
34 |
- phiA <- (physical(lM(object))[["phiA"]])[I, column ,drop=FALSE] |
|
35 |
- phiB <- (physical(lM(object))[["phiB"]])[I, column ,drop=FALSE] |
|
36 |
- tau2A <- (physical(lM(object))[["tau2A"]])[I, column ,drop=FALSE] |
|
37 |
- tau2B <- (physical(lM(object))[["tau2B"]])[I, column ,drop=FALSE] |
|
38 |
- sigma2A <- (physical(lM(object))[["sig2A"]])[I, column ,drop=FALSE] |
|
39 |
- sigma2B <- (physical(lM(object))[["sig2B"]])[I, column ,drop=FALSE] |
|
40 |
- corrAB <- (physical(lM(object))[["corrAB"]])[I, column ,drop=FALSE] |
|
41 |
- corrAA <- (physical(lM(object))[["corrAA"]])[I, column ,drop=FALSE] |
|
42 |
- corrBB <- (physical(lM(object))[["corrBB"]])[I, column ,drop=FALSE] |
|
15 |
+ nuA <- (physical(lM(object))[["nuA"]])[I, column , drop=TRUE] |
|
16 |
+ nuB <- (physical(lM(object))[["nuB"]])[I, column , drop=TRUE] |
|
17 |
+ phiA <- (physical(lM(object))[["phiA"]])[I, column ,drop=TRUE] |
|
18 |
+ phiB <- (physical(lM(object))[["phiB"]])[I, column ,drop=TRUE] |
|
19 |
+ tau2A <- (physical(lM(object))[["tau2A"]])[I, column ,drop=TRUE] |
|
20 |
+ tau2B <- (physical(lM(object))[["tau2B"]])[I, column ,drop=TRUE] |
|
21 |
+ sigma2A <- (physical(lM(object))[["sig2A"]])[I, column ,drop=TRUE] |
|
22 |
+ sigma2B <- (physical(lM(object))[["sig2B"]])[I, column ,drop=TRUE] |
|
23 |
+ corrAB <- (physical(lM(object))[["corrAB"]])[I, column ,drop=TRUE] |
|
24 |
+ corrAA <- (physical(lM(object))[["corrAA"]])[I, column ,drop=TRUE] |
|
25 |
+ corrBB <- (physical(lM(object))[["corrBB"]])[I, column ,drop=TRUE] |
|
43 | 26 |
} else { |
44 |
- nuA <- lM(object)[["nuA"]][I, column , drop=FALSE] |
|
45 |
- nuB <- lM(object)[["nuB"]][I, column , drop=FALSE] |
|
46 |
- phiA <- lM(object)[["phiA"]][I, column ,drop=FALSE] |
|
47 |
- phiB <- lM(object)[["phiB"]][I, column ,drop=FALSE] |
|
48 |
- tau2A <- lM(object)[["tau2A"]][I, column ,drop=FALSE] |
|
49 |
- tau2B <- lM(object)[["tau2B"]][I, column ,drop=FALSE] |
|
50 |
- sigma2A <- lM(object)[["sig2A"]][I, column ,drop=FALSE] |
|
51 |
- sigma2B <- lM(object)[["sig2B"]][I, column ,drop=FALSE] |
|
52 |
- corrAB <- lM(object)[["corrAB"]][I, column ,drop=FALSE] |
|
53 |
- corrAA <- lM(object)[["corrAA"]][I, column ,drop=FALSE] |
|
54 |
- corrBB <- lM(object)[["corrBB"]][I, column ,drop=FALSE] |
|
27 |
+ nuA <- lM(object)[["nuA"]][I, column , drop=TRUE] |
|
28 |
+ nuB <- lM(object)[["nuB"]][I, column , drop=TRUE] |
|
29 |
+ phiA <- lM(object)[["phiA"]][I, column ,drop=TRUE] |
|
30 |
+ phiB <- lM(object)[["phiB"]][I, column ,drop=TRUE] |
|
31 |
+ tau2A <- lM(object)[["tau2A"]][I, column ,drop=TRUE] |
|
32 |
+ tau2B <- lM(object)[["tau2B"]][I, column ,drop=TRUE] |
|
33 |
+ sigma2A <- lM(object)[["sig2A"]][I, column ,drop=TRUE] |
|
34 |
+ sigma2B <- lM(object)[["sig2B"]][I, column ,drop=TRUE] |
|
35 |
+ corrAB <- lM(object)[["corrAB"]][I, column ,drop=TRUE] |
|
36 |
+ corrAA <- lM(object)[["corrAA"]][I, column ,drop=TRUE] |
|
37 |
+ corrBB <- lM(object)[["corrBB"]][I, column ,drop=TRUE] |
|
55 | 38 |
} |
56 |
- for(i in seq(along=I)){ |
|
57 |
- if(plot.it){ |
|
58 |
- if(hasxlab(...)){ |
|
59 |
- plot(logA[i, ], logB[i, ], ...) |
|
60 |
- } else{ |
|
61 |
- plot(logA[i, ], logB[i, ], xlab="log(A)", ylab="log(B)", ...) |
|
62 |
- } |
|
63 |
- } |
|
64 |
- if(all(is.na(nuA))) { |
|
65 |
- message("Parameter estimates for batch ", b, " not available") |
|
66 |
- next() |
|
67 |
- } |
|
68 |
- for(CN in copynumber){ |
|
69 |
- for(CA in 0:CN){ |
|
70 |
- CB <- CN-CA |
|
71 |
- A.scale <- sqrt(tau2A[i, ]*(CA==0) + sigma2A[i, ]*(CA > 0)) |
|
72 |
- B.scale <- sqrt(tau2B[i, ]*(CB==0) + sigma2B[i, ]*(CB > 0)) |
|
73 |
- scale <- c(A.scale, B.scale) |
|
74 |
- if(CA == 0 & CB > 0) rho <- corrBB[i, ] |
|
75 |
- if(CA > 0 & CB == 0) rho <- corrAA[i, ] |
|
76 |
- if(CA > 0 & CB > 0) rho <- corrAB[i, ] |
|
77 |
- if(CA == 0 & CB == 0) rho <- 0 |
|
78 |
- lines(ellipse(x=rho, centre=c(log2(nuA[i, ]+CA*phiA[i, ]), |
|
79 |
- log2(nuB[i, ]+CB*phiB[i,])), |
|
39 |
+ if(all(is.na(nuA))) { |
|
40 |
+ message("Parameter estimates for batch ", b, " not available") |
|
41 |
+ next() |
|
42 |
+ } |
|
43 |
+ for(CN in copynumber){ |
|
44 |
+ for(CA in 0:CN){ |
|
45 |
+ CB <- CN-CA |
|
46 |
+ A.scale <- sqrt(tau2A*(CA==0) + sigma2A*(CA > 0)) |
|
47 |
+ B.scale <- sqrt(tau2B*(CB==0) + sigma2B*(CB > 0)) |
|
48 |
+ scale <- c(A.scale, B.scale) |
|
49 |
+ if(CA == 0 & CB > 0) rho <- corrBB |
|
50 |
+ if(CA > 0 & CB == 0) rho <- corrAA |
|
51 |
+ if(CA > 0 & CB > 0) rho <- corrAB |
|
52 |
+ if(CA == 0 & CB == 0) rho <- 0 |
|
53 |
+ if(x.axis=="A"){ |
|
54 |
+ lines(ellipse(x=rho, centre=c(log2(nuA+CA*phiA), |
|
55 |
+ log2(nuB+CB*phiB)), |
|
80 | 56 |
scale=scale), ...) |
57 |
+ } else { |
|
58 |
+ lines(ellipse(x=rho, centre=c(log2(nuB+CB*phiB), |
|
59 |
+ log2(nuA+CA*phiA)), |
|
60 |
+ scale=rev(scale)), ...) |
|
81 | 61 |
} |
82 |
- } |
|
83 |
- } |
|
62 |
+ } |
|
63 |
+ } |
|
84 | 64 |
} |