R/calc_ex.R

In ggruenhagen3/coselens: COnditional SELection on the Excess of NonSynonymous Substitutions

# Calculate the excess of missense and nonsense mutations per patient. Input dnds1 is a dndsout object generated by running dndscv on a list of known cancer genes.
calc_ex = function(dnds1) {
  
  mle_tcv = function(n_neutral, exp_rel_neutral, shape, scale) {
    tml = (n_neutral+shape-1)/(exp_rel_neutral+(1/scale))
    if (shape<=1) { # i.e. when theta<=1
      tml = max(shape*scale,tml) # i.e. tml is bounded to the mean of the gamma (i.e. y[9]) when theta<=1, since otherwise it takes meaningless values
    }
    return(tml)
  }
  
  calc_ex_single = function(j) {
    
    indneut = 1
    num_patients <- length(as.vector(unique(dnds1$annotmuts$sampleID)))
    
    y <- as.numeric(dnds1$genemuts[j,-1])
    exp_rel <- y[5:8]/y[5]
    shape <- dnds1$nbreg$theta
    scale <- y[9]/dnds1$nbreg$theta
    opt_t <- mle_tcv(n_neutral=sum(y[indneut]), exp_rel_neutral=sum(exp_rel[indneut]), shape=shape, scale=scale)
    mrfold = max(1e-10, opt_t/sum(y[5]))
    excess_mis <- (y[2]        - (y[6]       * mrfold)) / num_patients
    excess_non <- (sum(y[3:4]) - (sum(y[7:8])* mrfold)) / num_patients
    
    return(c(excess_mis,excess_non))
  }
  
  ex_mut = as.data.frame(t(sapply(1:nrow(dnds1$genemuts), calc_ex_single)))
  ex_mut = cbind(dnds1$genemuts[,1],ex_mut)
  colnames(ex_mut) = c("gene_name", "ex_mis", "ex_non")
  
  return(ex_mut)
}