1 Introduction

Pooled CRISPR perturbations screens employ a library of guide RNAs (gRNAs) that is transduced into a pool of cells with the aim to induce a single genetic perturbation in each cell. The perturbation effect is assessed by measuring the abundance of each gRNA after the screen selection phase and comparing it to its abundance in the plasmid library. The main goal of the following analysis is the detection of essential genes, i.e. genes whose knockout reduces the cell fitness. The package gscreend provides a method to rank genes based on count tables.

2 gscreend workflow

In order to identify essential genes starting from raw gRNA count data, gscreend performs the following analysis steps:

Input of raw gRNA counts at T0 (sequencing of library) and T1 (at the end of the screen). Normalization and calculation of log fold changes.
Split log fold changes into intervals dependent on the initial count at T0.
For every interval fit a skew-normal distribution to the data to model the null hypothesis (via least quantile regression).
Based on the null model calculate p-values for every gRNA.
Rank gRNAs according to p-value and perform robust ranking aggregation to calculate p-values on gene level.
Perform quality control of data and statistical model.

3 Installation

if(!requireNamespace("BiocManager", quietly = TRUE))
    install.packages("BiocManager")
BiocManager::install("gscreend")

library(gscreend)
library(SummarizedExperiment)

4 Analysis of simulated data with gscreend

4.1 Input data: gRNA counts

The simulated data used in this example has been generated using the simulation method available at https://github.com/imkeller/simulate_pooled_screen

Raw count data consists of gRNA counts in the library sequencing and different replicates after the screen proliferation phase. In order to estimate the effect of a specific gRNA on cell fitness, the relative abundances of the gRNA before and after the proliferation phase will be compares

raw_counts <- read.table(
    system.file("extdata", "simulated_counts.txt", package = "gscreend"),
    header=TRUE)

Generate a summarized experiment from the count data. gscreend currently uses SummarizedExperiment objects as an input format.

The count matrix contains raw gRNA counts. Each row represents one gRNA, each column represents one sample (T0, T1, replicates, …).

counts_matrix <- cbind(raw_counts$library0, 
                        raw_counts$R0_0, 
                        raw_counts$R1_0)

rowData <- data.frame(sgRNA_id = raw_counts$sgrna_id,
                    gene = raw_counts$Gene)

colData <- data.frame(samplename = c("library", "R1", "R2"),
                    # timepoint naming convention: 
                    # T0 -> reference, 
                    # T1 -> after proliferation
                    timepoint = c("T0", "T1", "T1"))

se <- SummarizedExperiment(assays=list(counts=counts_matrix),
                    rowData=rowData, colData=colData)

4.2 Run gscreend

In this step a gscreend experiment object is generated that will after the analysis contain all data related to gRNAs, genes and model parameters.

pse <- createPoolScreenExp(se)

## Creating PoolScreenExp object from a SummarizedExperiment object.

## References and samples are named correctly.

## Data concerning sgRNA and genes is provided.

Run gscreend with default parameters.

pse_an <- RunGscreend(pse)

## Size normalized count data.

## Calculated LFC.

## Fitted null distribution.

## Calculated p-values at gRNA level.

## Ranking genes...

## ... for positive fold changes

## ... for negative fold changes

## gscreend analysis has been completed.

4.3 Quality control

gscreend provides basic quality control functions for inspection of replicate correlation for example.

plotReplicateCorrelation(pse_an)

The plotModelParameters() function can be used to inspect the values of the parameters estimated for the skew normal distribution of the logarithmic fold change data.

plotModelParameters(pse_an)

5 Results

The ResultsTable function can be used to extract a table listing for each gene the p-value and fdr. These values correspond to the results from the statistical test indicting whether upon perturbation a specific gene reduces (direction = “negative”), or increasing (direction = “positive”) cell viability.

res <- ResultsTable(pse_an, direction = "negative")
head(res)

##                          Name          fdr    pval       lfc
## essential_0       essential_0 0.0003075031 0.00002 -1.421372
## essential_1       essential_1 0.0000000000 0.00000 -1.984253
## essential_10     essential_10 0.0001554726 0.00001 -1.605037
## essential_100   essential_100 0.0000000000 0.00000 -3.075151
## essential_1000 essential_1000 0.0000000000 0.00000 -2.023614
## essential_1001 essential_1001 0.0544843049 0.00486 -1.032688

6 Session Info

sessionInfo()

## R version 4.3.0 RC (2023-04-13 r84269)
## Platform: x86_64-pc-linux-gnu (64-bit)
## Running under: Ubuntu 22.04.2 LTS
## 
## Matrix products: default
## BLAS:   /home/biocbuild/bbs-3.17-bioc/R/lib/libRblas.so 
## LAPACK: /usr/lib/x86_64-linux-gnu/lapack/liblapack.so.3.10.0
## 
## locale:
##  [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C              
##  [3] LC_TIME=en_GB              LC_COLLATE=C              
##  [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       
## 
## time zone: America/New_York
## tzcode source: system (glibc)
## 
## attached base packages:
## [1] stats4    stats     graphics  grDevices utils     datasets  methods  
## [8] base     
## 
## other attached packages:
##  [1] SummarizedExperiment_1.30.0 Biobase_2.60.0             
##  [3] GenomicRanges_1.52.0        GenomeInfoDb_1.36.0        
##  [5] IRanges_2.34.0              S4Vectors_0.38.0           
##  [7] BiocGenerics_0.46.0         MatrixGenerics_1.12.0      
##  [9] matrixStats_0.63.0          gscreend_1.14.0            
## [11] BiocStyle_2.28.0           
## 
## loaded via a namespace (and not attached):
##  [1] sass_0.4.5              bitops_1.0-7            fGarch_4022.89         
##  [4] lattice_0.21-8          magrittr_2.0.3          digest_0.6.31          
##  [7] evaluate_0.20           grid_4.3.0              bookdown_0.33          
## [10] fastmap_1.1.1           jsonlite_1.8.4          Matrix_1.5-4           
## [13] BiocManager_1.30.20     gbutils_0.5             codetools_0.2-19       
## [16] jquerylib_0.1.4         Rdpack_2.4              cli_3.6.1              
## [19] timeSeries_4021.105     rlang_1.1.0             rbibutils_2.2.13       
## [22] XVector_0.40.0          cachem_1.0.7            DelayedArray_0.26.0    
## [25] yaml_2.3.7              cvar_0.5                tools_4.3.0            
## [28] parallel_4.3.0          BiocParallel_1.34.0     nloptr_2.0.3           
## [31] GenomeInfoDbData_1.2.10 R6_2.5.1                magick_2.7.4           
## [34] zlibbioc_1.46.0         bslib_0.4.2             Rcpp_1.0.10            
## [37] xfun_0.39               highr_0.10              knitr_1.42             
## [40] spatial_7.3-16          htmltools_0.5.5         fBasics_4022.94        
## [43] rmarkdown_2.21          timeDate_4022.108       compiler_4.3.0         
## [46] RCurl_1.98-1.12

Vignette illustrating the usage of gscreend on simulated data

2023-04-25

Contents