binomial thinning

Introduction

Investigate why effects with larger se are bigger.

Assume we have \(n\) samples and a fraction \(p\) of them belong to group 1 and the rest belong to group 2. So \(x = (1,1,1,...,1,0,0,0,...,0)^T\in R^n\) and \(\sum_ix_i=np\). Under this setting, in simple linear regression \(y = a+\beta x + \epsilon\), \(\epsilon\sim N(0,\sigma^2)\), the variance of \(\hat\beta\) is \(\hat s^2 = \frac{n\sigma^2}{n\sum_ix_i^2-(\sum_ix_i)^2}=\frac{\sigma^2}{np-np^2}\). For fixed \(n\) and \(p\), if \(\hat s\) is large, then this means \(\hat\sigma^2\) is large hence \(\sigma^2\) is large.

We now need to figure out the relationship between \(\beta\) and \(\sigma^2\).

Let’s assume we have RNA-Seq count data \(z_i\sim Poisson(\lambda)\) for \(i=1,2,...,n\). In binomial thinning, \(\beta\) is the log2 fold change between groups. Now assume \(\beta>0\), according to Gerard and Stephens(2017), the new(thinned) data vector is \(w_i\sim Poisson(\mu_i)\), where \(\mu_i=2^{-\beta(1-x_i)}\lambda\). The response \(y\) in the simple linear regression is the log transformation of \(w\), \(y_i=\log(w_i)\), \(i=1,2,...n\).

The Taylor series expansion of \(\log w_i\) around \(\mu_i\) is \(\log(w_i)\approx \log(\mu_i)+\frac{w_i-\mu_i}{\mu_i}\). So the mean of \(\log(w_i)\) is \(\log(\mu_i) = \lambda - \beta(1-x_i)\) and variance \(\frac{1}{\mu_i} = \frac{1}{2^{-\beta(1-x_i)}\lambda}\). So if \(\beta\) is large, then \(Var(\log(w_i))\) is large if \(x_i=0\). This explains the why effects with larger se are bigger.

Check

For non-null genes \(j\in Non.null.gene.set\), choose \(p_{1j}+p_{2j}=1\) and \(\frac{p_{1j}}{p_{2j}}=\exp(\beta_j)\), and thin the counts \(w_{ij}\sim Binomial(z_{ij},p_{\{group.of.i\}j})\), where \(z_{ij}\) is the observed counts for \(i\)th sample, and group of \(i\) is either 1 or 2.

For null genes \(j\in Null.gene.set\), \(w_{ij}\sim Binomial(z_{ij},0.5)\)

#'@param Z count matrix, sample by features
#'@param x 1 for group 1, 0 for group 2
#'@param beta effect of fearures,  0 for null.
#'@return W, thinned matrix
bi_thin = function(Z,x,beta){

  n=nrow(Z)
  p=ncol(Z)
  
  # group index
  g1 = which(x==1)
  g2 = which(x==0)
  
  
  p2 = 1/(1+exp(beta))
  p1 = 1-p2
  P = matrix(nrow = n,ncol = p)
  P[g1,] = t(replicate(length(g1),p1))
  P[g2,] = t(replicate(length(g2),p2))
  
  W = matrix(rbinom(n*p,Z,P),nrow=n)
  
  W
}

quiet <- function(x) { 
  sink(tempfile()) 
  on.exit(sink()) 
  invisible(force(x)) 
} 

Normal signal, sd=1.5. Run 50 reps.

library(sva)
load('data/scde/scCDT.RData')
Z = t(as.matrix(CDT))
rm.idx = which(colSums(Z!=0)<30)
Z = Z[,-rm.idx]
n = nrow(Z)
p = ncol(Z)

set.seed(12345)

nreps = 30

loglik=c()
roc_result = c()

for(rep in 1:nreps){
  
x = rbinom(n,1,0.5)
beta = rnorm(p,0,1.5)
beta[sample(1:p,p*0.9)]=0

W = bi_thin(Z,x,beta)

Wn = log(W+0.5)

X = model.matrix(~x)

sva_sva = quiet(sva(t(Wn),mod = X, mod0 = X[, -2, drop = FALSE], n.sv = 3))

X.sva = cbind(X, sva_sva$sv)
lmout = limma::lmFit(object = t(Wn), design = X.sva)
eout  = limma::eBayes(lmout)

svaout           <- list()
svaout$betahat   <- lmout$coefficients[, 2]
svaout$sebetahat <- lmout$stdev.unscaled[, 2] * sqrt(eout$s2.post)
svaout$pvalues   <- eout$p.value[, 2]

sva_limma_ash0 = ashr::ash(svaout$betahat,svaout$sebetahat,alpha=0)

sva_limma_ash1 = ashr::ash(svaout$betahat,svaout$sebetahat,alpha=1)

loglik = rbind(loglik,c(sva_limma_ash0$loglik,sva_limma_ash1$loglik))

#knitr::kable(cbind(sva_limma_ash0$loglik,sva_limma_ash1$loglik), 
#             col.names = c('alpha=0','alpha=1'), digits = 2,caption = 'log-lik')


which_null = ifelse(beta==0,1,0)
################
roc_out <- list(
  #pROC::roc(response = which_null, predictor = c(mout$result$lfsr)),
  #pROC::roc(response = which_null, predictor = c(mout1$result$lfsr)),
  pROC::roc(response = which_null, predictor = c(svaout$pvalues)),
  pROC::roc(response = which_null, predictor = c(sva_limma_ash0$result$lfsr)),
  pROC::roc(response = which_null, predictor = c(sva_limma_ash1$result$lfsr)))
#name_vec <- c("MOUTHWASH0","MOUTHWASH1","SVA-limma","SVA-limma-ash0","SVA-limma-ash1")

auc_vec <- sapply(roc_out, FUN = function(x) { x$auc })
roc_result = rbind(roc_result,auc_vec)

#knitr::kable(sort(auc_vec, decreasing = TRUE), col.names = "AUC", digits = 3)

}

name_vec <- c("sva-limma","sva-limma-ash0","sva-limma-ash1")
colnames(roc_result) <- name_vec

colnames(loglik) = c('alpha=0','alpha=1')

boxplot(loglik,ylab = 'loglik')

Past versions of unnamed-chunk-2-1.png

boxplot(roc_result,ylab='AUC')

Past versions of unnamed-chunk-2-2.png

library(sva)
load('data/scde/scCDT.RData')
Z = t(as.matrix(CDT))
rm.idx = which(colSums(Z!=0)<30)
Z = Z[,-rm.idx]
n = nrow(Z)
p = ncol(Z)

set.seed(12345)

nreps = 30

loglik=c()
roc_result = c()

beta = rnorm(p,0,1.5)
beta[sample(1:p,p*0.9)]=0

for(rep in 1:nreps){
  
x = rbinom(n,1,0.5)


W = bi_thin(Z,x,beta)

Wn = log(W+0.5)

X = model.matrix(~x)

sva_sva = quiet(sva(t(Wn),mod = X, mod0 = X[, -2, drop = FALSE], n.sv = 3))

X.sva = cbind(X, sva_sva$sv)
lmout = limma::lmFit(object = t(Wn), design = X.sva)
eout  = limma::eBayes(lmout)

svaout           <- list()
svaout$betahat   <- lmout$coefficients[, 2]
svaout$sebetahat <- lmout$stdev.unscaled[, 2] * sqrt(eout$s2.post)
svaout$pvalues   <- eout$p.value[, 2]

sva_limma_ash0 = ashr::ash(svaout$betahat,svaout$sebetahat,alpha=0)

sva_limma_ash1 = ashr::ash(svaout$betahat,svaout$sebetahat,alpha=1)

loglik = rbind(loglik,c(sva_limma_ash0$loglik,sva_limma_ash1$loglik))

#knitr::kable(cbind(sva_limma_ash0$loglik,sva_limma_ash1$loglik), 
#             col.names = c('alpha=0','alpha=1'), digits = 2,caption = 'log-lik')


which_null = ifelse(beta==0,1,0)
################
roc_out <- list(
  #pROC::roc(response = which_null, predictor = c(mout$result$lfsr)),
  #pROC::roc(response = which_null, predictor = c(mout1$result$lfsr)),
  pROC::roc(response = which_null, predictor = c(svaout$pvalues)),
  pROC::roc(response = which_null, predictor = c(sva_limma_ash0$result$lfsr)),
  pROC::roc(response = which_null, predictor = c(sva_limma_ash1$result$lfsr)))
#name_vec <- c("MOUTHWASH0","MOUTHWASH1","SVA-limma","SVA-limma-ash0","SVA-limma-ash1")

auc_vec <- sapply(roc_out, FUN = function(x) { x$auc })
roc_result = rbind(roc_result,auc_vec)

#knitr::kable(sort(auc_vec, decreasing = TRUE), col.names = "AUC", digits = 3)

}

name_vec <- c("sva-limma","sva-limma-ash0","sva-limma-ash1")
colnames(roc_result) <- name_vec

colnames(loglik) = c('alpha=0','alpha=1')

boxplot(loglik,ylab = 'loglik')

boxplot(roc_result,ylab='AUC')

Session information

sessionInfo()

R version 3.5.1 (2018-07-02)
Platform: x86_64-pc-linux-gnu (64-bit)
Running under: Scientific Linux 7.4 (Nitrogen)

Matrix products: default
BLAS/LAPACK: /software/openblas-0.2.19-el7-x86_64/lib/libopenblas_haswellp-r0.2.19.so

locale:
 [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C              
 [3] LC_TIME=en_US.UTF-8        LC_COLLATE=en_US.UTF-8    
 [5] LC_MONETARY=en_US.UTF-8    LC_MESSAGES=en_US.UTF-8   
 [7] LC_PAPER=en_US.UTF-8       LC_NAME=C                 
 [9] LC_ADDRESS=C               LC_TELEPHONE=C            
[11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C       

attached base packages:
[1] stats     graphics  grDevices utils     datasets  methods   base     

other attached packages:
[1] sva_3.30.0          BiocParallel_1.16.0 genefilter_1.64.0  
[4] mgcv_1.8-25         nlme_3.1-137       

loaded via a namespace (and not attached):
 [1] Biobase_2.42.0       bit64_0.9-7          splines_3.5.1       
 [4] foreach_1.4.4        assertthat_0.2.0     mixsqp_0.2-2        
 [7] stats4_3.5.1         blob_1.1.1           yaml_2.2.0          
[10] pillar_1.3.1         RSQLite_2.1.1        backports_1.1.2     
[13] lattice_0.20-38      glue_1.3.0           limma_3.38.2        
[16] pROC_1.13.0          digest_0.6.18        promises_1.0.1      
[19] colorspace_1.3-2     htmltools_0.3.6      httpuv_1.4.5        
[22] Matrix_1.2-15        plyr_1.8.4           XML_3.98-1.16       
[25] pkgconfig_2.0.2      purrr_0.3.2          xtable_1.8-3        
[28] scales_1.0.0         whisker_0.3-2        later_0.7.5         
[31] git2r_0.26.1         tibble_2.1.1         annotate_1.60.0     
[34] IRanges_2.16.0       ggplot2_3.1.1        ashr_2.2-39         
[37] BiocGenerics_0.28.0  lazyeval_0.2.1       survival_2.43-1     
[40] magrittr_1.5         crayon_1.3.4         memoise_1.1.0       
[43] evaluate_0.12        fs_1.3.1             doParallel_1.0.14   
[46] MASS_7.3-51.1        truncnorm_1.0-8      tools_3.5.1         
[49] matrixStats_0.54.0   stringr_1.3.1        S4Vectors_0.20.1    
[52] munsell_0.5.0        AnnotationDbi_1.44.0 compiler_3.5.1      
[55] rlang_0.4.0          grid_3.5.1           RCurl_1.95-4.11     
[58] iterators_1.0.10     bitops_1.0-6         rmarkdown_1.10      
[61] gtable_0.2.0         codetools_0.2-15     DBI_1.0.0           
[64] R6_2.3.0             knitr_1.20           dplyr_0.8.0.1       
[67] bit_1.1-14           workflowr_1.6.0      rprojroot_1.3-2     
[70] stringi_1.2.4        pscl_1.5.2           parallel_3.5.1      
[73] SQUAREM_2017.10-1    Rcpp_1.0.2           tidyselect_0.2.5