* collab:
bump version
made cnSetExample smaller. Fix notes
Trying to revert bad commit
remove cn-functions. update description
comment most of cn-functions.r
Resaved rdas
update data/cnSetExample.rda and data/cnSetExample2.rda
bump version
coercion method from CNSet to oligoSnpSet makes integer matrices of BAFs and lrr's
import ff_or_matrix from oligoClasses. bump dependency on oligoClasses version. Use library(oligoClasses) in some of the crlmm examples.
Cleaning pkg loading process: work still required
move Biobase and methods to imports
git-svn-id: file:///home/git/hedgehog.fhcrc.org/bioconductor/trunk/madman/Rpacks/crlmm@64324 bc3139a8-67e5-0310-9ffc-ced21a209358
... | ... |
@@ -1,21 +1,21 @@ |
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.13.12 |
|
5 |
-Date: 2010-12-10 |
|
6 |
-Author: Benilton S Carvalho <Benilton.Carvalho@cancer.org.uk>, Robert Scharpf <rscharpf@jhsph.edu>, Matt Ritchie <mritchie@wehi.edu.au>, Ingo Ruczinski <iruczins@jhsph.edu>, Rafael A Irizarry |
|
4 |
+Version: 1.13.14 |
|
5 |
+Author: Benilton S Carvalho, Robert Scharpf, Matt Ritchie, Ingo Ruczinski, Rafael A Irizarry |
|
7 | 6 |
Maintainer: Benilton S Carvalho <Benilton.Carvalho@cancer.org.uk>, Robert Scharpf <rscharpf@jhsph.edu>, Matt Ritchie <mritchie@wehi.EDU.AU> |
8 | 7 |
Description: Faster implementation of CRLMM specific to SNP 5.0 and 6.0 arrays, as well as a copy number tool specific to 5.0, 6.0, and Illumina platforms |
9 | 8 |
License: Artistic-2.0 |
10 |
-Depends: R (>= 2.13.0), |
|
11 |
- methods, |
|
12 |
- Biobase (>= 2.15.0), |
|
13 |
- oligoClasses (>= 1.17.34) |
|
14 |
-Imports: affyio (>= 1.19.2), |
|
9 |
+Depends: R (>= 2.14.0) |
|
10 |
+Imports: methods, |
|
11 |
+ Biobase (>= 2.15.4), |
|
12 |
+ oligoClasses (>= 1.17.36), |
|
13 |
+ BiocGenerics, |
|
14 |
+ affyio (>= 1.23.2), |
|
15 | 15 |
ellipse, |
16 |
- genefilter (>= 1.33.0), |
|
16 |
+ genefilter (>= 1.37.1), |
|
17 | 17 |
mvtnorm, |
18 |
- preprocessCore (>= 1.13.4), |
|
18 |
+ preprocessCore (>= 1.17.7), |
|
19 | 19 |
splines, |
20 | 20 |
stats, |
21 | 21 |
SNPchip, |
... | ... |
@@ -41,10 +41,8 @@ Collate: AllGenerics.R |
41 | 41 |
crlmmGT2.R |
42 | 42 |
crlmm-illumina.R |
43 | 43 |
snprma-functions.R |
44 |
- cn-functions.R |
|
45 | 44 |
utils.R |
46 | 45 |
zzz.R |
47 | 46 |
test_crlmm_package.R |
48 | 47 |
LazyLoad: yes |
49 | 48 |
biocViews: Microarray, Preprocessing, SNP, Bioinformatics, CopyNumberVariants |
50 |
-Packaged: 2011-04-30 17:10:59 UTC; biocbuild |
|
51 | 49 |
\ No newline at end of file |
... | ... |
@@ -1,56 +1,68 @@ |
1 | 1 |
useDynLib("crlmm", .registration=TRUE) |
2 | 2 |
|
3 |
-##--------------------------------------------------------------------------- |
|
4 |
-## Biobase |
|
5 |
-##--------------------------------------------------------------------------- |
|
6 |
-importClassesFrom(Biobase, AnnotatedDataFrame, AssayData, eSet, SnpSet, |
|
7 |
- NChannelSet, MIAME, Versioned, VersionedBiobase, |
|
8 |
- Versions) |
|
9 |
-importMethodsFrom(Biobase, annotation, "annotation<-", |
|
10 |
- annotatedDataFrameFrom, assayData, "assayData<-", |
|
11 |
- combine, dims, experimentData, "experimentData<-", |
|
12 |
- fData, featureData, "featureData<-", featureNames, |
|
13 |
- fvarMetadata, fvarLabels, pData, "pData<-", phenoData, |
|
14 |
- "phenoData<-", protocolData, "protocolData<-", |
|
15 |
- pubMedIds, rowMedians, sampleNames, snpCall, |
|
16 |
- snpCallProbability, |
|
17 |
- "snpCall<-", "snpCallProbability<-", storageMode, |
|
18 |
- "storageMode<-", updateObject, varLabels) |
|
19 |
-importFrom(Biobase, assayDataElement, assayDataElementNames, |
|
20 |
- assayDataElementReplace, assayDataNew, classVersion, |
|
21 |
- validMsg) |
|
22 |
- |
|
23 |
-##--------------------------------------------------------------------------- |
|
24 |
-## oligoClasses |
|
25 |
-##--------------------------------------------------------------------------- |
|
26 |
-importClassesFrom(oligoClasses, SnpSuperSet, AlleleSet, CNSet)##, ff_or_matrix) |
|
27 |
-importMethodsFrom(oligoClasses, allele, calls, "calls<-", confs, |
|
28 |
- "confs<-", cnConfidence, "cnConfidence<-", isSnp, |
|
29 |
- chromosome, position, A, B, |
|
30 |
- "A<-", "B<-", open, close, flags, |
|
31 |
- openff, closeff, |
|
32 |
- batchStatistics, "batchStatistics<-", updateObject, |
|
33 |
- checkOrder) |
|
34 |
- |
|
35 |
-importFrom(oligoClasses, chromosome2integer, celfileDate, list.celfiles, |
|
36 |
- copyNumber, initializeBigMatrix, initializeBigVector, isPackageLoaded) |
|
37 |
- |
|
38 |
- |
|
39 |
-importFrom(graphics, abline, axis, layout, legend, mtext, par, plot, |
|
40 |
- polygon, rect, segments, text, points, boxplot, lines) |
|
41 |
- |
|
42 |
-importFrom(lattice, xyplot, simpleKey, panel.grid, panel.xyplot, lrect, ltext, |
|
43 |
- lpoints, panel.number, lpolygon) |
|
44 |
- |
|
45 |
-importFrom(grDevices, grey) |
|
46 |
-importFrom(affyio, read.celfile.header, read.celfile) |
|
47 |
-importFrom(preprocessCore, normalize.quantiles.use.target, normalize.quantiles) |
|
48 |
-importFrom(utils, data, packageDescription, setTxtProgressBar, txtProgressBar) |
|
49 |
-importFrom(stats, coef, cov, dnorm, kmeans, lm, mad, median, quantile, sd, update) |
|
50 |
-importFrom(genefilter, rowSds) |
|
51 |
-importFrom(mvtnorm, dmvnorm) |
|
3 |
+importClassesFrom(Biobase, AssayData, eSet) |
|
4 |
+ |
|
5 |
+importClassesFrom(methods, ANY, character, formula, integer, list, |
|
6 |
+ matrix, oldClass) |
|
7 |
+importFrom(methods, setOldClass) |
|
8 |
+ |
|
9 |
+## importClassesFrom(oligoClasses, CNSet, CNSetLM, ff_matrix, |
|
10 |
+## ff_or_matrix, oligoSnpSet) |
|
11 |
+importClassesFrom(oligoClasses, CNSet, oligoSnpSet, ff_or_matrix) |
|
12 |
+##setOldClass(ff_or_matrix) |
|
13 |
+ |
|
14 |
+importMethodsFrom(Biobase, annotatedDataFrameFrom, annotation, |
|
15 |
+ assayData, experimentData, featureData, |
|
16 |
+ "featureData<-", featureNames, "featureNames<-", |
|
17 |
+ pData, "pData<-", phenoData, "phenoData<-", |
|
18 |
+ protocolData, "protocolData<-", rowMedians, |
|
19 |
+ sampleNames, snpCall, "snpCall<-", |
|
20 |
+ snpCallProbability, "snpCallProbability<-", |
|
21 |
+ storageMode, "storageMode<-", varLabels) |
|
22 |
+ |
|
23 |
+importMethodsFrom(BiocGenerics, cbind, colnames, Filter, get, |
|
24 |
+ intersect, lapply, ncol, NCOL, nrow, NROW, order, |
|
25 |
+ paste, pmax, pmin, rbind, rownames, sapply, setdiff, |
|
26 |
+ table, union, unique) |
|
27 |
+ |
|
28 |
+importMethodsFrom(genefilter, show) |
|
29 |
+ |
|
30 |
+importMethodsFrom(oligoClasses, A, "A<-", B, batch, batchNames, |
|
31 |
+ batchStatistics, "batchStatistics<-", calls, |
|
32 |
+ chromosome, close, confs, flags, isSnp, mean, nu, |
|
33 |
+ open, phi, "sampleNames<-") |
|
34 |
+ |
|
35 |
+importFrom(affyio, read.celfile, read.celfile.header) |
|
36 |
+ |
|
37 |
+importFrom(Biobase, assayDataElement, assayDataElementReplace, |
|
38 |
+ assayDataNew, copyEnv) |
|
39 |
+ |
|
52 | 40 |
importFrom(ellipse, ellipse) |
53 |
-##importFrom(ff, ffdf, physical.ff, physical.ffdf, ffrowapply) |
|
41 |
+ |
|
42 |
+importFrom(genefilter, rowSds) |
|
43 |
+ |
|
44 |
+importFrom(lattice, lpolygon, panel.grid, panel.number, panel.xyplot, |
|
45 |
+ xyplot) |
|
46 |
+ |
|
47 |
+importFrom(methods, "@<-", as, callNextMethod, is, new, validObject) |
|
48 |
+ |
|
49 |
+importFrom(mvtnorm, rmvnorm) |
|
50 |
+ |
|
51 |
+importFrom(oligoClasses, celfileDate, chromosome2integer, i2p, |
|
52 |
+ initializeBigMatrix, initializeBigVector, integerMatrix, |
|
53 |
+ isPackageLoaded, |
|
54 |
+ ldPath, ocLapply, ocProbesets, ocSamples, |
|
55 |
+ parStatus, |
|
56 |
+ splitIndicesByLength, splitIndicesByNode) |
|
57 |
+ |
|
58 |
+importFrom(preprocessCore, normalize.quantiles, |
|
59 |
+ normalize.quantiles.use.target) |
|
60 |
+ |
|
61 |
+importFrom(stats, coef, cov, dnorm, kmeans, lm, mad, median, quantile, |
|
62 |
+ sd) |
|
63 |
+ |
|
64 |
+importFrom(utils, packageDescription, setTxtProgressBar, |
|
65 |
+ txtProgressBar) |
|
54 | 66 |
|
55 | 67 |
##---------------------------------------------------------------------------- |
56 | 68 |
## export |
5 | 6 |
deleted file mode 100644 |
... | ... |
@@ -1,117 +0,0 @@ |
1 |
-copynumber <- function(filenames, |
|
2 |
- batch, |
|
3 |
- summaries=c("lrr", "baf", "ca", "cb", "gt", "gtconf"), |
|
4 |
- write=FALSE, |
|
5 |
- outdir=".", |
|
6 |
- onefile=FALSE, |
|
7 |
- rda=TRUE){ |
|
8 |
- fnamelist <- split(filenames, batch) |
|
9 |
- batchlist <- split(batch, batchlist) |
|
10 |
- ## for each element in list. Avoid |
|
11 |
- foreach(i=fnamelist) %do% simpleusage(filenames=fnameslist[[i]], batch=batchlist[[i]]) |
|
12 |
-} |
|
13 |
- |
|
14 |
-simpleusage <- function(filenames, batch, ...){ |
|
15 |
- object <- genotype(fnamelist, batch, ...) |
|
16 |
- ## run genotype on each element in fnamelist |
|
17 |
- object <- genotypeSummary(object) |
|
18 |
- ## different than currently implemented |
|
19 |
- ## - shrink to saved batch medians, etc. |
|
20 |
- ## - shinkage across markers (?) |
|
21 |
- object <- shrinkSummary(object) |
|
22 |
- results <- array(NA, nrow(object), ncol(object), length(summaries)) |
|
23 |
- if(c("ca", "cb") %in% summaries){ |
|
24 |
- ## estimateCnParameters |
|
25 |
- ##results <- estimatecnParameters(...) |
|
26 |
- } else results <- NULL |
|
27 |
- is.baf <- "baf" %in% summaries || "lrr" %in% summaries |
|
28 |
- if(is.baf){ |
|
29 |
- results2 <- calculateRtheta(object) ## returns list |
|
30 |
- ## keep as list, or coerce to array |
|
31 |
- } |
|
32 |
- if(write){ |
|
33 |
- ##write2(, onefile=onefile) |
|
34 |
- } |
|
35 |
- return(results) |
|
36 |
-} |
|
37 |
- |
|
38 |
- |
|
39 |
-imputeTheta <- function(ia, ib, theta, S=100){ |
|
40 |
-## y <- rbind(y1, y2) |
|
41 |
- y <- theta |
|
42 |
- N <- nrow(y); p <- ncol(y) |
|
43 |
- |
|
44 |
- mu0 <- apply(y, 2, mean, na.rm=TRUE) |
|
45 |
- sd0 <- mu0/5 |
|
46 |
- L0 <- matrix(.1, p, p); diag(L0) <- 1; L0 <- L0*outer(sd0, sd0) |
|
47 |
- nu0 <- p+2; S0 <- L0 |
|
48 |
- |
|
49 |
- ## starting values |
|
50 |
- Sigma <- S0 |
|
51 |
- Y.full <- y |
|
52 |
- O <- 1*(!is.na(y)) |
|
53 |
- if(all(O[1, ] == 0)){ |
|
54 |
- return(rep(NA, length(ia))) |
|
55 |
- } |
|
56 |
- for(j in seq_len(p)) Y.full[is.na(Y.full[, j]), j] <- mu0[j] |
|
57 |
- ## Gibbs sampler |
|
58 |
- ##THETA <- SIGMA <- Y.MISS <- NULL |
|
59 |
- ## problems: approx. 90 observations for the means in the other batches |
|
60 |
- ## -- only 3 observations for the mean in the current batch. |
|
61 |
-## THETA <- matrix(NA, iter, p) |
|
62 |
-## SIGMA <- matrix(NA, iter, p^2) |
|
63 |
- Y.MISS <- matrix(NA, S, sum(O[1,]==0)) |
|
64 |
- bafs <- matrix(NA, S, length(a)) |
|
65 |
- for(s in seq_len(S)){ |
|
66 |
- ## update lambda |
|
67 |
- ybar <- apply(Y.full, 2, mean) |
|
68 |
- Ln <- solve(solve(L0) + N*solve(Sigma)) |
|
69 |
- mun <- Ln%*%(solve(L0)%*%mu0 + N*solve(Sigma)%*%ybar) |
|
70 |
- lambda <- rmvnorm(1, mun, Ln) |
|
71 |
- |
|
72 |
- ##update sigma |
|
73 |
- Sn <- S0+(t(Y.full)-c(lambda))%*%t(t(Y.full)-c(lambda)) |
|
74 |
- Sigma <- solve(rwish(nu0+N, solve(Sn))) |
|
75 |
- |
|
76 |
- ## update missing data (only care about the first row.) |
|
77 |
- a <- O[1, ] == 1 |
|
78 |
- b <- O[1, ] == 0 |
|
79 |
- iSa <- solve(Sigma[a,a]) |
|
80 |
- beta.j <- Sigma[b,a]%*%iSa |
|
81 |
- Sigma.j <- Sigma[b,b]-Sigma[b,a]%*%iSa%*%Sigma[a,b] |
|
82 |
- lambda.j <- lambda[b]+beta.j%*%(t(Y.full[1,a, drop=FALSE]-lambda[a])) |
|
83 |
- Y.full[1,b] <- rmvnorm(1, lambda.j, Sigma.j) |
|
84 |
- ##SIGMA[s, ] <- c(Sigma) |
|
85 |
- Y.MISS[s, ] <- Y.full[1, O[1, ]==0] |
|
86 |
- } |
|
87 |
- na.cols <- which(is.na(y[1, ])) |
|
88 |
- THETA <- matrix(NA, S, 3) |
|
89 |
- THETA[, na.cols] <- Y.MISS |
|
90 |
- if(length(na.cols) < length(ia)) |
|
91 |
- THETA[, -na.cols] <- y[1, -na.cols] |
|
92 |
- |
|
93 |
- obs.theta <- atan2(ib, ia)*2/pi |
|
94 |
- theta.aa <- matrix(THETA[, 1], S, 3) |
|
95 |
- theta.ab <- matrix(THETA[, 2], S, 3) |
|
96 |
- theta.bb <- matrix(THETA[, 3], S, 3) |
|
97 |
- lessAA <- obs.theta < theta.aa |
|
98 |
- lessAB <- obs.theta < theta.ab |
|
99 |
- lessBB <- obs.theta < theta.bb |
|
100 |
- grAA <- !lessAA ## >= theta.aa |
|
101 |
- grAB <- !lessAB ## >= theta.ab |
|
102 |
- grBB <- !lessBB ## >= theta.bb |
|
103 |
- ##not.na <- !is.na(theta.aa) |
|
104 |
- I1 <- grAA & lessAB |
|
105 |
- I2 <- grAB & lessBB |
|
106 |
- ##mu <- apply(Y.MISS, 2, mean) |
|
107 |
- ##bf <- matrix(NA, S, length(ia)) |
|
108 |
- bf <- rep(NA, S*length(ia)) |
|
109 |
- I1 <- as.logical(I1); I2 <- as.logical(I2) |
|
110 |
- bf[I1] <- 0.5 * as.numeric(((obs.theta-theta.aa)/(theta.ab-theta.aa)))[I1] |
|
111 |
- bf[I2] <- as.numeric((.5 * (obs.theta - theta.ab) / (theta.bb - theta.ab)))[I2] + 0.5 |
|
112 |
- bf[as.logical(lessAA)] <- 0 |
|
113 |
- bf[as.logical(grBB)] <- 1 |
|
114 |
- bf <- matrix(bf, S, 3, byrow=TRUE) |
|
115 |
- pm.bf <- apply(bf, 2, mean) |
|
116 |
- return(pm.bf) |
|
117 |
-} |
... | ... |
@@ -1896,85 +1896,85 @@ isCall <- function(G, call){ |
1896 | 1896 |
G.call |
1897 | 1897 |
} |
1898 | 1898 |
|
1899 |
-computeSummary <- function(object, G.call, call, I, allele, Ns, index){ |
|
1900 |
- k <- match("grandMean", batchNames(object)) |
|
1901 |
- stats <- summaryStats(G.call, I, FUNS=c("rowMedians", "rowMAD")) |
|
1902 |
- Ns[, k] <- rowSums(G.call, na.rm=TRUE) |
|
1903 |
- updateStats(stats, Ns, object, call, allele, index) |
|
1904 |
- gc() |
|
1905 |
- return() |
|
1906 |
-} |
|
1899 |
+##computeSummary <- function(object, G.call, call, I, allele, Ns, index){ |
|
1900 |
+## k <- match("grandMean", batchNames(object)) |
|
1901 |
+## stats <- summaryStats(G.call, I, FUNS=c("rowMedians", "rowMAD")) |
|
1902 |
+## Ns[, k] <- rowSums(G.call, na.rm=TRUE) |
|
1903 |
+## updateStats(stats, Ns, object, call, allele, index) |
|
1904 |
+## gc() |
|
1905 |
+## return() |
|
1906 |
+##} |
|
1907 | 1907 |
|
1908 |
-updateTau <- function(object, tau2, G.call, call, I, allele, index){ |
|
1909 |
- k <- match("grandMean", batchNames(object)) |
|
1910 |
- logI <- log2(I) |
|
1911 |
- rm(I); gc() |
|
1912 |
- tau2[, k] <- summaryStats(G.call, logI, FUNS="rowMAD")^2 |
|
1913 |
- if(call==1 & allele=="A"){ |
|
1914 |
- tau2A.AA(object)[index, ] <- tau2 |
|
1915 |
- } |
|
1916 |
- if(call==1 & allele=="B"){ |
|
1917 |
- tau2B.AA(object)[index, ] <- tau2 |
|
1918 |
- } |
|
1919 |
- if(call==3 & allele=="A"){ |
|
1920 |
- tau2A.BB(object)[index, ] <- tau2 |
|
1921 |
- } |
|
1922 |
- if(call==3 & allele=="B"){ |
|
1923 |
- tau2B.BB(object)[index, ] <- tau2 |
|
1924 |
- } |
|
1925 |
- NULL |
|
1926 |
-} |
|
1908 |
+##updateTau <- function(object, tau2, G.call, call, I, allele, index){ |
|
1909 |
+## k <- match("grandMean", batchNames(object)) |
|
1910 |
+## logI <- log2(I) |
|
1911 |
+## rm(I); gc() |
|
1912 |
+## tau2[, k] <- summaryStats(G.call, logI, FUNS="rowMAD")^2 |
|
1913 |
+## if(call==1 & allele=="A"){ |
|
1914 |
+## tau2A.AA(object)[index, ] <- tau2 |
|
1915 |
+## } |
|
1916 |
+## if(call==1 & allele=="B"){ |
|
1917 |
+## tau2B.AA(object)[index, ] <- tau2 |
|
1918 |
+## } |
|
1919 |
+## if(call==3 & allele=="A"){ |
|
1920 |
+## tau2A.BB(object)[index, ] <- tau2 |
|
1921 |
+## } |
|
1922 |
+## if(call==3 & allele=="B"){ |
|
1923 |
+## tau2B.BB(object)[index, ] <- tau2 |
|
1924 |
+## } |
|
1925 |
+## NULL |
|
1926 |
+##} |
|
1927 | 1927 |
|
1928 |
-updateCors <- function(cors, G.call, I){ |
|
1929 |
- k <- match("grandMean", batchNames(object)) |
|
1930 |
- cors[, k] <- rowCors(I[[1]]*G.call, I[[2]]*G.call, na.rm=TRUE) |
|
1931 |
- if(call==1){ |
|
1932 |
- corrAA(object)[index, ] <- cors |
|
1933 |
- } |
|
1934 |
- if(call==2){ |
|
1935 |
- corrAB(object)[index, ] <- cors |
|
1936 |
- } |
|
1937 |
- if(call==3){ |
|
1938 |
- corrBB(object)[index, ] <- cors |
|
1939 |
- } |
|
1940 |
-} |
|
1928 |
+##updateCors <- function(cors, G.call, I){ |
|
1929 |
+## k <- match("grandMean", batchNames(object)) |
|
1930 |
+## cors[, k] <- rowCors(I[[1]]*G.call, I[[2]]*G.call, na.rm=TRUE) |
|
1931 |
+## if(call==1){ |
|
1932 |
+## corrAA(object)[index, ] <- cors |
|
1933 |
+## } |
|
1934 |
+## if(call==2){ |
|
1935 |
+## corrAB(object)[index, ] <- cors |
|
1936 |
+## } |
|
1937 |
+## if(call==3){ |
|
1938 |
+## corrBB(object)[index, ] <- cors |
|
1939 |
+## } |
|
1940 |
+##} |
|
1941 | 1941 |
|
1942 |
-updateStats <- function(stats, Ns, object, call, allele, tau2, index){ |
|
1943 |
- if(call==1){ |
|
1944 |
- Ns.AA(object)[index, ] <- Ns |
|
1945 |
- if(allele=="A"){ |
|
1946 |
- medianA.AA(object)[index, k] <- stats[, 1] |
|
1947 |
- madA.AA(object)[index, k] <- stats[, 2] |
|
1948 |
- updateTau(object, tau2, G.call, call, I, allele, index) |
|
1949 |
- } else { |
|
1950 |
- medianB.AA(object)[index, k] <- stats[, 1] |
|
1951 |
- madB.AA(object)[index, k] <- stats[, 2] |
|
1952 |
- updateTau(object, tau2, G.call, call, I, allele, index) |
|
1953 |
- } |
|
1954 |
- } |
|
1955 |
- if(call==2){ |
|
1956 |
- Ns.AB(object)[index, ] <- Ns |
|
1957 |
- if(allele=="A"){ |
|
1958 |
- medianA.AB(object)[index, k] <- stats[, 1] |
|
1959 |
- madA.AB(object)[index, k] <- stats[, 2] |
|
1960 |
- } else { |
|
1961 |
- medianB.AB(object)[index, k] <- stats[, 1] |
|
1962 |
- madB.AB(object)[index, k] <- stats[, 2] |
|
1963 |
- } |
|
1964 |
- } |
|
1965 |
- if(call==3){ |
|
1966 |
- Ns.BB(object)[index, ] <- Ns |
|
1967 |
- if(allele=="A"){ |
|
1968 |
- medianA.BB(object)[index, k] <- stats[, 1] |
|
1969 |
- madA.BB(object)[index, k] <- stats[, 2] |
|
1970 |
- updateTau(object, tau2, G.call, call, I, allele, index) |
|
1971 |
- } else { |
|
1972 |
- medianB.BB(object)[index, k] <- stats[, 1] |
|
1973 |
- madB.BB(object)[index, k] <- stats[, 2] |
|
1974 |
- updateTau(object, tau2, G.call, call, I, allele, index) |
|
1975 |
- } |
|
1976 |
- } |
|
1977 |
-} |
|
1942 |
+##updateStats <- function(stats, Ns, object, call, allele, tau2, index){ |
|
1943 |
+## if(call==1){ |
|
1944 |
+## Ns.AA(object)[index, ] <- Ns |
|
1945 |
+## if(allele=="A"){ |
|
1946 |
+## medianA.AA(object)[index, k] <- stats[, 1] |
|
1947 |
+## madA.AA(object)[index, k] <- stats[, 2] |
|
1948 |
+## updateTau(object, tau2, G.call, call, I, allele, index) |
|
1949 |
+## } else { |
|
1950 |
+## medianB.AA(object)[index, k] <- stats[, 1] |
|
1951 |
+## madB.AA(object)[index, k] <- stats[, 2] |
|
1952 |
+## updateTau(object, tau2, G.call, call, I, allele, index) |
|
1953 |
+## } |
|
1954 |
+## } |
|
1955 |
+## if(call==2){ |
|
1956 |
+## Ns.AB(object)[index, ] <- Ns |
|
1957 |
+## if(allele=="A"){ |
|
1958 |
+## medianA.AB(object)[index, k] <- stats[, 1] |
|
1959 |
+## madA.AB(object)[index, k] <- stats[, 2] |
|
1960 |
+## } else { |
|
1961 |
+## medianB.AB(object)[index, k] <- stats[, 1] |
|
1962 |
+## madB.AB(object)[index, k] <- stats[, 2] |
|
1963 |
+## } |
|
1964 |
+## } |
|
1965 |
+## if(call==3){ |
|
1966 |
+## Ns.BB(object)[index, ] <- Ns |
|
1967 |
+## if(allele=="A"){ |
|
1968 |
+## medianA.BB(object)[index, k] <- stats[, 1] |
|
1969 |
+## madA.BB(object)[index, k] <- stats[, 2] |
|
1970 |
+## updateTau(object, tau2, G.call, call, I, allele, index) |
|
1971 |
+## } else { |
|
1972 |
+## medianB.BB(object)[index, k] <- stats[, 1] |
|
1973 |
+## madB.BB(object)[index, k] <- stats[, 2] |
|
1974 |
+## updateTau(object, tau2, G.call, call, I, allele, index) |
|
1975 |
+## } |
|
1976 |
+## } |
|
1977 |
+##} |
|
1978 | 1978 |
|
1979 | 1979 |
crlmmCopynumber <- function(object, |
1980 | 1980 |
MIN.SAMPLES=10, |
... | ... |
@@ -14,7 +14,7 @@ getCor <- function(object, gt){ |
14 | 14 |
getTau2 <- function(object, gt){ |
15 | 15 |
bs <- batchStatistics(object) |
16 | 16 |
nma <- paste("tau2A", gt, sep=".") |
17 |
- nma <- paste("tau2B", gt, sep=".") |
|
17 |
+ nmb <- paste("tau2B", gt, sep=".") |
|
18 | 18 |
tau2.a <- assayDataElement(bs, nma) |
19 | 19 |
tau2.b <- assayDataElement(bs, nmb) |
20 | 20 |
cbind(tau2.a, tau2.b) |
... | ... |
@@ -11,6 +11,8 @@ cnSet2oligoSnpSet <- function(object){ |
11 | 11 |
is.lds <- ifelse(is(calls(object), "ff_matrix") | is(calls(object), "ffdf"), TRUE, FALSE) |
12 | 12 |
if(is.lds) stopifnot(isPackageLoaded("ff")) |
13 | 13 |
b.r <- calculateRBaf(object) |
14 |
+ b <- integerMatrix(b.r[[1]], 1000) |
|
15 |
+ r <- integerMatrix(b.r[[2]], 100) |
|
14 | 16 |
## if(is.lds){ |
15 | 17 |
## ## initialize a big matrix for raw copy number |
16 | 18 |
## message("creating an ff object for storing total copy number") |
... | ... |
@@ -438,36 +440,36 @@ setMethod("calculatePosteriorMean", signature(object="CNSet"), |
438 | 440 |
return(pm) |
439 | 441 |
}) |
440 | 442 |
|
441 |
- .bivariateCenter <- function(nu, phi){ |
|
442 |
- ## lexical scope for mus, CA, CB |
|
443 |
- if(CA <= 2 & CB <= 2 & (CA+CB) < 4){ |
|
444 |
- mus[,1, ] <- log2(nu[, 1, ] + CA * |
|
445 |
- phi[, 1, ]) |
|
446 |
- mus[,2, ] <- log2(nu[, 2, ] + CB * |
|
447 |
- phi[, 2, ]) |
|
448 |
- } else { ## CA > 2 |
|
449 |
- if(CA > 2){ |
|
450 |
- theta <- pi/4*Sigma[,2,] |
|
451 |
- shiftA <- CA/4*phi[, 1, ] * cos(theta) |
|
452 |
- shiftB <- CA/4*phi[, 1, ] * sin(theta) |
|
453 |
- mus[, 1, ] <- log2(nu[, 1, ] + 2 * phi[,1,]+shiftA) |
|
454 |
- mus[, 2, ] <- log2(nu[, 2, ] + CB *phi[,2,]+shiftB) |
|
455 |
- } |
|
456 |
- if(CB > 2){ |
|
457 |
- ## CB > 2 |
|
458 |
- theta <- pi/2-pi/4*Sigma[,2,] |
|
459 |
- shiftA <- CB/4*phi[, 2, ] * cos(theta) |
|
460 |
- shiftB <- CB/4*phi[, 2, ] * sin(theta) |
|
461 |
- mus[, 1, ] <- log2(nu[, 1, ] + CA*phi[,1,]+shiftA) |
|
462 |
- mus[, 2, ] <- log2(nu[, 2, ]+ 2*phi[,2,]+shiftB) |
|
463 |
- } |
|
464 |
- if(CA == 2 & CB == 2){ |
|
465 |
- mus[, 1, ] <- log2(nu[, 1, ] + 1/2*CA*phi[,1,]) |
|
466 |
- mus[, 2, ] <- log2(nu[, 2, ]+ 1/2*CB*phi[,2,]) |
|
467 |
- } |
|
468 |
- } |
|
469 |
- mus |
|
470 |
- } |
|
443 |
+##.bivariateCenter <- function(nu, phi){ |
|
444 |
+## ## lexical scope for mus, CA, CB |
|
445 |
+## if(CA <= 2 & CB <= 2 & (CA+CB) < 4){ |
|
446 |
+## mus[,1, ] <- log2(nu[, 1, ] + CA * |
|
447 |
+## phi[, 1, ]) |
|
448 |
+## mus[,2, ] <- log2(nu[, 2, ] + CB * |
|
449 |
+## phi[, 2, ]) |
|
450 |
+## } else { ## CA > 2 |
|
451 |
+## if(CA > 2){ |
|
452 |
+## theta <- pi/4*Sigma[,2,] |
|
453 |
+## shiftA <- CA/4*phi[, 1, ] * cos(theta) |
|
454 |
+## shiftB <- CA/4*phi[, 1, ] * sin(theta) |
|
455 |
+## mus[, 1, ] <- log2(nu[, 1, ] + 2 * phi[,1,]+shiftA) |
|
456 |
+## mus[, 2, ] <- log2(nu[, 2, ] + CB *phi[,2,]+shiftB) |
|
457 |
+## } |
|
458 |
+## if(CB > 2){ |
|
459 |
+## ## CB > 2 |
|
460 |
+## theta <- pi/2-pi/4*Sigma[,2,] |
|
461 |
+## shiftA <- CB/4*phi[, 2, ] * cos(theta) |
|
462 |
+## shiftB <- CB/4*phi[, 2, ] * sin(theta) |
|
463 |
+## mus[, 1, ] <- log2(nu[, 1, ] + CA*phi[,1,]+shiftA) |
|
464 |
+## mus[, 2, ] <- log2(nu[, 2, ]+ 2*phi[,2,]+shiftB) |
|
465 |
+## } |
|
466 |
+## if(CA == 2 & CB == 2){ |
|
467 |
+## mus[, 1, ] <- log2(nu[, 1, ] + 1/2*CA*phi[,1,]) |
|
468 |
+## mus[, 2, ] <- log2(nu[, 2, ]+ 1/2*CB*phi[,2,]) |
|
469 |
+## } |
|
470 |
+## } |
|
471 |
+## mus |
|
472 |
+## } |
|
471 | 473 |
|
472 | 474 |
## for a given copy number, return a named list of bivariate normal prediction regions |
473 | 475 |
## - elements of list are named by genotype |
... | ... |
@@ -214,7 +214,7 @@ setMethod("annotatedDataFrameFrom", "ffdf", Biobase:::annotatedDataFrameFromMatr |
214 | 214 |
|
215 | 215 |
## Document this... |
216 | 216 |
getBAF <- function(theta, canonicalTheta) |
217 |
- .Call('normalizeBAF', theta, ct) |
|
217 |
+ .Call('normalizeBAF', theta, canonicalTheta) |
|
218 | 218 |
|
219 | 219 |
|
220 | 220 |
validCEL <- function(celfiles){ |
... | ... |
@@ -11,7 +11,9 @@ |
11 | 11 |
|
12 | 12 |
} |
13 | 13 |
\details{ |
14 |
- This object was created from the copynumber vignette in inst/scripts. |
|
14 |
+ This object was created from the copynumber vignette in |
|
15 |
+ inst/scripts. A subset of markers was selected to keep the package |
|
16 |
+ size small. |
|
15 | 17 |
} |
16 | 18 |
\usage{ |
17 | 19 |
|