LLR

LLR.R

pval.to.qual <- function(p) {
    return(-10*log10(p))
}

qual.to.pval <- function(q) {
    return(10^(-q/10))
}

LLR <- function(p, n, x) {
    ## p is editing ratio
    ## n is number of reads
    ## x is number of edited bases
    
    return(1/(n+1) * 1/beta(x+1, n-x+1) * p^x * (1-p)^(n-x))
}

LLR.binorm <- function(p, n, x) {
    return(choose(n, x) * p^n * (1-p)^(n-x))
}

GetObservedBases <- function(ref, alt) {
    ## Get observed bases with sequence error
    ref.bases <- 0
    alt.bases <- 0
    for (a in alt) {
        alt.bases <- alt.bases + a
    }
    for (ref in ref) {
        ref.bases <- ref.bases + (1-ref)
    }
    return(alt.bases + ref.bases)
}

## Observed data in the candidate editing site
##G.quals <- c(0.001)
##A.quals <- c(0.001, 0.001, 0.01, 0.01)
G.quals <- c(0.0001)
A.quals <- c(0.001, 0.001, 0.001, 0.001)
total.reads <- length(G.quals) + length(A.quals)

##G.obs <- (1-0.001) + 0.001 + 0.01 + 0.01
##A.obs <- 0.001 + (1-0.001) + (1-0.001) + (1-0.01) + (1-0.01)
G.obs <- GetObservedBases(G.quals, A.quals)
A.obs <- GetObservedBases(A.quals, G.quals)

## optimized editing ratio, r, error in obs count
r <- G.obs/total.reads

## errors in editing ratio
r2 <- (1-0.0001) / (1-0.001 + 1-0001 + 1-0.001 + 1-0.001 + 1-0.0001)

LLR.unedits <- LLR(sum(G.quals), total.reads, G.obs)
LLR.unedits2 <- LLR(sum(G.quals), total.reads, G.obs)

message("LLR 1: ",log10(r/LLR.unedits))
message("LLR 2: ",log10(r2/LLR.unedits2))