R/quantification-model.R

In djhshih/orient-bias-filter: Orientation Bias Filter

#' Estimate global phi using a consistent estimator.
#'
#' For oxoG damage, we restrict the analysis to G>T and C>A variant sites as
#' well as G and C non-variant sites.  We consider G>T on read 1 and C>A on
#' read 2 as possibly oxoG artifact or genuine variant, and we consider all
#' oxoG-consistent variants as genuine variants:
#' C>A on read 1 and G>T on read 2.
#'
#' Note that we must examine observed original nucleotide of the
#' read, so reads mapping to the minus reference strand should be
#' reverse-complemented back to the original sequence.
#' (Reads mapping to the minus reference strand has been
#' reverse-complemeneted by the aligner.)
#' 
#' `xc = xr + xd`
#' where xr is the unobserved count of damage-consistent but real variants, 
#' and xd is the unobserved count of damage-induced artificial variants.
#' 
#' @param xc   total alt count of reads that are consistent with damage
#'             (e.g. for oxoG, G>T on read 1; C>A on read 2)
#' @param xi   total alt count of reads that are inconsistent with damage
#'             (e.g. for oxoG, C>A on read 2; G>T on read 2)
#' @param nc   total count of reads that are consistent
#'             (e.g. for oxoG, G>G and G>T on read 1; C>C and C>A on read 2)
#' @param ni   total count of reads that are inconsistent
#'             (e.g. for oxoG, C>C and C>A on read 2; G>G and G>T on read 2)
#' @return estimated value of global damage (\code{phi}) under a model in which 
#'         every locus has the same damage probability
#' @export
quantify_global_damage <- function(xc, xi, nc, ni) {
	theta_hat <- quant_estimate_global_theta(xi, ni);
	quant_estimate_global_phi(xc, nc, theta_hat)
}

# Estimate global theta using a consistent estimator.
# 
# @param xi   total alt count of reads that are inconsistent with damage
#             (e.g. for oxoG, C>A on read 2; G>T on read 2)
# @param ni   total count of reads that are inconsistent
#             (e.g. for oxoG, C>C and C>A on read 2; G>G and G>T on read 2)
# @return estimated value of theta
# @export
quant_estimate_global_theta <- function(xi, ni) {
	xi / ni
}

quant_estimate_global_phi <- function(xc, nc, theta_hat) {
	phi_hat <- (xc/nc - theta_hat) / (1 - theta_hat);
	phi_hat[phi_hat < 0] <- 0;
	phi_hat
}