ashr on sctransformed CROP-seq data

Data property

nlocus <- colSums(exp.per.enhancer>0)
gene.percent <- rowMeans(gene.exp>0)
hist(gene.percent,xlab = '% of cells a gene is in',ylab = 'Count')

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
cat('Number of cells uniquely targeted by each gRNA:')

Number of cells uniquely targeted by each gRNA:

print(ncell_uniq)

        BAG5_1         BAG5_2         BAG5_3       BCL11B_1       BCL11B_2 
            36             33             31             13             33 
      BCL11B_3       CHRNA3_1       CHRNA3_2       CHRNA3_3      GALNT10_1 
            13              3             23             13             38 
     GALNT10_2      GALNT10_3       KCTD13_1       KCTD13_2       KCTD13_3 
            37             50             30             32             15 
       KMT2E_1        KMT2E_2        KMT2E_3 LOC105376975_1 LOC105376975_2 
            11             19             18             35             50 
LOC105376975_3       miR137_1       miR137_2       miR137_3         NAB2_1 
             8             44             30             54             36 
        NAB2_2         NAB2_3     PCDHA123_1     PCDHA123_2     PCDHA123_3 
            24             25             53             64             62 
        RERE_1         RERE_2         RERE_3     SETD1A_1st     SETD1A_1st 
            48             46             25             23             28 
    SETD1A_1st        STAT6_1        STAT6_2        STAT6_3       TRANK1_1 
             4             24             38             53              3 
      TRANK1_2       TRANK1_3     UBE2Q2P1_1     UBE2Q2P1_2     UBE2Q2P1_3 
            19             65             11              8             11 
       VPS45_1        VPS45_2        VPS45_3         NGEF_1         NGEF_2 
            14             79              9             11              2 
        NGEF_3         DPYD_1         DPYD_2         DPYD_3 LOC100507431_1 
             0             37             68             28             15 
LOC100507431_2 LOC100507431_3    LINC00637_1    LINC00637_2    LINC00637_3 
            29             49             45             34             53 
       PBRM1_1        PBRM1_2        PBRM1_3     PPP1R16B_1     PPP1R16B_2 
             4             12              4             44             31 
    PPP1R16B_3    TCF4-ITF2_1    TCF4-ITF2_2    TCF4-ITF2_3     pos_SNAP91 
            35             50             45             12              3 
    pos_SNAP91       neg_EGFP       neg_EGFP       neg_EGFP       neg_CTRL 
            26              5             29             35             17 
      neg_CTRL 
            26 

Categorical regression on sctransformed data

We transformed the single-cell raw count data using sctransform to Pearson residuals that supposedly follow \(N(0,1)\) distribution per gene.

We filtered the genes to be those present in at least 20% of all cells (8117 in total).

Then, 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. We conducted categorical regression on the Pearson residuals (~ 2 conditions) to obtain the summary statistics reflecting how differentially expressed each gene is under the 2 conditions. 70 (gRNA vs neg ctrl) pairs were tested.

categoric <- function(mtx,tg.cells,neg.cells){
  pval <- rep(NA,nrow(mtx))
  beta <- rep(NA,nrow(mtx))
  se <- rep(NA,nrow(mtx))
  subset <- mtx[,c(tg.cells,neg.cells)]
  for (i in 1:nrow(subset)){
    df <- data.frame(row.names = c(tg.cells,neg.cells),
                     residual = as.numeric(subset[i,]),
                     condition=c(rep(1,length(tg.cells)),
                                 rep(0,length(neg.cells))))
    cat.lm <- summary(lm(residual~factor(condition), data = df))
    beta[i] <- cat.lm$coefficients[2,1]
    se[i] <- cat.lm$coefficients[2,2]
    pval[i] <- cat.lm$coefficients[2,4]
  }
  return(list(beta=beta,se=se,pval=pval))
}
vst_out <- readRDS(paste0(wkdir,'data/sctransformed_object.rds'))
sctransform_data <- vst_out$y
## Filter the genes to be present in > 20% cells:
subset.genes <- names(gene.percent)[gene.percent>0.2]
subset.sctransform_data <- sctransform_data[subset.genes,]
dim(subset.sctransform_data)
## Categorical regression on sctransformed data 
summary_per_gene <- list()
neg.cells <- colnames(exp.per.enhancer)[exp.per.enhancer["neg",]>0 & nlocus==1]
# btw each gRNA condition and neg ctrl
for (glocus in row.names(exp.per.gRNA)[-c(51,72:76)]){
  print(paste('gRNA target:',glocus))
  tg.cells <- colnames(exp.per.gRNA)[exp.per.gRNA[glocus,]>0 & nlocus==1]
  print(paste('# of targeted cells:',length(tg.cells)))
  summary_per_gene[[glocus]] <- categoric(subset.sctransform_data,tg.cells,neg.cells)
}
saveRDS(summary_per_gene,file = 'data/categoric_summary.sctrans_gene0.2.rds')

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).

summary_per_gene <- readRDS(paste0(wkdir,'data/categoric_summary.sctrans_gene0.2.rds'))
beta.mtx <- sapply(summary_per_gene, function(x){x$beta})
se.mtx <- sapply(summary_per_gene, function(x){x$se})
pval.mtx <- sapply(summary_per_gene, function(x){x$pval})
subset.genes <- names(gene.percent)[gene.percent>0.2]
rownames(beta.mtx) <- subset.genes
rownames(se.mtx) <- subset.genes
rownames(pval.mtx) <- subset.genes

lfsr.mtx <- data.frame(matrix(nrow = nrow(beta.mtx),ncol = ncol(beta.mtx)),
                       row.names = rownames(beta.mtx))
names(lfsr.mtx) <- colnames(beta.mtx)
betaest.mtx <- lfsr.mtx
for (i in 1:nrow(beta.mtx)){
  ash.tmp <- ash(beta.mtx[i,],se.mtx[i,])
  lfsr.mtx[i,] <- get_lfsr(ash.tmp)
  betaest.mtx[i,] <- get_pm(ash.tmp)
}
save(betaest.mtx,lfsr.mtx,file = 'data/ashr_estimation.Rdata')

load(paste0(wkdir,'data/ashr_estimation.Rdata'))
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.')

There are 462 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:
kable(compact.genes.df) %>%
  kable_styling() %>%
  scroll_box(width = "100%", height = "400px")

num_gRNA_pergene <- sapply(signif_genes.gRNA.lst,FUN = function(x){length(x)})
hist(num_gRNA_pergene,xlab = 'Number of gRNA targets per gene',ylab = 'Count',
     main = 'Distribution of 462 significant genes')

library(gridExtra)
plot_beta <- function(gene){
  num <- num_gRNA_pergene[gene]
  log10_lfsr <- log10(as.numeric(lfsr.mtx[gene,]))
  beta.plot <- qplot(as.numeric(beta.mtx[gene,]),as.numeric(betaest.mtx[gene,]),
                     color=log10_lfsr, main = paste0(gene,' (',num,')'),
                     xlab = 'beta from regression',
                     ylab = 'beta after shrinkage') + 
               geom_abline(intercept = 0,slope = 1,linetype = "dotted")
  return(beta.plot)
}

plot1 <- plot_beta('ARAF')
plot2 <- plot_beta('CWF19L2')
plot3 <- plot_beta('THG1L')
plot4 <- plot_beta('ZNF48')
grid.arrange(plot1,plot2,plot3,plot4,ncol=2)

plot1 <- plot_beta('A2M')
plot2 <- plot_beta('SUMO3')
plot3 <- plot_beta('DNAJC9')
plot4 <- plot_beta('APC')
grid.arrange(plot1,plot2,plot3,plot4,ncol=2)

plot_zscore <- function(gene){
  num <- num_gRNA_pergene[gene]
  zscore <- as.numeric(beta.mtx[gene,])/as.numeric(se.mtx[gene,])
  zscore.df <- data.frame(score=zscore)
  zscore.plot <- ggplot(zscore.df,aes(zscore)) + geom_histogram(aes(y=..density..),bins = 20) +
    labs(x = 'z-score from regression',y='frequency',title = paste0(gene,' (',num,')')) +
    stat_function(fun = dnorm, args = list(mean = 0, sd = 1),col=2)
  return(zscore.plot)
}
plot1 <- plot_zscore('ARAF')
plot2 <- plot_zscore('CWF19L2')
plot3 <- plot_zscore('THG1L')
plot4 <- plot_zscore('ZNF48')
grid.arrange(plot1,plot2,plot3,plot4,ncol=2)

plot1 <- plot_zscore('A2M')
plot2 <- plot_zscore('SUMO3')
plot3 <- plot_zscore('DNAJC9')
plot4 <- plot_zscore('APC')
grid.arrange(plot1,plot2,plot3,plot4,ncol=2)

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:
kable(compact.gRNAs.df) %>%
  kable_styling() %>%
  scroll_box(width = "100%", height = "400px")

num_gene_pergRNA <- sapply(gRNA.signif_genes.lst,FUN = function(x){length(x)})
hist(num_gene_pergRNA,xlab = 'Number of DE genes per gRNA target',ylab = 'Count',
     main = 'Distribution of 70 gRNAs')

gene_info <- read.delim(paste0(wkdir,'data/genes.tsv'),header = F,
                        sep = '\t',col.names = c('ID','name'))
gene_info$ID <- sapply(strsplit(gene_info$ID,split = '[.]'),
                       FUN = function(x){x[[1]]})
tmp.indx <- match(signif_genes,gene_info$name)
ensembl.signif_genes <- gene_info$ID[tmp.indx]

library(org.Hs.eg.db)
library(clusterProfiler)
ego.bp.signif_genes <- enrichGO(gene       = ensembl.signif_genes,
                                universe      = gene_info$ID,
                                OrgDb         = org.Hs.eg.db,
                                keyType       = 'ENSEMBL',
                                ont           = "BP",
                                pAdjustMethod = "BH",
                                pvalueCutoff  = 0.01,
                                qvalueCutoff  = 0.05,
                                readable      = TRUE)
barplot(ego.bp.signif_genes, drop=TRUE, showCategory=6)

kable(head(ego.bp.signif_genes),row.names = F) %>%
  kable_styling() %>%
  scroll_box(width = "100%", height = "300px")

ego.cc.signif_genes <- enrichGO(gene       = ensembl.signif_genes,
                                universe      = gene_info$ID,
                                OrgDb         = org.Hs.eg.db,
                                keyType       = 'ENSEMBL',
                                ont           = "CC",
                                pAdjustMethod = "BH",
                                pvalueCutoff  = 0.01,
                                qvalueCutoff  = 0.05,
                                readable      = TRUE)
barplot(ego.cc.signif_genes, drop=TRUE, showCategory=6)

kable(head(ego.cc.signif_genes),row.names = F) %>%
  kable_styling() %>%
  scroll_box(width = "100%", height = "300px")

gene.num <- sapply(signif_genes.gRNA.lst,FUN = function(x){length(x)})
names(gene.num) <- ensembl.signif_genes
gene.num <- gene.num[order(gene.num,decreasing = T)]

gsea.signif_genes <- gseGO(geneList     = gene.num,
                           OrgDb        = org.Hs.eg.db,
                           keyType      = 'ENSEMBL',
                           ont          = "BP",
                           nPerm        = 1000,
                           minGSSize    = 100,
                           maxGSSize    = 500,
                           pvalueCutoff = 0.05,
                           verbose      = FALSE)
dotplot(gsea.signif_genes)

# convert the results back to gene names
convert <- function(x){
  indx <- match(x,gene_info$ID)
  names <- gene_info$name[indx]
  return(paste(names,collapse = '/'))
}
tmp <- strsplit(gsea.signif_genes@result$core_enrichment,split = '/')
tmp1 <- sapply(tmp,FUN = convert)
gsea.signif_genes@result$core_enrichment <- tmp1
kable(head(gsea.signif_genes),row.names = F) %>%
  kable_styling() %>%
  scroll_box(width = "100%", height = "400px")

Permutation test

Please see this page for results of permutation tests.

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: /Library/Frameworks/R.framework/Versions/3.5/Resources/lib/libRblas.0.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] parallel  stats4    stats     graphics  grDevices utils     datasets 
[8] methods   base     

other attached packages:
 [1] clusterProfiler_3.10.1 org.Hs.eg.db_3.7.0     AnnotationDbi_1.44.0  
 [4] IRanges_2.16.0         S4Vectors_0.20.1       Biobase_2.42.0        
 [7] BiocGenerics_0.28.0    gridExtra_2.3          kableExtra_1.0.1      
[10] knitr_1.22             ggplot2_3.1.0          ashr_2.2-32           

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