zzh
/
gatk-3.8
1
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

Calculate the phasing probabilities for particular intra-het distances 

git-svn-id: file:///humgen/gsa-scr1/gsa-engineering/svn_contents/trunk@4838 348d0f76-0448-11de-a6fe-93d51630548a
This commit is contained in:
fromer 2010-12-14 18:44:59 +00:00
parent 526ae92093
commit c167c6f9eb
2 changed files with 48 additions and 0 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

1
R/phasing/RBP_theoretical.R
View File

@ -122,6 +122,7 @@ pPhaseHetPairAtDistanceUsingDepthAndWindow <- function(windowDistances, phaseInd
if (phaseIndex < 2 || phaseIndex > n) { if (phaseIndex < 2 || phaseIndex > n) {
stop("phaseIndex < 2 || phaseIndex > n") stop("phaseIndex < 2 || phaseIndex > n")
} }
#print(paste("windowDistances= (", paste(windowDistances, collapse=", "), ")", sep=""))
# A. Pre-compute the upper diagonal of square matrix of n CHOOSE 2 values of: # A. Pre-compute the upper diagonal of square matrix of n CHOOSE 2 values of:
# pDirectlyPhaseHetPairAtDistanceUsingDepth(meanDepth, nReadsToPhase, L, dist(i,j), Im, Is) # pDirectlyPhaseHetPairAtDistanceUsingDepth(meanDepth, nReadsToPhase, L, dist(i,j), Im, Is)

47
R/phasing/calcPhasingProbsForWindowDistances.R 100644
View File

@ -0,0 +1,47 @@
calcPhasingProbsForWindowDistances <- function(distances, MAX_WINDOW_SIZE, FILE_NAME = NULL) {
WINDOW_SIZES = 2:MAX_WINDOW_SIZE
phaseProbsPositionWindow = matrix(data = NA, nrow=length(distances), ncol=length(WINDOW_SIZES))
for (i in 1:length(distances)) {
# Try to phase (i+1)-st position [relative to i] using varying window sizes:
for (j in 1:length(WINDOW_SIZES)) {
windowSize = WINDOW_SIZES[j]
remainingSize = windowSize - 2 # exlcude i, i+1
numOnLeft = i - 1
numOnRight = (length(distances) + 1) - (i + 2) + 1
if (numOnLeft <= numOnRight) {
halfToUse = floor(remainingSize / 2) # skimp on the left [floor], and be generous with the right side
useOnLeft = min(halfToUse, numOnLeft)
useOnRight = min(remainingSize - useOnLeft, numOnRight)
}
else {
halfToUse = ceiling(remainingSize / 2) # be generous with the right side [ceiling]
useOnRight = min(halfToUse, numOnRight)
useOnLeft = min(remainingSize - useOnRight, numOnLeft)
}
startInd = i - useOnLeft # go left from position i
stopInd = i + 1 + useOnRight # go right from position i + 1
usePositionRange = seq(from=startInd, to=stopInd, by=1)
useDistancesRange = seq(from=startInd, to=stopInd-1, by=1) # since there are N-1 distances between N consecutive positions
phaseIndex = which(usePositionRange == i+1)
if (length(phaseIndex) != 1) stop("NO phaseIndex!")
windowDistances = distances[useDistancesRange]
print(paste("Try to phase position ", i+1, " [relative to ", i, "] using positions: (", paste(usePositionRange, collapse=", "), "), windowDistances= (", paste(windowDistances, collapse=", "), "), [phaseIndex= ", phaseIndex, ", i=", i, "]", sep=""))
p = pPhaseHetPairAtDistanceUsingDepthAndWindow(windowDistances, phaseIndex, meanDepth, nReadsToPhase, L, Im, Is)
print(paste("phase prob: ", p, sep=""))
phaseProbsPositionWindow[i, j] = p
}
if (!is.null(FILE_NAME)) {
save(list = ls(all=TRUE), file = paste(FILE_NAME, ".RData", sep=""))
}
}
phaseProbsPositionWindow
}