Browse code

crlmm on github

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

Rob Scharp authored on 04/10/2011 16:11:02
Showing10 changed files

... ...
@@ -30,12 +30,14 @@ importMethodsFrom(oligoClasses, allele, calls, "calls<-", confs,
30 30
 		  openff, closeff,
31 31
 		  batchStatistics, "batchStatistics<-", updateObject,
32 32
 		  order, checkOrder)
33
+
33 34
 importFrom(oligoClasses, chromosome2integer, celfileDate, list.celfiles,
34 35
            copyNumber, initializeBigMatrix, initializeBigVector, isPackageLoaded)
35 36
 
36 37
 
37 38
 importFrom(graphics, abline, axis, layout, legend, mtext, par, plot,
38 39
            polygon, rect, segments, text, points, boxplot, lines)
40
+
39 41
 importFrom(lattice, xyplot, simpleKey, panel.grid, panel.xyplot, lrect, ltext,
40 42
 	   lpoints, panel.number, lpolygon)
41 43
 
... ...
@@ -48,6 +50,7 @@ importFrom(genefilter, rowSds)
48 50
 importFrom(mvtnorm, dmvnorm)
49 51
 importFrom(ellipse, ellipse)
50 52
 importFrom(ff, ffdf, physical.ff, physical.ffdf, ffrowapply)
53
+
51 54
 importClassesFrom(oligoClasses, ff_matrix, ffdf)
52 55
 exportMethods(lines)
53 56
 exportMethods(CA, CB)
... ...
@@ -221,6 +221,7 @@ genotype <- function(filenames,
221 221
 	if(!(is(A(cnSet), "ff") || is(A(cnSet), "ffdf"))){
222 222
 		stop("The ff package is required for this function.")
223 223
 	}
224
+
224 225
 	mixtureParams <- snprmaAffy(cnSet, filenames=filenames,
225 226
 				    mixtureSampleSize=mixtureSampleSize,
226 227
 				    eps=eps,
... ...
@@ -271,7 +271,7 @@ readIDAT <- function(idatFile){
271 271
                  MostlyNull <- readChar(tempCon, nChars)
272 272
                  MostlyNull
273 273
              },
274
-             "Barcode" = {                 
274
+             "Barcode" = {
275 275
                  seek(tempCon, fields["Barcode", "Byte Offset"])
276 276
                  nChars <- readBin(tempCon, "integer", n=1, size=1, signed=FALSE)
277 277
                  Barcode <- readChar(tempCon, nChars)
... ...
@@ -334,7 +334,7 @@ readIDAT <- function(idatFile){
334 334
 
335 335
   readFields <- setdiff(rownames(fields), "nSNPsRead")
336 336
   names(readFields) <- readFields
337
-  
337
+
338 338
   allFields <- lapply(readFields, readBlock)
339 339
 
340 340
   close(tempCon)
... ...
@@ -350,7 +350,7 @@ readIDAT <- function(idatFile){
350 350
 
351 351
   InfoNames <- c("MidBlock", "RunInfo", "RedGreen", "Barcode", "ChipType")
352 352
   Info <- allFields[intersect(names(allFields), InfoNames)]
353
-  
353
+
354 354
   idatValues <-
355 355
     list(fileSize=fileSize,
356 356
          versionNumber=versionNumber,
... ...
@@ -20,7 +20,6 @@ linearParamElementReplace <- function(obj, elt, value) {
20 20
 }
21 21
 
22 22
 ## parameters
23
-
24 23
 ## allele A
25 24
 ##   autosome SNPs
26 25
 ##   autosome NPs
... ...
@@ -215,3 +215,5 @@ setMethod("annotatedDataFrameFrom", "ffdf", Biobase:::annotatedDataFrameFromMatr
215 215
 ## Document this...
216 216
 getBAF <- function(theta, canonicalTheta)
217 217
     .Call('normalizeBAF', theta, ct)
218
+
219
+
... ...
@@ -99,8 +99,6 @@ if (require(genomewidesnp6Crlmm) & require(hapmapsnp6)){
99 99
   gender <- c("female", "female", "male")
100 100
   gender[gender == "female"] <- 2
101 101
   gender[gender == "male"] <- 1
102
-  \dontrun{  (crlmmOutput <- crlmm(cels, gender=gender))}
103
-
104 102
   \dontrun{
105 103
 	  ## we need the A and B intensities
106 104
 	  stopifnot(require("ff"))
... ...
@@ -36,14 +36,12 @@ data(sample.CNSet2)
36 36
       marker.index <- marker.index[1:60e3]
37 37
       sample.CNSet <- hapmapSet[marker.index, c(1168:1169)]
38 38
       ## 2 samples, many markers
39
-      save(sample.CNSet, file="~/Software/crlmm/data/sample.CNSet.rda")
40 39
 
41 40
       ## all samples, a few markers
42 41
       snp.index <- which(isSnp(hapmapSet))[1:100]
43 42
       np.index <- which(!isSnp(hapmapSet))[1:100]
44 43
       marker.index <- c(snp.index, np.index)
45 44
       sample.CNSet2 <- hapmapSet[marker.index, ]
46
-      save(sample.CNSet2, file="~/Software/crlmm/data/sample.CNSet2.rda")
47 45
 }
48 46
 data(sample.CNSet)
49 47
 ## --------------------------------------------------
... ...
@@ -92,10 +90,6 @@ phi(sample.CNSet, "A")[1:5, ]
92 90
 (ca <- CA(sample.CNSet, i=1:5, j=1:2))
93 91
 ## copy number for allele B, first 5 markers, first 2 samples
94 92
 (cb <- CB(sample.CNSet, i=1:5, j=1:2))
95
-## total copy number for first 5 markers, first 2 samples
96
-(cn1 <- ca+cb)
97
-
98
-## total copy number at first 5 nonpolymorphic loci
99 93
 index <- which(!isSnp(sample.CNSet))[1:5]
100 94
 cn2 <- CA(sample.CNSet, i=index, j=1:2)
101 95
 ## note, cb is NA at nonpolymorphic loci
... ...
@@ -112,5 +106,23 @@ all.equal(cn4, cn2)
112 106
 cn5 <- totalCopynumber(sample.CNSet, i=1:5)
113 107
 ## all markers, samples 1-2
114 108
 cn6 <- totalCopynumber(sample.CNSet, j=1:2)
109
+=======
110
+index <- which(!isSnp(cnSet))[1:5]
111
+cn2 <- CA(cnSet, i=index, j=1:2)
112
+## note, cb is NA at nonpolymorphic loci
113
+(cb <- CB(cnSet, i=index, j=1:2))
114
+## note, ca+cb will give NAs at nonpolymorphic loci
115
+CA(cnSet, i=index, j=1:2) + cb
116
+## A shortcut for total copy number
117
+cn3 <- totalCopynumber(cnSet, i=1:5, j=1:2)
118
+all.equal(cn3, cn1)
119
+cn4 <- totalCopynumber(cnSet, i=index, j=1:2)
120
+all.equal(cn4, cn2)
121
+
122
+## markers 1-5, all samples
123
+cn5 <- totalCopynumber(cnSet, i=1:5)
124
+## all markers, samples 1-5
125
+cn6 <- totalCopynumber(cnSet, j=1:5)
126
+>>>>>>> crlmm on github
115 127
 }
116 128
 \keyword{datasets}
... ...
@@ -12,3 +12,4 @@ SEXP gtypeCallerPart2(SEXP *, SEXP *, SEXP *, SEXP *,
12 12
                       SEXP *, SEXP *, SEXP *);
13 13
 
14 14
 SEXP normalizeBAF(SEXP *, SEXP *);
15
+
... ...
@@ -10,5 +10,5 @@ static const R_CallMethodDef CallEntries[] = {
10 10
 };
11 11
 
12 12
 void R_init_crlmm(DllInfo *dll){
13
-    R_registerRoutines(dll, NULL, CallEntries, NULL, NULL);    
13
+    R_registerRoutines(dll, NULL, CallEntries, NULL, NULL);
14 14
 }
... ...
@@ -76,7 +76,7 @@ void trimmed_stats(double *data, double *m1, double *m2, double *m3, int *class,
76 76
     n1=0;
77 77
     n2=0;
78 78
     n3=0;
79
-    
79
+
80 80
     for (j=0; j < cols; j++){
81 81
       if (class[j*rows + i] == 1){
82 82
 	datvec[j]=data[j*rows + i];
... ...
@@ -201,7 +201,7 @@ static void mad_stats(double *data, double *m1, double *m2, double *m3, int *cla
201 201
     n1=0;
202 202
     n2=0;
203 203
     n3=0;
204
-    
204
+
205 205
     for (j=0; j < cols; j++){
206 206
       if (class[j*rows + i] == 1){
207 207
 	datvec[j]=data[j*rows + i];
... ...
@@ -239,7 +239,7 @@ SEXP test_mad_median(SEXP X, SEXP Y, SEXP trim){
239 239
   PROTECT(dim1 = getAttrib(X,R_DimSymbol));
240 240
   rows = INTEGER(dim1)[0];
241 241
   cols = INTEGER(dim1)[1];
242
-  
242
+
243 243
   Xptr = NUMERIC_POINTER(AS_NUMERIC(X));
244 244
   Yptr = INTEGER_POINTER(AS_INTEGER(Y));
245 245
   Tptr = NUMERIC_POINTER(AS_NUMERIC(trim));
... ...
@@ -247,11 +247,11 @@ SEXP test_mad_median(SEXP X, SEXP Y, SEXP trim){
247 247
   PROTECT(estimates1 = allocMatrix(REALSXP, rows, 3));
248 248
   PROTECT(estimates2 = allocMatrix(REALSXP, rows, 3));
249 249
   PROTECT(estimates3 = allocMatrix(REALSXP, rows, 3));
250
-  
250
+
251 251
   Mptr1 = NUMERIC_POINTER(estimates1);
252 252
   Mptr2 = NUMERIC_POINTER(estimates2);
253 253
   Mptr3 = NUMERIC_POINTER(estimates3);
254
-  
254
+
255 255
   mad_stats(Xptr, Mptr1, Mptr2, Mptr3, Yptr, rows, cols, Tptr);
256 256
 
257 257
   PROTECT(output = allocVector(VECSXP,3));
... ...
@@ -260,6 +260,6 @@ SEXP test_mad_median(SEXP X, SEXP Y, SEXP trim){
260 260
   SET_VECTOR_ELT(output, 2, estimates3);
261 261
 
262 262
   UNPROTECT(5);
263
-  
263
+
264 264
   return output;
265 265
 }