git-svn-id: file:///home/git/hedgehog.fhcrc.org/bioconductor/trunk/madman/Rpacks/crlmm@48936 bc3139a8-67e5-0310-9ffc-ced21a209358
... | ... |
@@ -57,6 +57,7 @@ importFrom(ff, ffdf, physical.ff, physical.ffdf) |
57 | 57 |
|
58 | 58 |
exportClasses(CNSetLM, ffdf, list) |
59 | 59 |
exportMethods(open, "[", show, lM, lines, nu, phi, corr, sigma2, tau2) |
60 |
+exportMethods(CA, CB) |
|
60 | 61 |
export(crlmm, |
61 | 62 |
crlmmCopynumber, |
62 | 63 |
crlmmIllumina, |
... | ... |
@@ -73,8 +74,9 @@ export(crlmm, |
73 | 74 |
crlmmCopynumber2, crlmmCopynumberLD) |
74 | 75 |
export(constructIlluminaCNSet) |
75 | 76 |
export(linesCNSetLM) |
76 |
-export(linearModel.noweights, computeCN, fit.lm1, fit.lm2, fit.lm3, fit.lm4, construct, |
|
77 |
+export(computeCN, fit.lm1, fit.lm2, fit.lm3, fit.lm4, construct, |
|
77 | 78 |
dqrlsWrapper, nuphiAllele) |
79 |
+export(computeCopynumber) |
|
78 | 80 |
|
79 | 81 |
|
80 | 82 |
|
... | ... |
@@ -406,32 +406,6 @@ fit.wls <- function(allele, Ystar, W, Ns, autosome=TRUE){ |
406 | 406 |
return(betahat) |
407 | 407 |
} |
408 | 408 |
|
409 |
-## linear regression without weights -- design matrix is same for all snps |
|
410 |
-linearModel.noweights <- function(allele, Ystar, W, Ns){ |
|
411 |
- complete <- rowSums(is.na(W)) == 0 |
|
412 |
- if(any(!is.finite(W))){## | any(!is.finite(V))){ |
|
413 |
- i <- which(rowSums(!is.finite(W)) > 0) |
|
414 |
- stop("Possible zeros in the within-genotype estimates of the spread (vA, vB). ") |
|
415 |
- } |
|
416 |
- NOHET <- mean(Ns[, 2], na.rm=TRUE) < 0.05 |
|
417 |
- if(missing(allele)) stop("must specify allele") |
|
418 |
- if(allele == "A") X <- cbind(1, 2:0) else X <- cbind(1, 0:2) |
|
419 |
- if(NOHET) X <- X[-2, ] ##more than 1 X chromosome, but all homozygous |
|
420 |
- ##How to quickly generate Xstar, Xstar = diag(W) %*% X |
|
421 |
- Xstar <- apply(W, 1, generateX, X) |
|
422 |
- IXTX <- apply(Xstar, 2, generateIXTX, nrow=nrow(X)) |
|
423 |
- betahat <- matrix(NA, 2, nrow(Ystar)) |
|
424 |
- ses <- matrix(NA, 2, nrow(Ystar)) |
|
425 |
- for(i in 1:nrow(Ystar)){ |
|
426 |
- betahat[, i] <- crossprod(matrix(IXTX[, i], ncol(X), ncol(X)), crossprod(matrix(Xstar[, i], nrow=nrow(X)), Ystar[i, ])) |
|
427 |
- ssr <- sum((Ystar[i, ] - matrix(Xstar[, i], nrow(X), ncol(X)) %*% matrix(betahat[, i], ncol(X), 1))^2) |
|
428 |
- ses[, i] <- sqrt(diag(matrix(IXTX[, i], ncol(X), ncol(X)) * ssr)) |
|
429 |
- } |
|
430 |
- nu <- betahat[1, ] |
|
431 |
- phi <- betahat[2, ] |
|
432 |
- return(list(nu, phi)) |
|
433 |
-} |
|
434 |
- |
|
435 | 409 |
##nuphiAlleleX <- function(allele, Ystar, W, Ns, chrX=FALSE){ |
436 | 410 |
## complete <- rowSums(is.na(W)) == 0 |
437 | 411 |
## if(any(!is.finite(W))){## | any(!is.finite(V))){ |
... | ... |
@@ -1865,7 +1839,7 @@ nonpolymorphic <- function(object, cnOptions, tmp.objects){ |
1865 | 1839 |
nuA <- getParam(object, "nuA", batch) |
1866 | 1840 |
phiA <- getParam(object, "phiA", batch) |
1867 | 1841 |
} |
1868 |
- CA(object)[!isSnp(object), ] <- 1/phiA[!isSnp(object)]*(A - nuA[!isSnp(object)]) |
|
1842 |
+ ##CA(object)[!isSnp(object), ] <- 1/phiA[!isSnp(object)]*(A - nuA[!isSnp(object)]) |
|
1869 | 1843 |
sig2A <- getParam(object, "sig2A", batch) |
1870 | 1844 |
sig2A[!isSnp(object)] <- rowMAD(log2(A*normal), na.rm=TRUE)^2 |
1871 | 1845 |
object <- pr(object, "sig2A", batch, sig2A) |
... | ... |
@@ -93,7 +93,7 @@ setMethod("computeCopynumber", "CNSet", |
93 | 93 |
THR.NU.PHI=THR.NU.PHI, |
94 | 94 |
thresholdCopynumber=thresholdCopynumber) |
95 | 95 |
bias.adj <- cnOptions[["bias.adj"]] |
96 |
- if(bias.adj & all(is.na(nu(object, "A")[, 1])){ |
|
96 |
+ if(bias.adj & all(is.na(fData(object)$nuA_1))){ |
|
97 | 97 |
cnOptions[["bias.adj"]] <- FALSE |
98 | 98 |
} |
99 | 99 |
object <- computeCopynumber.CNSet(object, cnOptions) |
... | ... |
@@ -237,3 +237,11 @@ setMethod("corr", c("CNSetLM", "character"), function(object, allele){ |
237 | 237 |
return(res) |
238 | 238 |
}) |
239 | 239 |
|
240 |
+setMethod("CB", "CNSet", function(object){ |
|
241 |
+ ##assayDataElement(object, "CB") |
|
242 |
+ browser() |
|
243 |
+}) |
|
244 |
+setMethod("CA", "CNSet", function(object) { |
|
245 |
+ browser() |
|
246 |
+ ##assayDataElement(object, "CA") |
|
247 |
+}) |