Single cell cytokines data, filter out genes with < 10 non-zero elements

Introduction

See if sva works when comparing two control groups in single cell cytokine study.

Focus on 8 types of cells, B_cells, CD4_T_cells, CD8_T_cells, NK_cells, Dendritic-cells, Ly6C+, Ly6C-,Neutrolphils.

Number of cells in each control group:

load('data/cytokine/DE_controls_output_filter10_addlimma.RData')
n.sample = c()
for(cell in names(output)){
  n.sample = rbind(n.sample,c(sum(output[[cell]]$group_idx),sum(1-output[[cell]]$group_idx)))
}
rownames(n.sample) = names(output)
colnames(n.sample) = c('Ctrl-1','Ctrl-2')
knitr::kable(n.sample,caption = 'Number of samples')

Number of genes with >= 10 non-zero elements in two control groups.

n.gene = c(14853)
for(cell in names(output)){
  n.gene = rbind(n.gene,c(14853-length(output[[cell]]$rm.idx)))
}
rownames(n.gene) = c('Total',names(output))
knitr::kable(n.gene,caption = 'Number of genes considered',col.names = '#genes')

Correlations between PCs and groups

Plot of correlations between groups and first 20 principle components.

par(mfrow=c(3,3))
for(cell in names(output)){
  plot(output[[cell]]$pc.cor[1:20],xlab='PCs',ylab='corr',main=paste(cell),ylim = c(-0.5,0.6),pch=20)
  abline(h=0,lty=3)
}

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

Plot of principle component that has maximum absolute correlation with groups for each cell. Vertical line separates two groups.

par(mfrow=c(3,3))
for(cell in names(output)){
  plot(output[[cell]]$pc.cor.max,xlab='',ylab='',
       main=paste(cell,', PC:',which.max(abs(output[[cell]]$pc.cor)),', corr:',round(max(abs(output[[cell]]$pc.cor)),2),sep=''),pch=1,col='grey50')
  abline(v=sum(output[[cell]]$group_idx),lty=3)
}

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

Compare p-value distributions from t-test, limma and sva-limma

Now compare distributions of p-values from t-test, limma and sva-limma:

par(mfrow=c(8,3))
for(cell in names(output)){
  hist(output[[cell]]$pvalue_t,main=paste(cell, ',t test'),xlab='')
  hist(output[[cell]]$pvalue_limma,main=paste(cell, ',limma'),xlab='')
  hist(output[[cell]]$pvalue_sva_limma,main=paste(cell, ',sva-limma'),xlab='')
}

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

Number of surrogate variables: based on the default method in sva - a permutation procedure originally prooposed by Buja and Eyuboglu 1992

n.sv = c()
for(cell in names(output)){
  n.sv = rbind(n.sv,output[[cell]]$sva_sva$n.sv)
}
rownames(n.sv) = names(output)
colnames(n.sv) = '#sv'
knitr::kable(n.sv,caption = 'Number of surrogate variables')

Compare the number of significant genes

The number of significant genes at \(fdr=0.05\) by BH procedure.

par(mfrow=c(1,1))
n.sig = c()
for(cell in names(output)){
  n.sig = rbind(n.sig,c(length(output[[cell]]$rej.idx.ttest),length(output[[cell]]$rej.idx.limma),length(output[[cell]]$rej.idx.sva)))
}
rownames(n.sig) = names(output)
colnames(n.sig) = c('t-test','limma','sva-limma')
knitr::kable(n.sig,caption = 'Number of significant genes at fdr=0.05')

The number of significant genes at \(fdr=0.01\) by BH procedure.

BH = function(p,alpha=0.05){
  n=length(p)
  i=rank(p)
  idx = which(p<=(i/n*alpha))
  if(length(idx)==0){
    NULL
  }else{
    i0= max(i[idx])
    rej.idx = which(i<=i0)
    rej.idx
  }
}

n.sig = c()
for(cell in names(output)){
  n.sig = rbind(n.sig,c(length(BH(output[[cell]]$pvalue_t,0.01)),length(BH(output[[cell]]$pvalue_limma,0.01)),
                        length(BH(output[[cell]]$pvalue_sva_limma,0.01))))
}
rownames(n.sig) = names(output)
colnames(n.sig) = c('t-test','limma','sva-limma')
knitr::kable(n.sig,caption = 'Number of significant genes at fdr=0.01')

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     

loaded via a namespace (and not attached):
 [1] workflowr_1.6.0 Rcpp_1.0.2      digest_0.6.18   later_0.7.5    
 [5] rprojroot_1.3-2 R6_2.3.0        backports_1.1.2 git2r_0.26.1   
 [9] magrittr_1.5    evaluate_0.12   highr_0.7       stringi_1.2.4  
[13] fs_1.3.1        promises_1.0.1  whisker_0.3-2   rmarkdown_1.10 
[17] tools_3.5.1     stringr_1.3.1   glue_1.3.0      httpuv_1.4.5   
[21] yaml_2.2.0      compiler_3.5.1  htmltools_0.3.6 knitr_1.20