Quality control, modeling, and visualization of CRISPR screens with MAGeCK-VISPR

Quality control measurements for CRISPR screening experiments

Apart from the determination of essential genes with MAGeCK, a central purpose of MAGeCK-VISPR is to collect quality control (QC) measurements at various levels (http://www.bioinformatics.babraham.ac.uk/projects/fastqc). The proposed measurements (Table 1) can be divided into four categories: sequence level, read count level, sample level, and gene level (Fig. 2).

QC term	Description	Expected
GC content	GC content distribution of the sequencing reads	Similar distribution for all samples from same library
Base quality	Base quality distribution of the sequencing reads	Single-peak distribution with median base quality at least 25
Sequencing reads	Total number of sequencing reads	Varies depending on sequencing platform
Mapped reads	Total number of reads mapped to the sgRNA library	300 * (number of sgRNAs)
% Mapped reads	Percentage of mapped reads to the total number of sequencing reads	At least 65 %
Zero sgRNAs	Number of sgRNAs with zero read counts	At most 1 % of total sgRNAs
Gini index	Gini index of log-scaled read count distributions	At most 0.1 for plasmid or initial state samples, and at most 0.2 for negative selection samples
Sample correlation	Pearson correlation coefficient between samples	At least 0.8 for replicates
Correlation clustering or PCA clustering	Hierarchical clustering of samples or first three PCA components	Samples with similar conditions should cluster together
Ribosomal gene selection	Negative selection enrichment statistics of ribosomal genes	Significant P values (<0.001) for ribosomal subunit (GO:0044391) in negative selection experiments

Sequence level QC measurements aim to detect problems with the sequencing, similar as in other next-generation sequencing (NGS) experiments. Two measurements are reported: sample GC content distribution (Fig. 2a) and the base quality distribution of sequencing reads (Fig. 2b, c). Ideally, sequencing reads should have reasonable base qualities (median value >25), and samples from the same experiment should have similar GC content distributions.

The second level of QC measurements is based on the sgRNA read counts collected from MAGeCK. Raw sequencing reads are first mapped to sgRNA sequences in the library with no mismatches tolerated. After that, the number of sequencing reads, mapped reads (and thereof the percentage of mapped reads), sgRNAs with zero read count, and the Gini index of read count distribution are reported for each sample (Fig. 2d-f). The percentage of mapped reads is a good indicator of sample quality, and low mappability could be due to sequencing error, oligonucleotide synthesis error, or sample contamination. Good statistical power of downstream analysis relies on sufficient reads (preferably over 300 reads) for each sgRNA, with low number of zero-count sgRNAs in the plasmid library or early time points. Gini index, a common measure of income inequality in economics, can measure the evenness of sgRNA read counts [14]. It is perfectly normal for later time points in positive selection experiments to have higher Gini index since a few surviving clones (a few sgRNA with extreme high counts) could dominate the final pool while most of the other cells die (more sgRNAs with zero-count). In contrast, high Gini index in plasmid library, in early time points, or in negative selection experiments may indicate CRISPR oligonucleotide synthesis unevenness, low viral transfection efficiency, and over selection, respectively.

Sample level QC (Fig. 2g-j) checks the consistency between samples. MAGeCK-VISPR reports the distributions of normalized read counts by box plots and cumulative distribution functions. It also calculates pairwise Pearson correlations of sample log read counts, and draws the samples on the first three components of a Principle Component Analysis (PCA). Biological replicates or samples with similar conditions should have similar read count distributions and higher correlations, and appear closer to each other in the PCA plot. PCA plots can also identify potential batch effects if the screens are conducted under different batches.

Finally, gene level QC determines the extent of negative selection in the screens. Since knocking out ribosomal genes lead to a strong negative selection phenotype [1, 2], the significance of negative selection on ribosomal genes can be evaluated in MAGeCK-VISPR by Gene Ontology (GO) enrichment analysis using GOrilla [15]. A working negative selection experiment should have a significant P value (<0.001), although many good experiments could have much smaller P values (<1e-10, see Section A of Additional file 1).

Calling essential genes under multiple conditions with MAGeCK-MLE

The effects of r different conditions are represented as the score ‘β _gr ’, a measurement of gene selections similar to the term of ‘log fold change’ in differential expression analysis. The presence or absence of each condition on each sample is encoded into binary elements of the design matrix d _jr , and can be obtained from experiment designs. ‘β’ scores reflect the extent of selection in each condition: β _gr >0 (or <0) means g is positively (or negatively) selected in condition r. μ _ij is also dependent on β _i0, the initial sgRNA abundance which is usually measured in plasmid or the day 0 of the experiment.

The values of β, together with the information whether an sgRNA is efficient, can be estimated by maximizing the joint log-likelihood of observing all sgRNA read counts of g on all different samples, and are optimized using an Expectation-Maximization (EM) algorithm. In the EM algorithm, MAGeCK-MLE iteratively determines the knockout efficiency of each sgRNA based on the current estimation of ‘β’ scores (the E step), and uses the updated knockout efficiency information to re-calculate ‘β’ scores (the M step). By examining the patterns of read counts of each sgRNA across all samples, the EM algorithm minimizes the effect of inefficient sgRNAs. A detailed description of the method is presented in the Methods section.

We tested MAGeCK algorithms on four public datasets. The first two datasets (the ‘ESC’ and ‘leukemia’ dataset) correspond to negative selection experiments on mouse embryonic stem cells (ESCs) and two human leukemia cell lines (KBM7 and HL-60), respectively (Fig. 3a and b) [1, 4]. In both datasets, cells were grown with their natural growing condition and negative selections occurred in cells after CRISPR/Cas9 is activated. The other two datasets (‘melanoma’ knockout and activation dataset) are different CRISPR screens on the human melanoma cell line A375 that harbors a BRAF V600E mutation (Figs. 4 and 5). The cells were treated with BRAF inhibitor vemurafenib (PLX) or dimethyl sulfoxide (DMSO) control, and screened either with GeCKO [2] or with CRISPR/dCas9 Synergistic Activation Mediator (SAM) libraries [5]. These two datasets include multiple experimental conditions that are difficult to compare directly using the original MAGeCK-RRA algorithm. In the melanoma knockout dataset, cells were under 7-day or 14-day selection [2]. In the melanoma activation dataset, two different drugs (puromycin and zeocin) were used to select cells with lentiviral infection, and both DMSO and PLX treatments were profiled under 3-day or 21-day selection [5].

In two-condition comparisons, MAGeCK-MLE gives similar results with existing methods such as MAGeCK-RRA, RSA, and RIGER. All the algorithms identified genes that are commonly essential to different cell types [16], as well as known positively selected genes in PLX treated conditions in two melanoma datasets (Fig. 3; also see Section A and B of Additional file 1). In the leukemia dataset, two-condition comparison algorithms (like MAGeCK-RRA) identified genes that are differentially selected in two cell lines by a direct comparison of HL60 and KBM7 (Fig. 3a) [10]. However, not all of these genes are equally biologically interesting, as MAGeCK-MLE further distinguished them into two groups: genes having little effect in one (β scores close to zero) but strong selection effect in the other cell line (large absolute β scores), and genes having weak and opposite effects in two cell lines (Fig. 3c). The first group of genes are often more biologically interesting as they are cell type-specific genes. This includes some well-known driver genes (like BCR in KBM7) as well as genes that may be functional in only one cell type: CDK6 and TRIB1 in HL60 [17, 18], and RUNX1 in KBM7 [19].

One of the advantages of MAGeCK-MLE over other methods is that it enables accurate comparisons of gene essentialities across multiple conditions and experiments in one run (Fig. 4 and Section C of Additional file 1). In the melanoma knockout dataset, a k-means clustering of the β scores of top selected genes demonstrated that these genes have various essentialities across conditions (Fig. 4a). Some of the genes are universally positively or negatively selected in all conditions (cluster 3), while others have different essentiality across different conditions (clusters 1, 2, and 4). Genes in cluster 4 are particularly interesting as they show strong positive selection in 14-day PLX treated condition. Indeed, genes whose knockout leads to strong positive selection in PLX-treated cells are in cluster 4, including NF1, NF2, MED12, CUL3 [2]. In contrast, the k-means clusters of measurements from other algorithms did not reveal the strong effect of genes in cluster 4 (Section C of Additional file 1). This is because their score distributions are similar across different conditions (Fig. 4b), and do not reflect the fact that the one condition (PLX 14-day treatment) induces much stronger positive selection than other conditions [2]. This is partly because MAGeCK-RRA, RIGER, and RSA all use a rank-based method to compare sgRNA between two conditions, which may lose quantitative information.

Another example of using MAGeCK-MLE on multiple conditions is demonstrated in the melanoma activation dataset, where cells underwent different selection methods (using puromycin or zeocin), drug treatments (DMSO or PLX), and durations (3-day or 21-day treatment) (Fig. 5). Similar to the melanoma knockout dataset, we performed k-means clustering of the top-selected gene β scores. Many positively selected genes are dependent on the selection method, which might not be biologically interesting. For example, genes in clusters 2 and 4 correspond to positively selected genes that are specific to puromycin or zeocin selection, respectively. A small set of genes (cluster 5) are consistently selected in both zeocin and puromycin, including genes that are validated in the original study, for example, EGFR, GPR35, LPAR1/5 [5]. We further examined the genes in cluster 5 (Fig. 5b), and focused on genes positively selected in the CRISPR activation experiment but strongly negatively selected in the knockout experiment. These genes include EGFR and BRAF, two known kinases that drive melanoma progression and PLX resistance [20, 21], and CRKL, a protein kinase that activates RAS and JUN pathway. CRKL amplification is reported to lead to drug resistance against EGFR inhibitors by activating EGFR downstream pathways [22], implying its potential role in PLX drug resistance.

Visualization of QC measurements and gene essentiality with VISPR

VISPR (VISualization of crisPR screens) is a web-based frontend for interactive visualization of CRISPR screen QC and comparison results. Interactive access is provided by an HTML5 based browser interface, while visualizations are realized with Vega [23], a declarative visualization grammar on top of Data-Driven Documents (D3) [24]. VISPR provides three types of views for interactive exploration of CRISPR screening: a quality control view, a result view, and an experiment comparison view. The quality control view shows the QC measurements described before (Fig. 2).

In the result view, screening results can be interactively explored. It contains a table showing the comparison results of each gene (Fig. 6a). The table can be sorted by different columns and filtered (from ‘Search’) via gene names or regular expressions. Further, the distribution of P values is displayed as cumulative distribution function (CDF) (Fig. 6b) and as a histogram (Fig. 6c). For each gene, the normalized sgRNA counts in all samples can be displayed in a parallel coordinate visualization (Fig. 6d). If available, knockout efficiency predictions [11] and gene coordinates of each sgRNA are displayed as separate axes. Axes can be reordered or toggled on or off, and sgRNAs can be highlighted by selecting ranges on each axis. Genes selected in the table are highlighted in the CDF, allowing to assess their occurrence within the P value distribution of all genes.

VISPR provides various ways to further explore the analysis results. Individual genes can be viewed in Ensembl [25] and IGV [26]. Selected genes can be visualized in terms of their interaction network and function via GeneMANIA [27]. Functional analysis can be performed with GOrilla [15], an online Gene Ontology (GO) enrichment analysis tool. GOrilla takes a ranked list of genes (here based on the P values reported by MAGeCK) to perform a threshold-free enrichment analysis. The resulting GO term enrichments can be further used for gene-level quality control.

The comparison view of VISPR can compare different experiments by visualizing the common and exclusive significant genes via Euler diagrams (Fig. 6e). Clicking on segments of the Euler diagram opens the result views of the corresponding experiments. For example, clicking on the intersection between two experiments will open ‘restricted’ result views for each experiment, where only the common significant genes are displayed. These views provide the same features as the unconstrained result views described above. However, in this case, GO enrichment analysis with GOrilla is performed with the shown genes (that is, the genes from the intersection) as foreground and the other genes of the experiment as background.

The visualizations displayed in VISPR can be downloaded as publication-ready SVG files. In addition, a command line interface is provided to store visualizations as Vega specifications. This format allows users to modify and style the output of VISPR programmatically.

Implementation of the MAGeCK-VISPR workflow with Snakemake

We implemented the MAGeCK-VISPR workflow with the workflow management system Snakemake [13], allowing an automatic execution of some or all of the MAGeCK-VISPR functions: quality control, essential gene analysis, and visualization. Choosing a workflow management system like Snakemake has several advantages. First, the workflow steps can be automatically parallelized and executed on workstations, servers, and compute clusters without the modification of the workflow. Second, Snakemake tracks metadata (like creation date, input, and log files) for all generated result and intermediate files. This way, used data, methods, and parameters are documented comprehensively for each analysis (also called data provenance), an important requirement of reproducible science. MAGeCK-VISPR provides a command line interface to initialize the workflow in a given work directory. This installs the workflow definition as a so-called Snakefile, along with a configuration file and documentation. The configuration file is used to define locations of raw data and additional parameters for MAGeCK-VISPR. Once configured, the Snakefile can be executed with Snakemake. Since the Snakefile is installed into the given work directory, it can be easily modified or extended by the user.

We provide all components of the workflow as Conda packages [28], such that MAGeCK-VISPR can be installed with a single command. Optionally, the Conda package manager can create isolated environments for the workflow to, for example, freeze or compare different software versions or publish snapshots of a MAGeCK-VISPR workflow instance along with all data and used software. This further increases the reproducibility of the generated results.