From 66931d433cd5e570133bd5d286cc5e2c98eafdce Mon Sep 17 00:00:00 2001
From: depristo <depristo@348d0f76-0448-11de-a6fe-93d51630548a>
Date: Wed, 30 Jun 2010 16:38:49 +0000
Subject: [PATCH] useful routines for R

git-svn-id: file:///humgen/gsa-scr1/gsa-engineering/svn_contents/trunk@3685 348d0f76-0448-11de-a6fe-93d51630548a
---
 R/titvFPEst.R | 43 +++++++++++++++++++++++++++++++++++++++++++
 1 file changed, 43 insertions(+)
 create mode 100755 R/titvFPEst.R

diff --git a/R/titvFPEst.R b/R/titvFPEst.R
new file mode 100755
index 000000000..19d6d2f65
--- /dev/null
+++ b/R/titvFPEst.R
@@ -0,0 +1,43 @@
+titvFPEst <- function(titvExpected, titvObserved) { 1 - (titvObserved - 0.5) / (titvExpected - 0.5) }
+
+calcHet <- function(nknown, knownTiTv, nnovel, novelTiTv, callable) {
+  TP <- nknown + (1-titvFPEst(knownTiTv, novelTiTv)) * nnovel
+  2 * TP / 3 / callable
+}
+
+marginalTiTv <- function( nx, titvx, ny, titvy ) {
+    tvx = nx / (titvx + 1)
+    tix = nx - tvx
+    tvy = ny / (titvy + 1)
+    tiy = ny - tvy
+    tiz = tix - tiy
+    tvz = tvx - tvy
+    return(tiz / tvz)
+}
+marginaldbSNPRate <- function( nx, dbx, ny, dby ) {
+    knownx = nx * dbx / 100
+    novelx = nx - knownx
+    knowny = ny * dby / 100
+    novely = ny - knowny
+    knownz = knownx - knowny
+    novelz = novelx - novely
+    return(knownz / ( knownz + novelz ) * 100)
+}
+
+numExpectedCalls <- function(L, theta, calledFractionOfRegion, nIndividuals, dbSNPRate) {
+    nCalls <- L * theta * calledFractionOfRegion * sum(1 / seq(1, 2 * nIndividuals))
+    return(list(nCalls = nCalls, nKnown = dbSNPRate * nCalls, nNovel = (1-dbSNPRate) * nCalls))
+}
+
+normalize <- function(x) {
+    x / sum(x)
+}
+
+normcumsum <- function(x) {
+    cumsum(normalize(x))
+}
+
+cumhist <- function(d) {
+    h <- hist(d)
+    #plot(h$mids, cumsum(h$count), type="b", col="orange", lwd=2)
+}