R/sharingProbabilityCalculations.R In RVS: Computes estimates of the probability of related individuals sharing a rare variant

Documented in denomProbexactSharingProbnumerProboneFounderSharingProbtwoFounderSharingProb

#' numerator of sharing probability
#' @keywords internal
#'
#' @description calculates the numerator of the sharing probability
#' outline in section 2.1 of Bureau et al.
#' @param gRain bayesian network
#' @param procPed pedigree object that has been process with processPedigree
#' @return numerator value
numerProb <- function(net, procPed)
{
rvInCarriers <- sapply(simplify=FALSE, FUN=function(dummy) 1:2,
X=as.character(procPed\$carriers))
noRvInNonCarriers <- sapply(simplify=FALSE, FUN=function(dummy) 0,
X=as.character(setdiff(procPed\$affected, procPed\$carriers)))
return(marginalProb(net, c(rvInCarriers, noRvInNonCarriers)))
}

#' denominator of sharing probability
#' @keywords internal
#'
#' @description calculates the denominator of the sharing probability
#' outline in section 2.1 of Bureau et al.
#' @param gRain bayesian network
#' @param procPed pedigree object that has been process with processPedigree
#' @return denominator value
denomProb <- function(net, procPed)
{
noRvInAny <- sapply(simplify=FALSE, FUN=function(dummy) 0,
X=as.character(procPed\$affected))
return(1 - marginalProb(net, noRvInAny))
}

#' calculate sharing probability in basic case
#' @keywords internal
#'
#' @description Assume that only one founder can introduce the variant to
#' the pedigree. Condition on each founder and sum over all resulting
#' probabilities.
#' @param procPed pedigree that has been through processPedigree()
#' @return sharing probability
oneFounderSharingProb <- function(procPed)
{
# set all founders to 0 (no variant)
net <- try(createNetwork(procPed))
if (class(net)=="try-error") stop("Creation of Bayesian network for the pedigree failed.")
net <- try(gRain::setEvidence(net, as.character(procPed\$founders),
rep('0', length(procPed\$founders))))
if (class(net)=="try-error") stop("Setting founder genotypes of the pedigree failed.")

# sum over probs, conditioning on each founder introducing variant
numer <- denom <- 0
for (f in procPed\$founders) #TODO: use sapply here
{
# condition on founder and calculate distribution
condNet <- gRain::retractEvidence(net, as.character(f))
condNet <- gRain::setEvidence(condNet, as.character(f), '1')

# compute probability
denom <- denom + denomProb(condNet, procPed)
numer <- numer + numerProb(condNet, procPed)
}
return(numer/denom)
}

#' sharing probability when founder pair introduces variant
#' @keywords internal
#'
#' @description In the case of relatedness among founders, assume that up
#' to two founders could introduce the variant and condition on all possible
#' pairs.
#' @param procPed pedigree that has been through processPedigree()
#' @param kinshipCoeff mean kinship coefficient among the founders
#' @param kinshipOrder order of the polynomial approximation to the distribtion
#' of the number of distinct alleles in the founders (d in Bureau et al.).
#' Must be <= 5
#' @return sharing probability
twoFounderSharingProb <- function(procPed, kinshipCoeff, kinshipOrder)
{
# set all founders to 0 (no variant)
net <- try(createNetwork(procPed))
if (class(net)=="try-error") stop("Creation of Bayesian network for the pedigree failed.")
net <- try(gRain::setEvidence(net, as.character(procPed\$founders),
rep('0', length(procPed\$founders))))
if (class(net)=="try-error") stop("Setting founder genotypes of the pedigree failed.")

# calculate kinship correction and initialize variables
numer <- denom <- 0
remainingFounders <- procPed\$founders
cor <- relatedFoundersCorrection(length(procPed\$founders),
kinshipCoeff, kinshipOrder)
for (f1 in procPed\$founders) #TODO: use sapply here
{
# condition on founder
net1 <- gRain::retractEvidence(net, as.character(f1))
net1 <- gRain::setEvidence(net1, as.character(f1), '1')

for (f2 in remainingFounders)
{
# condition on founder
net2 <- gRain::retractEvidence(net1, as.character(f2))
net2 <- gRain::setEvidence(net2, as.character(f2), '1')

# calculate (weighted) conditional probability
# Correction to formula 2 of Bioinformatics paper
w <- ifelse(f1==f2, cor, 2*(1-cor) / (length(procPed\$founders)-1))
numer <- numer + w * numerProb(net2, procPed)
denom <- denom + w * denomProb(net2, procPed)
}
remainingFounders <- setdiff(remainingFounders, f1)
}
return(numer/denom)
}

#' exact sharing probability calculation
#' @keywords internal
#'
#' @description Calculate the exact sharing probability given the minor allele
#' frequency among the founders (population).
#' @param procPed pedigree that has been through processPedigree()
#' @param alleleFreq allele frequency among the founders
#' @return sharing probability
exactSharingProb <- function(procPed, alleleFreq)
{
# create network with founders having a prior based on the allele freq
prior <- c((1-alleleFreq)^2, 2*alleleFreq*(1-alleleFreq), alleleFreq^2)
net <- createNetwork(procPed, prior)

# compute probability
return (numerProb(net, procPed) / denomProb(net, procPed))
}

Try the RVS package in your browser

Any scripts or data that you put into this service are public.

RVS documentation built on Nov. 8, 2020, 6:57 p.m.