About & methods

The data

The explorer shows the K562 essential-gene screen from Replogle et al. (2022, Cell): 2,176 CRISPRi knockdowns covering 2,062 genes, each compared with unperturbed control cells across 8,563 measured genes.

The experiment. Perturb-seq pairs a pooled CRISPR screen with single-cell RNA-seq. Each cell receives a guide RNA that knocks down one gene, and its whole transcriptome is read out together with the identity of that guide. One experiment therefore measures the downstream effects of thousands of knockdowns at once.

This screen. K562 is a chronic myeloid leukaemia cell line. The library targeted about 2,000 genes that are essential for growth across many cell lines, using CRISPR interference (CRISPRi), which silences a gene without cutting DNA. Some genes were targeted at two alternative start sites; these appear as separate entries, such as GAPDH_P1 and GAPDH_P2. The same study also ran a genome-wide K562 screen and an essential-gene screen in RPE1 cells, which are not shown here yet.

What the explorer starts from. Not raw cells, but the authors' published per-perturbation summary: for each knockdown, the average expression of every gene across its cells, scaled against control cells. A typical perturbation averages 116 cells (range 5 to 1,996).

Reading the results

• Knockdown is not knockout. The median knockdown is 83%, but 238 perturbations reduced their target by less than half. Check the knockdown figure before reading anything downstream.
• Many knockdowns change little. 78 perturbations have no significantly changed gene. That is often the true answer, not a fault of the page.
• Few cells, less trust. Perturbations measured in fewer than about 50 cells are noisier, and the significance calls for them are probably too generous (see Checks).
• Downstream is not direct. A changed gene may sit several regulatory steps away from the one knocked down.
• One cancer cell line, one time point. Effects may differ in other cell types, or at other times after knockdown.

Methods

Everything below is computed once, offline, by pipeline/03_summarise.py. The website itself does no statistics: it only draws the small files that script writes.

1. Starting point: the authors' normalised pseudobulk

The input file has one row per perturbation and one column per measured gene. Each value is the mean, over that perturbation's cells, of expression z-normalised against non-targeting control cells. A value of +1 means "on average, one control standard deviation above the control mean". The file also gives, for each gene, the control mean and standard deviation, and, for each perturbation, the number of cells and the target gene's remaining expression.

Why: this is the authors' own normalisation, so the explorer starts from the same numbers as the paper, and the file is small enough to process on a laptop.

2. Controls and perturbations

The 109 non-targeting control rows are set aside as the reference; the remaining 2,176 rows are perturbations. Genes targeted at two start sites keep one row per site rather than being merged, because the two guides can knock down the gene by different amounts.

3. Knockdown of the target

Knockdown = 100 × (1 − remaining expression), the reduction in the target gene relative to controls. The target's own column is then removed before any testing, so a perturbation never looks significant just because its own target went down.

4. Three ways to rank affected genes

The same change can be large on one scale and small on another, so the explorer stores three rankings and lets you switch between them.

Fold and absolute change are reconstructed by undoing the z-scaling with each gene's control mean \(\mu_g\) and standard deviation \(\sigma_g\), using a pseudocount \(c = 0.1\):

\[ \text{after}_g = \max(0,\ \mu_g + z_g\,\sigma_g), \qquad \text{lfc}_g = \log_2 \frac{\text{after}_g + c}{\mu_g + c}, \qquad \text{abs}_g = \text{after}_g - \mu_g \]

These are approximations: they assume one mean and one standard deviation describe the control cells for each gene, while the authors normalised each sequencing batch separately.

5. Which changes are significant

Each row is an average over \(n\) cells, so its noise shrinks as \(1/n\). The control rows, each also an average of a known number of cells, give an estimate of the single-cell variance of every gene \(g\):

\[ \hat{s}^2_g = \frac{1}{C-1} \sum_{c=1}^{C} n_c\,(x_{cg} - \mu_g)^2 \]

where \(C = 109\) is the number of control rows, \(n_c\) the number of cells in control row \(c\) and \(\mu_g\) the mean over control rows. A perturbation \(p\) measured in \(n_p\) cells then gets a test statistic

\[ t_{pg} = \frac{x_{pg} - \mu_g}{\sqrt{\hat{s}^2_g / n_p}} \]

compared with a standard normal distribution (two-sided). P-values are corrected with Benjamini-Hochberg across genes within each perturbation, and a gene is called changed at a 5% false discovery rate. A floor at the 1st percentile of the variances stops near-constant genes from producing huge statistics. On the explorer, significant genes are drawn as solid bars and the rest as faded bars.

Why this test: it needs only the published summary file, not the single cells, and it takes the number of cells into account, so a small shift seen in many cells and a large shift seen in a few are judged on the same scale.

6. What is sent to the browser

For each perturbation, the top 30 genes up and the top 30 down under each of the three rankings, merged, with all three values and a significance flag. Each file is a few kilobytes, and the page can re-rank without downloading anything else.

Checks

Controls should show nothing. Each control row in turn was tested against a reference built from the other 108. Most controls get no genes called (median 0), but the worst gets 41.

Agreement with the paper. The paper counted changed genes per perturbation with its own test on single cells (Anderson-Darling). Across all perturbations, the two counts rank almost identically (Spearman correlation 0.96), but ours are higher (median 84 genes against 62). The gap is concentrated in perturbations with few cells:

The likely cause: averaged over only a few cells, a sparsely expressed gene's value is far from normally distributed, so the normal test underestimates how often extreme values arise by chance. This is being checked (pipeline/04_check_null.py); until then, treat calls for perturbations with fewer than about 50 cells with caution.

Known limitations

• The test ignores batch structure, and any change in a gene's variance caused by the perturbation itself.
• With hundreds of cells, significant is not the same as large: check fold change and absolute change before reading biology into a significant gene.
• Fold changes are reconstructed from the summary file, not measured from counts.

Source and licence

Data: Replogle & Weissman, Figshare+, doi:10.25452/figshare.plus.20029387 (CC BY 4.0), file K562_essential_normalized_bulk_01.h5ad. Please cite Replogle et al. 2022 if you use these data.

Code: github.com/annacuomo/perturbseq-explorer (MIT licence).