Single cell cytokines data

Introduction

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

Take a look at the annotation file:

data_ann = read.csv('/project2/mstephens/chevrier-stephens/data/whole_cyto_annot.csv.gz')
data_ann[1:5,]

                  X0 sample well n_counts log_counts n_genes     mt_frac
1 AAACCTGAGAATAGGG-1  CCL20    1     1611   7.384610     852 0.030415890
2 AAACCTGAGAGTAAGG-1  CXCL1    1     6482   8.776784    2340 0.032860227
3 AAACCTGAGTGGCACA-1  CCL20    1     1940   7.570443    1058 0.048453607
4 AAACCTGCAAGCTGGA-1  CXCL1    1     2186   7.689829    1014 0.030192131
5 AAACCTGCAAGGACAC-1  CCL20    1      849   6.744059     469 0.001177856
  louvain louvain_sub_5 louvain_sub_9 subclustered   cell_type
1       0             0             0            0     B_cells
2       0             0             0            0     B_cells
3       1             1             1            1     B_cells
4       0             0             0            0     B_cells
5       3             3             3            3 Neutrophils

sample: treatment or cytokines, including two controls.

well? There are 25 wells, each of which contains two samples.

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

#read data matrix
# library(hdf5r)
# library(Matrix)
# f <- H5File$new("whole_cyto_normalized.h5ad",mode = "r")
# print(names(f))
# out <- f[["X"]]
# print(h5attributes(out))
# i <- out[["indices"]][]
# j <- out[["indptr"]][]
# x <- out[["data"]][]
#library(SparseM)
# X.csr = new('matrix.csr',ra=x,ja=as.integer(i+1),ia=as.integer(j+1),dimension=h5attributes(out)$h5sparse_shape)
#load("data/cytokine_normalized.RData")
load('data/cytokine/DE_controls_output.RData')

Number of cells in each control group:

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 that have at least one measurement 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-5-1.png

Plot of principle component that has maximum absolute correlation with groups for each cell. Vertical line seperates 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-6-1.png

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

Note: We will see a lot of p-values from t-test around 0.3-0.4. In single cell DE study, some genes are only measured once in one group while have no observation in another group. For example, gene expression in group 1 \(= (0,0,...,0,0)\) and gene expression in group 2 \(= (x,0,...,0,0)\). So in this case, unequal variance two-sample t-test always gives t-statistic = \(1\) with \(df=n_2-1\), where \(n_2\) is the number of samples in group 2. Let’s plot p-value vs df.

Suppose we have at least 5 samples in group 2, then p-value starts at 0.3739(df=4) and converges to 0.3173 as df goes to infinite.

par(mfrow=c(1,1))
plot(4:1e3,(1-pt(1,4:1e3))*2,xlab='df',ylab='p-value',main='t-statistics = 1',pch=20)

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

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

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

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

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.sva)))
}
rownames(n.sig) = names(output)
colnames(n.sig) = c('t-test','sva')
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_sva_limma,0.01))))
}
rownames(n.sig) = names(output)
colnames(n.sig) = c('t-test','sva')
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