Browse code

move code for checking genotypes from crlmm doc to genotype doc

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

Rob Scharp authored on 11/10/2011 12:18:27
Showing 4 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.58
4
+Version: 1.11.60
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>
... ...
@@ -35,8 +35,7 @@ data(cnSetExample2)
35 35
       data(hapmapSet, package="CnvScripts")
36 36
       marker.index <- which(chromosome(hapmapSet) == 8)
37 37
       marker.index <- marker.index[1:60e3]
38
-      sample.CNSet <- hapmapSet[marker.index, c(1168:1169)]
39
-      cnSetExample <- sample.CNSet
38
+      cnSetExample <- hapmapSet[marker.index, c(1168:1169)]
40 39
       save(cnSetExample, file="~/Software/crlmm/data/cnSetExample.rda")
41 40
       ## 2 samples, many markers
42 41
 
... ...
@@ -44,11 +43,10 @@ data(cnSetExample2)
44 43
       snp.index <- which(isSnp(hapmapSet))[1:100]
45 44
       np.index <- which(!isSnp(hapmapSet))[1:100]
46 45
       marker.index <- c(snp.index, np.index)
47
-      sample.CNSet2 <- hapmapSet[marker.index, ]
48
-      sample.CNSet2$gender <- sample.CNSet2$gender[]
49
-      sample.CNSet2$SKW <- sample.CNSet2$SKW[]
50
-      sample.CNSet2$SNR <- sample.CNSet2$SNR[]
51
-      cnSetExample2 <- sample.CNSet2
46
+      cnSetExample2 <- hapmapSet[marker.index, ]
47
+      cnSetExample2$gender <- cnSetExample2$gender[]
48
+      cnSetExample2$SKW <- cnSetExample2$SKW[]
49
+      cnSetExample2$SNR <- cnSetExample2$SNR[]
52 50
       save(cnSetExample2, file="~/Software/crlmm/data/cnSetExample2.rda")
53 51
 }
54 52
 }
... ...
@@ -99,41 +99,7 @@ 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{
103
-	  ## we need the A and B intensities
104
-     if(require("ff")){
105
-	  path <- system.file("celFiles", package="hapmapsnp6")
106
-	  cels <- list.celfiles(path, full.names=TRUE)
107
-	  gtSet <- genotype(cels, batch=rep(1L,3))
108 102
 
109
-          gt1 <- calls(crlmmOutput)
110
-          gt2 <- calls(gtSet)[isSnp(gtSet), ]
111
-          gt2 <- gt2[match(rownames(gt2), rownames(gt1)), ]
112
-	  stopifnot(all.equal(gt1, gt2))
113
-
114
-	  XIndex <- which(chromosome(gtSet)==23 & isSnp(gtSet))
115
-	  YIndex <- which(chromosome(gtSet)==24 & isSnp(gtSet))
116
-	  A.X <- log2(A(gtSet)[XIndex,,drop=FALSE])
117
-	  B.X <- log2(B(gtSet)[XIndex,,drop=FALSE])
118
-	  meds.X <- apply(A.X+B.X, 2, median)/2
119
-	  A.Y <- log2(A(gtSet)[YIndex,,drop=FALSE])
120
-	  B.Y <- log2(B(gtSet)[YIndex,,drop=FALSE])
121
-	  meds.Y <- apply(A.Y+B.Y, 2, median)/2
122
-	  R <- meds.X - meds.Y
123
-	  SNR <- gtSet$SNR[]
124
-	  SNRMin <- 5
125
-	  gender <- kmeans(R, c(min(R[SNR[]>SNRMin]), max(R[SNR[]>SNRMin])))[["cluster"]]
126
-	  plot(seq_along(R), R, pch=21, bg=c("royalblue", "red")[gender],
127
-	       xlim=c(0, 4), xaxt='n', xlab="sample index")
128
-	  legend("topright", fill=c("royalblue", "red"), legend=c("male", "female"))
129
-	  delete(A(gtSet))
130
-	  delete(B(gtSet))
131
-	  delete(calls(gtSet))
132
-	  delete(assayDataElement(gtSet, "snpCallProbability"))
133
-	  ff.filenames <- list.files(".", pattern=".ff")
134
-	  unlink(ff.filename)
135
-     } ## end if(require("ff"))
136
-  }
137 103
 }
138 104
 
139 105
 \dontrun{
... ...
@@ -108,6 +108,45 @@ if (require(ff) & require(genomewidesnp6Crlmm) & require(hapmapsnp6)){
108 108
   gender[gender == "male"] <- 1
109 109
   dim(cnSet)
110 110
   table(isSnp(cnSet))
111
+  \dontrun{
112
+          ##---------------------------------------------------------------------------
113
+          ##
114
+	  ## Genotype Imputation
115
+          ##
116
+          ##---------------------------------------------------------------------------
117
+     if(require("ff")){
118
+	  path <- system.file("celFiles", package="hapmapsnp6")
119
+	  cels <- list.celfiles(path, full.names=TRUE)
120
+	  gtSet <- genotype(cels, batch=rep(1L,3))
121
+
122
+          gt1 <- calls(crlmmOutput)
123
+          gt2 <- calls(gtSet)[isSnp(gtSet), ]
124
+          gt2 <- gt2[match(rownames(gt2), rownames(gt1)), ]
125
+	  stopifnot(all.equal(gt1, gt2))
126
+
127
+	  XIndex <- which(chromosome(gtSet)==23 & isSnp(gtSet))
128
+	  YIndex <- which(chromosome(gtSet)==24 & isSnp(gtSet))
129
+	  A.X <- log2(A(gtSet)[XIndex,,drop=FALSE])
130
+	  B.X <- log2(B(gtSet)[XIndex,,drop=FALSE])
131
+	  meds.X <- apply(A.X+B.X, 2, median)/2
132
+	  A.Y <- log2(A(gtSet)[YIndex,,drop=FALSE])
133
+	  B.Y <- log2(B(gtSet)[YIndex,,drop=FALSE])
134
+	  meds.Y <- apply(A.Y+B.Y, 2, median)/2
135
+	  R <- meds.X - meds.Y
136
+	  SNR <- gtSet$SNR[]
137
+	  SNRMin <- 5
138
+	  gender <- kmeans(R, c(min(R[SNR[]>SNRMin]), max(R[SNR[]>SNRMin])))[["cluster"]]
139
+	  plot(seq_along(R), R, pch=21, bg=c("royalblue", "red")[gender],
140
+	       xlim=c(0, 4), xaxt='n', xlab="sample index")
141
+	  legend("topright", fill=c("royalblue", "red"), legend=c("male", "female"))
142
+	  delete(A(gtSet))
143
+	  delete(B(gtSet))
144
+	  delete(calls(gtSet))
145
+	  delete(assayDataElement(gtSet, "snpCallProbability"))
146
+	  ff.filenames <- list.files(".", pattern=".ff")
147
+	  unlink(ff.filename)
148
+     } ## end if(require("ff"))
149
+  }
111 150
 }
112 151
 }
113 152
 \keyword{ classif }