Browse code

bringing changes from GitHub; from now on, keeping code on BioC servers

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

Benilton Carvalho authored on 03/10/2011 18:12:55
Showing8 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.11.43
4
+Version: 1.11.46
5 5
 Date: 2010-12-10
6 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
7 7
 Maintainer: Benilton S Carvalho <Benilton.Carvalho@cancer.org.uk>, Robert Scharpf <rscharpf@jhsph.edu>, Matt Ritchie <mritchie@wehi.EDU.AU>
... ...
@@ -47,7 +47,3 @@ Collate: AllGenerics.R
47 47
 	 test_crlmm_package.R
48 48
 LazyLoad: yes
49 49
 biocViews: Microarray, Preprocessing, SNP, Bioinformatics, CopyNumberVariants
50
-<<<<<<< HEAD
51
-=======
52
-Packaged: 2011-04-30 17:10:59 UTC; biocbuild
53
->>>>>>> crlmm on github
... ...
@@ -118,11 +118,6 @@ construct <- function(filenames,
118 118
 constructAffy <- function(filenames, sns, cdfName, batch, verbose=TRUE){
119 119
 	stopifnot(!missing(filenames))
120 120
 	if(missing(sns)) sns <- basename(filenames)
121
-	stopifnot(length(sns)==length(batch))
122
-	is.lds <- ifelse(isPackageLoaded("ff"), TRUE, FALSE)
123
-	if(!is.lds) stop("this function now requires that the ff package be loaded")
124
-	if(missing(cdfName)) stop("must specify cdfName")
125
-	if(!isValidCdfName(cdfName)) stop("cdfName not valid.  see validCdfNames")
126 121
 	stopifnot(!missing(batch))
127 122
 	##---------------------------------------------------------------------------
128 123
 	##
... ...
@@ -219,6 +214,7 @@ genotype <- function(filenames,
219 214
 		     gender=NULL,
220 215
 		     returnParams=TRUE,
221 216
 		     badSNP=0.7){
217
+	stopifnot(require("ff"))
222 218
 	cnSet <- constructAffy(filenames=filenames,
223 219
 			       cdfName=cdfName,
224 220
 			       batch=batch, verbose=verbose)
... ...
@@ -211,3 +211,7 @@ loadObject <- function(filename, load.it){
211 211
 
212 212
 setMethod("annotatedDataFrameFrom", "ff_matrix", Biobase:::annotatedDataFrameFromMatrix)
213 213
 setMethod("annotatedDataFrameFrom", "ffdf", Biobase:::annotatedDataFrameFromMatrix)
214
+
215
+## Document this...
216
+getBAF <- function(theta, canonicalTheta)
217
+    .Call('normalizeBAF', theta, ct)
... ...
@@ -43,6 +43,8 @@ data(sample.CNSet)
43 43
 baf.lrr <- calculateRBaf(sample.CNSet, "SHELF")
44 44
 hist(baf.lrr[["baf"]], breaks=100)
45 45
 hist(baf.lrr[["lrr"]], breaks=100)
46
+
47
+
46 48
 }
47 49
 % Add one or more standard keywords, see file 'KEYWORDS' in the
48 50
 % R documentation directory.
... ...
@@ -1,49 +1,51 @@
1
-\name{readIDAT}
2
-\alias{readIDAT}
3
-\alias{readIDAT}
4
-
5
-\title{Low-level function to read idat files}
6
-
7
-\description{
8
-  Reads intensity information for each bead type from a single .idat file 
9
-  from Infinium II platforms.  This is a low-level function which \code{readIdatFiles} 
10
-  is a wrapper to}
11
-
12
-\usage{
13
-readIDAT(idatFile)
14
-}
15
-
16
-\arguments{
17
-  \item{idatFile}{character string specifying idat file to be read in}
18
-}
19
-
20
-\details{
21
-This function returns a list containing summarised intensities and other 
22
-information extracted from a single .idat file.
23
-
24
-The \code{readIdatFiles} function makes use of \code{readIDAT} to read 
25
-the paired Cy3 and Cy5 idats from one or more arrays.
26
-
27
-Thanks to Keith Baggerly who providing this code.
28
-}
29
-
30
-\value{
31
-  list which includes item \code{Quants} which contains average intensity (\code{Mean}),
32
-  number of beads (\code{NBeads}) and a measure of variability (\code{SD}) for 
33
-  each bead type on the array.
34
-}
35
-
36
-\references{
37
-  Ritchie ME, Carvalho BS, Hetrick KN, Tavar\'{e} S, Irizarry RA.
38
-  R/Bioconductor software for Illumina's Infinium whole-genome 
39
-  genotyping BeadChips. Bioinformatics. 2009 Oct 1;25(19):2621-3.
40
-}
41
-
42
-\author{Keith Baggerly, with modifications by Matt Ritchie}
43
-
44
-\examples{
45
-#idatdata = readIDAT("4019585367_A_Grn.idat")
46
-#names(idatdata
47
-#idatdata$Quants[1:5,]
48
-}
49
-\keyword{IO}
1
+\name{readIDAT}
2
+\alias{readIDAT}
3
+\alias{readIDAT}
4
+
5
+\title{Low-level function to read idat files}
6
+
7
+\description{
8
+  Reads intensity information for each bead type from a single .idat file
9
+  from Infinium II platforms.  This is a low-level function which \code{readIdatFiles}
10
+  is a wrapper to}
11
+
12
+\usage{
13
+readIDAT(idatFile)
14
+}
15
+
16
+\arguments{
17
+  \item{idatFile}{character string specifying idat file to be read in}
18
+}
19
+
20
+\details{
21
+
22
+This function returns a list containing summarised intensities and other
23
+information extracted from a single .idat file.
24
+
25
+The \code{readIdatFiles} function makes use of \code{readIDAT} to read
26
+the paired Cy3 and Cy5 idats from one or more arrays.
27
+
28
+Thanks to Keith Baggerly who providing this code.
29
+}
30
+
31
+\value{
32
+  list which includes item \code{Quants} which contains average intensity (\code{Mean}),
33
+  number of beads (\code{NBeads}) and a measure of variability (\code{SD}) for
34
+  each bead type on the array.
35
+}
36
+
37
+\references{
38
+  Ritchie ME, Carvalho BS, Hetrick KN, Tavar\'{e} S, Irizarry RA.
39
+  R/Bioconductor software for Illumina's Infinium whole-genome
40
+  genotyping BeadChips. Bioinformatics. 2009 Oct 1;25(19):2621-3.
41
+}
42
+
43
+\author{Keith Baggerly, with modifications by Matt Ritchie}
44
+
45
+\examples{
46
+#idatdata = readIDAT("4019585367_A_Grn.idat")
47
+#names(idatdata
48
+#idatdata$Quants[1:5,]
49
+}
50
+\keyword{IO}
51
+
... ...
@@ -10,3 +10,5 @@ SEXP gtypeCallerPart2(SEXP *, SEXP *, SEXP *, SEXP *,
10 10
                       SEXP *, SEXP *, SEXP *, SEXP *,
11 11
 		      SEXP *, SEXP *, SEXP *, SEXP *,
12 12
                       SEXP *, SEXP *, SEXP *);
13
+
14
+SEXP normalizeBAF(SEXP *, SEXP *);
... ...
@@ -5,6 +5,7 @@
5 5
 static const R_CallMethodDef CallEntries[] = {
6 6
     {"gtypeCallerPart1", (DL_FUNC)&gtypeCallerPart1, 17},
7 7
     {"gtypeCallerPart2", (DL_FUNC)&gtypeCallerPart2, 19},
8
+    {"normalizeBAF", (DL_FUNC)&normalizeBAF, 2},
8 9
     {NULL, NULL, 0}
9 10
 };
10 11
 
... ...
@@ -142,6 +142,54 @@ void mad_median(double *datavec, int *classvec, int class, double trim, int cols
142 142
   Free(buffer);
143 143
 }
144 144
 
145
+SEXP normalizeBAF(SEXP theta, SEXP cTheta){
146
+  /*
147
+    ARGUMENTS:
148
+    theta.: N x C matrix with estimated \theta
149
+    cTheta: N x 3 matrix with canonical \thetas (AA, AB, BB)
150
+
151
+    VALUE:
152
+    baf: N x C matrix with normalized \theta
153
+  */
154
+
155
+  SEXP baf;
156
+  double *p2baf, *p2theta, *p2ctheta;
157
+  int i, j, idx, rowsT, rowsCT, colsT, colsCT;
158
+  rowsT = INTEGER(getAttrib(theta, R_DimSymbol))[0];
159
+  rowsCT = INTEGER(getAttrib(cTheta, R_DimSymbol))[0];
160
+  if (rowsT != rowsCT)
161
+    error("Number of rows of 'theta' must match number of rows of 'cTheta'\n");
162
+  colsCT = INTEGER(getAttrib(cTheta, R_DimSymbol))[1];
163
+  if (colsCT != 3)
164
+    error("'cTheta' must have 3 columns: AA, AB and BB\n");
165
+  colsT = INTEGER(getAttrib(theta, R_DimSymbol))[1];
166
+
167
+  PROTECT(baf = allocMatrix(REALSXP, rowsT, colsT));
168
+  p2baf = NUMERIC_POINTER(baf);
169
+  p2theta = NUMERIC_POINTER(theta);
170
+  p2ctheta = NUMERIC_POINTER(cTheta);
171
+  for (i=0; i < rowsT; i++){
172
+    for (j=0; j < colsT; j++){
173
+      idx = i + j*rowsT;
174
+      if (ISNA(p2theta[idx]) || ISNA(p2ctheta[i]) || ISNA(p2ctheta[i+rowsT]) || ISNA(p2ctheta[i+2*rowsT])){
175
+	p2baf[idx] = NA_REAL;
176
+      }else if (p2theta[idx] < p2ctheta[i]){
177
+	p2baf[idx] = 0;
178
+      }else if (p2theta[idx] >= p2ctheta[i] & p2theta[idx] < p2ctheta[i + rowsT]){
179
+	p2baf[idx] = .5*(p2theta[idx]-p2ctheta[i])/(p2ctheta[i+rowsT]-p2ctheta[i]);
180
+      }else if(p2theta[idx] >= p2ctheta[i+rowsT] & p2theta[idx] < p2ctheta[i + 2*rowsT]){
181
+	p2baf[idx] = .5+.5*(p2theta[idx]-p2ctheta[i+rowsT])/(p2ctheta[i+2*rowsT]-p2ctheta[i+rowsT]);
182
+      }else{
183
+	p2baf[idx] = 1;
184
+      }
185
+    }
186
+  }
187
+
188
+  UNPROTECT(1);
189
+  return(baf);
190
+}
191
+
192
+
145 193
 /* Pieces below are for testing */
146 194
 
147 195
 static void mad_stats(double *data, double *m1, double *m2, double *m3, int *class, int rows, int cols, double *trim){