Permutation test – trans-gene analysis (t-test + ash)

library(ashr)
library(ggplot2)
library(knitr)
library(kableExtra)
library(gridExtra)
options(stringsAsFactors = F)
wkdir <- '~/Downloads/ASoC/singlecell/'
load(paste0(wkdir,'data/cropseq_expression.Rd'))
nlocus <- colSums(exp.per.enhancer>0)
ncell_uniq <- rep(NA,nrow(exp.per.gRNA))
for (i in 1:nrow(exp.per.gRNA)){
  glocus <- row.names(exp.per.gRNA)[i]
  ncell_uniq[i]=sum(exp.per.gRNA[glocus,]>0 & nlocus==1)
}
tmp <- sapply(strsplit(row.names(exp.per.gRNA),split = '_'), 
              function(x){paste(x[1],x[2],sep = '_')})
names(ncell_uniq) <- tmp

DE analysis on permutated sctransformed data

We still use the sctransformed Pearson residuals, and we filtered the genes to be those present in at least 20% of all cells (8117 in total).

For each target gRNA, we selected 2 groups of cells: cells that only contain the target gRNAs, and cells that only contain the negative control gRNAs. Then, we permute the group labels and conducted categorical regression of the Pearson residuals (against the 2 permuted conditions) to obtain the summary statistics. As before, 70 (gRNA vs neg ctrl) pairs were tested.

ashr meta-analysis on regression summary statistics

For each gene, we combined the 70 DE effect sizes and standard errors estimated in the previous step, and used an adaptive shrinkage method, ashr, to discern which effect sizes are truly non-zero (significant). Here, we used a cutoff of 0.1 for lfsr (the local false sign rate).

get.result <- function(betaest.mtx,lfsr.mtx){
  
  tmp.gRNA <- colSums(lfsr.mtx<0.1)
  tmp.gene <- rowSums(lfsr.mtx<0.1)
  signif_genes <- row.names(betaest.mtx)[tmp.gene>0]
  cat('There are',length(signif_genes),'significant genes in total.')
  
  signif_genes.gRNA.lst <- list()
  for(i in signif_genes){
    gRNA.group <- colnames(betaest.mtx)[lfsr.mtx[i,]<0.1]
    signif_genes.gRNA.lst[[i]] <- sapply(strsplit(gRNA.group,split = '_'), 
                                         function(x){paste(x[1],x[2],sep = '_')})
  }
  # Make table:
  compact.genes.lst <- sapply(signif_genes.gRNA.lst,FUN = function(x){paste(x,collapse = ', ')})
  compact.genes.df <- as.data.frame(compact.genes.lst)
  tmp.indx <- order(rownames(compact.genes.df))
  compact.genes.df <- data.frame(signif_gene=rownames(compact.genes.df)[tmp.indx],
                                 gRNA_targets=compact.genes.df[tmp.indx,])
  # Present table:
  print(kable(compact.genes.df, caption = 'gRNA targets per significant gene') %>% 
          kable_styling() %>% scroll_box(width = "100%", height = "300px"))
  num_gRNA_pergene <- sapply(signif_genes.gRNA.lst,FUN = function(x){length(x)})

  ### Significant genes per gRNA target
  gRNA.signif_genes.lst <- list()
  for (i in signif_genes){
    cells <- signif_genes.gRNA.lst[[i]]
    for (c in cells){
      if (is.na(match(c, names(gRNA.signif_genes.lst)))){
        gRNA.signif_genes.lst[[c]] <- i
      } else {
        gRNA.signif_genes.lst[[c]] <- c(gRNA.signif_genes.lst[[c]],i)
      }
    }
  }
  # Make table:
  compact.gRNAs.lst <- sapply(gRNA.signif_genes.lst,FUN = function(x){paste(x,collapse = ', ')})
  compact.gRNAs.df <- as.data.frame(compact.gRNAs.lst)
  tmp.indx <- order(rownames(compact.gRNAs.df))
  compact.gRNAs.df <- data.frame(gRNA_target=rownames(compact.gRNAs.df)[tmp.indx],
                                 signif_genes=compact.gRNAs.df[tmp.indx,])
  compact.gRNAs.df$gRNA_target <- sapply(strsplit(compact.gRNAs.df$gRNA_target,split = '_'),
                                         function(x){paste(x[1],x[2],sep = '_')})
  compact.gRNAs.df$num_cells <- ncell_uniq[compact.gRNAs.df$gRNA_target]
  compact.gRNAs.df <- compact.gRNAs.df[,c(3,1,2)]
  # Present table:
  print(kable(compact.gRNAs.df, caption = 'Significant genes per gRNA target') %>% 
          kable_styling() %>% scroll_box(width = "100%", height = "300px"))
  num_gene_pergRNA <- sapply(gRNA.signif_genes.lst,FUN = function(x){length(x)})
  
  # Plots:
  plot1 <- qplot(num_gRNA_pergene,geom = 'histogram',xlab = 'Number of gRNA targets per gene',ylab = 'Count',
        main = paste('Distribution of', length(num_gRNA_pergene), 'significant genes'),bins=10)
  plot2 <- qplot(num_gene_pergRNA,geom = 'histogram',xlab = 'Number of DE genes per gRNA target',
        ylab = 'Count', main = paste('Distribution of', length(num_gene_pergRNA), 'gRNAs'),bins=20)
  plot3 <- qplot(compact.gRNAs.df$num_cells,num_gene_pergRNA[compact.gRNAs.df$gRNA_target],
      xlab = 'Number of cells in the target condition',ylab = 'Number of significant genes')
  grid.arrange(plot1,plot2,plot3,ncol=2)
}

Permutation 1

Here are all the significant genes (\(lfsr<0.1\)) and the gRNA groups where they are differentially expressed:

load(paste0(wkdir,'data/perm.ashr_estimation01.Rdata'))
get.result(betaest.mtx,lfsr.mtx)

There are 418 significant genes in total.

• We notice that gRNA targets with a large number of significant genes are the ones with relatively few cells (compared to the 111 negative control cells).

---

Permutation 2

load(paste0(wkdir,'data/perm.ashr_estimation02.Rdata'))
get.result(betaest.mtx,lfsr.mtx)

There are 374 significant genes in total.

---

Permutation 3

load(paste0(wkdir,'data/perm.ashr_estimation03.Rdata'))
get.result(betaest.mtx,lfsr.mtx)

There are 302 significant genes in total.

---

Permutation 4

load(paste0(wkdir,'data/perm.ashr_estimation04.Rdata'))
get.result(betaest.mtx,lfsr.mtx)

There are 357 significant genes in total.

---

Permutation 5

load(paste0(wkdir,'data/perm.ashr_estimation05.Rdata'))
get.result(betaest.mtx,lfsr.mtx)

There are 401 significant genes in total.

---

Session Information

sessionInfo()

R version 3.5.2 (2018-12-20)
Platform: x86_64-apple-darwin15.6.0 (64-bit)
Running under: macOS Mojave 10.14

Matrix products: default
BLAS: /System/Library/Frameworks/Accelerate.framework/Versions/A/Frameworks/vecLib.framework/Versions/A/libBLAS.dylib
LAPACK: /Library/Frameworks/R.framework/Versions/3.5/Resources/lib/libRlapack.dylib

locale:
[1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8

attached base packages:
 [1] grid      parallel  stats4    stats     graphics  grDevices utils    
 [8] datasets  methods   base     

other attached packages:
 [1] MAST_1.8.2                  SingleCellExperiment_1.4.1 
 [3] SummarizedExperiment_1.12.0 DelayedArray_0.8.0         
 [5] BiocParallel_1.16.6         matrixStats_0.54.0         
 [7] GenomicRanges_1.34.0        GenomeInfoDb_1.18.2        
 [9] edgeR_3.24.3                limma_3.38.3               
[11] lattice_0.20-38             sctransform_0.2.0          
[13] Matrix_1.2-16               clusterProfiler_3.10.1     
[15] org.Hs.eg.db_3.7.0          AnnotationDbi_1.44.0       
[17] IRanges_2.16.0              S4Vectors_0.20.1           
[19] Biobase_2.42.0              BiocGenerics_0.28.0        
[21] dplyr_0.8.0.1               gridExtra_2.3              
[23] kableExtra_1.0.1            knitr_1.22                 
[25] ggplot2_3.1.0               ashr_2.2-32                

loaded via a namespace (and not attached):
 [1] fgsea_1.8.0            colorspace_1.4-0       ggridges_0.5.1        
 [4] qvalue_2.14.1          XVector_0.22.0         rstudioapi_0.9.0      
 [7] listenv_0.7.0          farver_1.1.0           urltools_1.7.2        
[10] ggrepel_0.8.0          bit64_0.9-7            xml2_1.2.0            
[13] codetools_0.2-16       splines_3.5.2          pscl_1.5.2            
[16] doParallel_1.0.14      GOSemSim_2.8.0         polyclip_1.9-1        
[19] jsonlite_1.6           GO.db_3.7.0            ggforce_0.2.1         
[22] readr_1.3.1            compiler_3.5.2         httr_1.4.0            
[25] rvcheck_0.1.3          assertthat_0.2.0       lazyeval_0.2.1        
[28] tweenr_1.0.1           htmltools_0.3.6        prettyunits_1.0.2     
[31] tools_3.5.2            igraph_1.2.4           gtable_0.2.0          
[34] glue_1.3.0             GenomeInfoDbData_1.2.0 reshape2_1.4.3        
[37] DO.db_2.9              fastmatch_1.1-0        Rcpp_1.0.0            
[40] enrichplot_1.2.0       iterators_1.0.10       ggraph_1.0.2          
[43] xfun_0.5               stringr_1.4.0          globals_0.12.4        
[46] rvest_0.3.2            future_1.13.0          DOSE_3.8.2            
[49] zlibbioc_1.28.0        europepmc_0.3          MASS_7.3-51.1         
[52] scales_1.0.0           hms_0.4.2              RColorBrewer_1.1-2    
[55] yaml_2.2.0             memoise_1.1.0          UpSetR_1.3.3          
[58] triebeard_0.3.0        stringi_1.3.1          RSQLite_2.1.1         
[61] SQUAREM_2017.10-1      highr_0.7              foreach_1.4.4         
[64] truncnorm_1.0-8        rlang_0.3.1            pkgconfig_2.0.2       
[67] bitops_1.0-6           evaluate_0.13          purrr_0.3.1           
[70] labeling_0.3           cowplot_0.9.4          bit_1.1-14            
[73] tidyselect_0.2.5       plyr_1.8.4             magrittr_1.5          
[76] R6_2.4.0               DBI_1.0.0              pillar_1.3.1          
[79] withr_2.1.2            abind_1.4-5            RCurl_1.95-4.12       
[82] mixsqp_0.1-97          tibble_2.0.1           future.apply_1.2.0    
[85] crayon_1.3.4           rmarkdown_1.11         viridis_0.5.1         
[88] progress_1.2.0         locfit_1.5-9.1         data.table_1.12.0     
[91] blob_1.1.1             digest_0.6.18          webshot_0.5.1         
[94] tidyr_0.8.3            gridGraphics_0.3-0     munsell_0.5.0         
[97] viridisLite_0.3.0      ggplotify_0.0.3