Vignette illustrating the usage of gscreend on simulated data
2019-10-29
Contents
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)## Loading required package: GenomicRanges## Loading required package: stats4## Loading required package: BiocGenerics## Loading required package: parallel##
## Attaching package: 'BiocGenerics'## The following objects are masked from 'package:parallel':
##
## clusterApply, clusterApplyLB, clusterCall, clusterEvalQ,
## clusterExport, clusterMap, parApply, parCapply, parLapply,
## parLapplyLB, parRapply, parSapply, parSapplyLB## The following objects are masked from 'package:stats':
##
## IQR, mad, sd, var, xtabs## The following objects are masked from 'package:base':
##
## Filter, Find, Map, Position, Reduce, anyDuplicated, append,
## as.data.frame, basename, cbind, colnames, dirname, do.call,
## duplicated, eval, evalq, get, grep, grepl, intersect,
## is.unsorted, lapply, mapply, match, mget, order, paste, pmax,
## pmax.int, pmin, pmin.int, rank, rbind, rownames, sapply,
## setdiff, sort, table, tapply, union, unique, unsplit, which,
## which.max, which.min## Loading required package: S4Vectors##
## Attaching package: 'S4Vectors'## The following object is masked from 'package:base':
##
## expand.grid## Loading required package: IRanges## Loading required package: GenomeInfoDb## Loading required package: Biobase## Welcome to Bioconductor
##
## Vignettes contain introductory material; view with
## 'browseVignettes()'. To cite Bioconductor, see
## 'citation("Biobase")', and for packages 'citation("pkgname")'.## Loading required package: DelayedArray## Loading required package: matrixStats##
## Attaching package: 'matrixStats'## The following objects are masked from 'package:Biobase':
##
## anyMissing, rowMedians## Loading required package: BiocParallel##
## Attaching package: 'DelayedArray'## The following objects are masked from 'package:matrixStats':
##
## colMaxs, colMins, colRanges, rowMaxs, rowMins, rowRanges## The following objects are masked from 'package:base':
##
## aperm, apply, rowsum4 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.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...## 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.0006134969 0.000040 -1.421372
## essential_1 essential_1 0.0001347709 0.000008 -1.984253
## essential_10 essential_10 0.0004987531 0.000032 -1.605037
## essential_100 essential_100 0.0000000000 0.000000 -3.075151
## essential_1000 essential_1000 0.0000000000 0.000000 -2.023614
## essential_1001 essential_1001 0.0547085202 0.004880 -1.0326886 Session Info
sessionInfo()## R version 3.6.1 (2019-07-05)
## Platform: x86_64-pc-linux-gnu (64-bit)
## Running under: Ubuntu 18.04.3 LTS
##
## Matrix products: default
## BLAS: /home/biocbuild/bbs-3.10-bioc/R/lib/libRblas.so
## LAPACK: /home/biocbuild/bbs-3.10-bioc/R/lib/libRlapack.so
##
## locale:
## [1] LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C
## [3] LC_TIME=en_US.UTF-8 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
##
## attached base packages:
## [1] parallel stats4 stats graphics grDevices utils datasets
## [8] methods base
##
## other attached packages:
## [1] SummarizedExperiment_1.16.0 DelayedArray_0.12.0
## [3] BiocParallel_1.20.0 matrixStats_0.55.0
## [5] Biobase_2.46.0 GenomicRanges_1.38.0
## [7] GenomeInfoDb_1.22.0 IRanges_2.20.0
## [9] S4Vectors_0.24.0 BiocGenerics_0.32.0
## [11] gscreend_1.0.0 BiocStyle_2.14.0
##
## loaded via a namespace (and not attached):
## [1] Rcpp_1.0.2 nloptr_1.2.1 compiler_3.6.1
## [4] BiocManager_1.30.9 XVector_0.26.0 bitops_1.0-6
## [7] tools_3.6.1 timeSeries_3042.102 zlibbioc_1.32.0
## [10] digest_0.6.22 evaluate_0.14 lattice_0.20-38
## [13] rlang_0.4.1 fGarch_3042.83.1 Matrix_1.2-17
## [16] yaml_2.2.0 xfun_0.10 GenomeInfoDbData_1.2.2
## [19] stringr_1.4.0 knitr_1.25 grid_3.6.1
## [22] spatial_7.3-11 rmarkdown_1.16 bookdown_0.14
## [25] magrittr_1.5 fBasics_3042.89 htmltools_0.4.0
## [28] timeDate_3043.102 stringi_1.4.3 RCurl_1.95-4.12