#' prior functions
# P(m|M,C) Given a context (including reference), and the fact that a mutation has
# occurred what is the probability that the mutation is to each of the three possible alternate bases.
# 96 rows, each (context,cref) combination (n=32) has three elements, and the sum of their probabilities
# is 1
.compute_priors <- function(vr,reference,prior_lod){
prior_vars <- vr %>% dplyr::filter(mutect_odds >= prior_lod & pass_all)
mutation_prior <- .compute_mutation_prior(prior_vars)
context_prior <- .compute_context_prior(prior_vars)
global_prior <- .compute_global_prior(reference)
#P(M) = 3e-6 = mu and P(C) = 1/96 (there are 96 context) P(M)/P(C) = 96*3e-6 = 2.88e-4
# full prior= mutation_prior * context_prior * P(M)/P(C)
# 4/23/2018: Add the real global context prior. replacing P(C) = 1/96
joint_prior <- mutation_prior %>%
dplyr::left_join(context_prior,by=c('context','cref')) %>%
dplyr::left_join(global_prior,by=c('context','cref')) %>%
dplyr::mutate(joint_prior=((mutation_prior*context_prior)/global_prior)*3e-6) %>%
dplyr::select(cref,calt,context,mutation_prior,joint_prior,global_prior)
vr <- vr %>% dplyr::left_join(joint_prior,by=c('context','cref','calt')) %>%
tidyr::replace_na(list(joint_prior=1e-6)) %>%
dplyr::mutate(logit_prior=log10(joint_prior/(1-joint_prior)),
prior_odds=10**(TLOD+logit_prior))
vr
}
.compute_mutation_prior <- function(vars) {
vars %>% dplyr::group_by(context,cref,calt) %>%
dplyr::summarise(n=n()) %>%
dplyr::mutate(mutation_prior=n/sum(n)) %>%
dplyr::select(-n)
}
.compute_context_prior <- function(vars){
vars %>% dplyr::group_by(context,cref) %>%
dplyr::summarise(n=n()) %>% dplyr::ungroup() %>%
dplyr::mutate(context_prior=n/sum(n)) %>%
dplyr::select(-n)
}
.compute_global_prior <- function(reference){
seqs <- SomaticSignatures::kmerFrequency(reference,n=1e6,k=3)
fr <- data.frame(ctxt=as.character(names(seqs)),global_prior=as.numeric(seqs))
fr <- fr %>% dplyr::mutate(context=paste0(stringr::str_sub(ctxt,1L,1L),'.',stringr::str_sub(ctxt,-1L,-1L)),
cref=stringr::str_sub(ctxt,2L,2L))
fr
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.