Browse code

Add Lynn's calculateBAFandLRR function -- not exported

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

Rob Scharp authored on 01/10/2011 04:49:26
Showing2 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.31
4
+Version: 1.11.32
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>
... ...
@@ -2741,3 +2741,61 @@ ABpanel <- function(x, y, predictRegion,
2741 2741
 		}
2742 2742
 	}
2743 2743
 }
2744
+
2745
+calculateBAFandLRR <- function(cnSet) {
2746
+  ## Calculates B allele frequency and log2 R ratio for all snps for a given cnSet
2747
+  ## Returns a 3-D array of size Num Of Snps x Num of Samples x 2
2748
+  ## where result[,,1] is B allele frequency and result[,,2] is log2 R ratio
2749
+
2750
+  snp.index <- which(isSnp(cnSet))
2751
+  b <- (B(cnSet))[snp.index, ]
2752
+  a <- (A(cnSet))[snp.index, ]
2753
+  centers<-medians(cnSet, i=snp.index, j=1:length(unique(batch(cnSet))))
2754
+  bafAndLrr <- array(NA_integer_, dim=c(length(snp.index), length(sampleNames(cnSet)), 2), dimnames=list(featureNames(cnSet)[snp.index], sampleNames(cnSet),c("baf", "lrr")))
2755
+
2756
+  for (batch in unique(batch(cnSet))) {
2757
+
2758
+    ## get batch- and snp-specific centers
2759
+    center.aa.x <- centers[, 1, 1, batch]
2760
+    center.aa.y <- centers[, 2, 1, batch]
2761
+    center.ab.x <- centers[, 1, 2, batch]
2762
+    center.ab.y <- centers[, 2, 2, batch]
2763
+    center.bb.x <- centers[, 1, 3, batch]
2764
+    center.bb.y <- centers[, 2, 3, batch]
2765
+    theta.aa <- atan2(center.aa.y, center.aa.x) * 2 / pi
2766
+    r.aa <- center.aa.x + center.aa.y
2767
+    theta.ab <- atan2(center.ab.y, center.ab.x) * 2 / pi
2768
+    r.ab <- center.ab.x + center.ab.y
2769
+    theta.bb <- atan2(center.bb.y, center.bb.x) * 2 / pi
2770
+    r.bb <- center.bb.x + center.bb.y
2771
+
2772
+    index.sample <- which(batch(cnSet) %in% batch)
2773
+    for (i in index.sample) {
2774
+      theta <- atan2(b[, i], a[, i]) * 2 / pi
2775
+      r <- a[, i] + b[, i]
2776
+
2777
+      baf <- vector(mode="numeric", length=length(theta))
2778
+      baf[theta < theta.aa] <- 0
2779
+      baf[theta >= theta.bb] <- 1
2780
+      idx.ab <- which(theta >= theta.aa & theta < theta.ab)
2781
+      baf[idx.ab] <- .5 * (theta[idx.ab] - theta.aa[idx.ab]) / (theta.ab[idx.ab] - theta.aa[idx.ab])
2782
+      idx.bb <- which(theta >= theta.ab & theta < theta.bb)
2783
+      baf[idx.bb] <- .5 * (theta[idx.bb] - theta.ab[idx.bb]) / (theta.bb[idx.bb] - theta.ab[idx.bb]) + .5
2784
+      bafAndLrr[, i ,1] <- baf
2785
+
2786
+      r.expected <- vector(mode="numeric", length=length(r))
2787
+      r.expected[theta < theta.aa] <- r.aa[theta < theta.aa]
2788
+      r.expected[theta >= theta.bb] <- r.bb[theta >= theta.bb]
2789
+
2790
+      idx.ab <- which(theta >= theta.aa & theta < theta.ab)
2791
+      r.expected[idx.ab] <- (theta[idx.ab] - theta.aa[idx.ab]) * (r.ab[idx.ab]-r.aa[idx.ab]) /
2792
+        (theta.ab[idx.ab] - theta.aa[idx.ab]) + r.aa[idx.ab]
2793
+      idx.bb <- which(theta >= theta.ab & theta < theta.bb)
2794
+      r.expected[idx.bb] <- (theta[idx.bb] - theta.ab[idx.bb]) * (r.bb[idx.bb]-r.ab[idx.bb])/
2795
+        (theta.bb[idx.bb] - theta.ab[idx.bb]) + r.ab[idx.bb]
2796
+      bafAndLrr[, i, 2] <- log2(r / r.expected)
2797
+
2798
+    }
2799
+  }
2800
+  return(bafAndLrr)
2801
+}