README.md

metadata

title: Tutorial - Autorep sgRNA strength from pooled CRISPRi-seq and RNA-seq
emoji: 🧮
colorFrom: green
colorTo: blue
sdk: gradio
sdk_version: 6.30.0
app_file: app.py
pinned: true
license: mit
short_description: Inferring sgRNA strength from NGS autorep experiments
tags:
  - biology
  - sequencing
  - pooled-screen
  - fitness
  - gradio

Tutorial – Autorep sgRNA strength from pooled CRISPRi-seq and RNA-seq

This is the repository for the interactive tutorial, acccessed here. Non-interactive text from the tutorial is below.

The autorep experiment

CRISPRi sgRNAs differ in how strongly they repress their target, and the difference depends on the spacer sequence. How do we rank sgRNAs by strength, in a pool?

In the autorep design, every strain carries a fluorescent reporter (mWasabi) with the sgRNA's own target site in its 5'UTR. Each strain has one sgRNA, so the spacer is also the strain's barcode. The strains are grown together in a pool. When dCas9 is induced, each sgRNA represses its own reporter, and we measure two things by sequencing:

• $c_i(t)$: gDNA spacer-amplicon counts for strain $i$ at time $t$. These count cells.
• $x_i(t)$: RNA-seq counts of the barcoded reporter transcript for strain $i$ at time $t$. These count transcripts.

A dCas9 bound at the target site blocks transcription. If a single site is at equilibrium occupancy, the fraction of time the promoter is free (unrepressed) is

$$\frac{1}{1+K_{i}R}\backslash frac\{1\}\{1\; +\; K\_i\; R\}$$

• $K_i$: affinity of the dCas9–sgRNA complex for its site. This is the sgRNA strength we want to rank.
• $R$: concentration of free dCas9. It is the same for every strain.

So production of reporter transcript from strain $i$ is

$$E_i=\frac{E_{max}}{1+K_{i}R}E\_i\; =\; \backslash frac\{E\_\{max\}\}\{1\; +\; K\_i\; R\}$$

where $E_{max}$ is production with no dCas9 bound. Because $R$ is shared, we can only identify the dimensionless product $K_i R$. That is enough for ranking, since multiplying every $K_i$ by the same constant does not change the order.

Transcript dynamics in growing cells

Let $M_i$ be the total number of reporter transcripts from strain $i$ in the pool, and $n_i$ the number of cells. Transcripts are made at rate $E_i$ per cell and decay at rate $d$:

$$\frac{d{M}_i}{dt}=E_i{n}_i-d{M}_i\backslash frac\{dM\_i\}\{dt\}\; =\; E\_i\; n\_i\; -\; d\; M\_i$$

This is awkward to work with because $M_i$ also grows whenever the population grows. We want the transcript level per cell, $m_i = M_i / n_i$. Differentiating $M_i / n_i$ with the quotient rule, and substituting $\frac{dM_i}{dt}$ from above, $$\frac{d}{dt}\frac{M_i}{n_i}=\frac{n_i(E_i{n}_i-d{M}_i)-M_i\frac{d{n}_i}{dt}}{n_i^2}=E_i-d{m}_i-\left(\frac{1}{n_i}\frac{d{n}_i}{dt}\right)m_i\backslash frac\{d\}\{dt\}\; \backslash frac\{M\_i\}\{n\_i\}\; =\; \backslash frac\{n\_i\; \backslash left(\; E\_i\; n\_i\; -\; d\; M\_i\; \backslash right)\; -\; M\_i\; \backslash frac\{dn\_i\}\{dt\}\}\{n\_i^2\}\; =\; E\_i\; -\; d\; \backslash ,\; m\_i\; -\; \backslash left(\; \backslash frac\{1\}\{n\_i\}\; \backslash frac\{dn\_i\}\{dt\}\; \backslash right)\; m\_i$$

The bracketed term is the instantaneous per-capita growth rate of the strain. Call it $r_i$:

$$r_i=\frac{1}{n_i}\frac{d{n}_i}{dt}=\frac{d\log n_i}{dt}r\_i\; =\; \backslash frac\{1\}\{n\_i\}\; \backslash frac\{dn\_i\}\{dt\}\; =\; \backslash frac\{d\; \backslash log\; n\_i\}\{dt\}$$

Then

$$\frac{d{m}_i}{dt}=E_i-(d+r_i)m_i\backslash frac\{dm\_i\}\{dt\}\; =\; E\_i\; -\; (d\; +\; r\_i)\; m\_i$$

The extra term $r_i m_i$ is dilution. Every division splits a cell's transcripts between the two daughters, so growth removes transcripts per cell in the same way decay does.

• $m_i(t)$: reporter transcripts per cell in strain $i$.
• $d$: mRNA decay rate. For bacterial mRNA, half-lives of 2–5 minutes correspond to $d \approx$ 0.15–0.35 per minute.
• $r_i(t)$: instantaneous per-capita growth rate of strain $i$. In exponential phase it equals the intrinsic growth rate $w_i$ of the previous tutorial. In the Lotka–Volterra model it is $w_i$ multiplied by $\left(1 - \frac{\Sigma_j n_j}{K}\right)$, a factor shared by all strains. It can be measured from the DNA counts as $d \log n_i / dt$, with the spike-in putting $n_i(t)$ on an absolute scale (see the previous tutorial).

This equation has a closed-form solution if $r_i$ is constant over the relaxation. If $m_i$ starts at $m_i(0)$,

$$m_i(t)=m_i^{\ast }+(m_i(0)-m_i^{\ast })e^{-(d+r_i)t},m_i^{\ast }=\frac{E_i}{d+r_i}m\_i(t)\; =\; m\_i^*\; +\; \backslash left(\; m\_i(0)\; -\; m\_i^*\; \backslash right)\; e^\{-(d\; +\; r\_i)\; t\},\; \backslash qquad\; m\_i^*\; =\; \backslash frac\{E\_i\}\{d\; +\; r\_i\}$$

The transcript level relaxes to its steady state $m_i^*$ at rate $d + r_i$. With a decay rate of 0.15–0.35 per minute, this takes a few minutes.

Steady state

Samples taken tens of minutes or hours after induction are on the plateau, so we can use $m_i = m_i^*$ (the quasi-steady state, QSS). This holds when $r_i(t)$ changes slowly compared with the relaxation rate $d + r_i$, so that $m_i$ tracks its instantaneous steady state $E_i / (d + r_i(t))$. Below, $r_i$ means its value at the sampling time. Compare strain $i$ to a reference strain 0 with no repression ($K_0 = 0$, e.g. a non-targeting sgRNA):

$$\frac{m_i}{m_0}=\frac{d+r_0}{d+r_i}\cdot \frac{1}{1+K_{i}R}\backslash frac\{m\_i\}\{m\_0\}\; =\; \backslash frac\{d\; +\; r\_0\}\{d\; +\; r\_i\}\; \backslash cdot\; \backslash frac\{1\}\{1\; +\; K\_i\; R\}$$

The first factor is a growth term, evaluated with the instantaneous growth rates at the sampling time. It is 1 when $r_i = r_0$, and then

$$\frac{m_i}{m_0}=\frac{1}{1+K_{i}R}\backslash frac\{m\_i\}\{m\_0\}\; =\; \backslash frac\{1\}\{1\; +\; K\_i\; R\}$$

So transcript level per cell falls with sgRNA strength, with a known shape, and no time dependence under steady state (late timepoint) consitions.

We assume the growth term is 1 for the analysis below. This holds near carrying capacity, where every strain grows slowly and similarly. It fails for guides that have a strong effect on an essential chromosomal gene. Those strains grow slower, dilute their transcripts less, and look weaker than they are.

It might seem better to fit the time course after induction. But dCas9 does not switch on instantly: $R(t)$ rises over a timescale comparable to or slower than mRNA relaxation. The early transient therefore mostly reports how fast dCas9 accumulates, not $d$ or $K_i$. Fitting it with a model that assumes a step change in $R$ misattributes the rates. For ranking, use late timepoints only.

Read counts from next-generation sequencing

But we don't measure $m_i$ directly. Each sample is a random sampling of the pool, followed by molecular biology handling and uneven sequencing depth. Each guide also has its own capture and amplification efficiency, which differs between strains because the 5'UTR sequence differs.

For one guide in sample $s$, relative to its own $t = 0$ sample,

$$\frac{c_i(t)}{c_i(0)}={\varphi }_s\frac{n_i(t)}{n_i(0)}\frac{x_i(t)}{x_i(0)}=k_s\frac{M_i(t)}{M_i(0)}=k_s\frac{m_i(t)}{m_i(0)}\frac{n_i(t)}{n_i(0)}\backslash frac\{c\_i(t)\}\{c\_i(0)\}\; =\; \backslash phi\_s\; \backslash frac\{n\_i(t)\}\{n\_i(0)\}\; \backslash qquad\; \backslash frac\{x\_i(t)\}\{x\_i(0)\}\; =\; k\_s\; \backslash frac\{M\_i(t)\}\{M\_i(0)\}\; =\; k\_s\; \backslash frac\{m\_i(t)\}\{m\_i(0)\}\; \backslash frac\{n\_i(t)\}\{n\_i(0)\}$$

Variables:

• $c_i(t)$: gDNA read count of strain $i$ at time $t$.
• $x_i(t)$: RNA read count of strain $i$ at time $t$.
• $\phi_s$: ratio of DNA sampling depth at time $t$ to that at time 0 for sample $s$.
• $k_s$: the same for RNA-seq.

Dividing by $t = 0$ cancels each guide's constant capture efficiency and any per-guide mRNA stability effect that is the same before and after induction. Dividing the RNA ratio by the DNA ratio cancels the number of cells, $\frac{n_i(t)}{n_i(0)}$, which contains all the population dynamics. Define

$$A_i=\frac{x_i(t)\mathrm{/}{c}_i(t)}{x_i(0)\mathrm{/}{c}_i(0)}=\frac{k_s}{{\varphi }_s}\frac{m_i(t)}{m_i(0)}A\_i\; =\; \backslash frac\{x\_i(t)\; /\; c\_i(t)\}\{x\_i(0)\; /\; c\_i(0)\}\; =\; \backslash frac\{k\_s\}\{\backslash phi\_s\}\; \backslash frac\{m\_i(t)\}\{m\_i(0)\}$$

$A_i$ is the normalised RNA/DNA ratio. It is proportional to transcripts per cell, relative to before induction.

Removing the sample-level nuisance parameter

The factor $\frac{k_s}{\phi_s}$ is the RNA-to-DNA depth ratio of sample $s$ relative to $t = 0$. It is one unknown number per sample, and is shared by every guide. We get rid of it the same way as in the previous tutorial: take a reference.

For the non-targeting reference guides ($K_0 R = 0$, and $m_0$ constant),

$$A_0=\frac{k_s}{{\varphi }_s}A\_0\; =\; \backslash frac\{k\_s\}\{\backslash phi\_s\}$$

Before induction, dCas9 is off, so $R = 0$ and every guide is at steady state with $m_i(0) = \frac{E_{max}}{d + r_i}$. If the guides do not change growth before induction ($r_i = r_0$ at $t = 0$), then $m_i(0) = m_0(0)$ and the ratio of the two normalised ratios is the ratio of late steady states:

$$\frac{A_i}{A_0}=\frac{m_i}{m_0}=\frac{d+r_0}{d+r_i}\cdot \frac{1}{1+K_{i}R}\backslash frac\{A\_i\}\{A\_0\}\; =\; \backslash frac\{m\_i\}\{m\_0\}\; =\; \backslash frac\{d\; +\; r\_0\}\{d\; +\; r\_i\}\; \backslash cdot\; \backslash frac\{1\}\{1\; +\; K\_i\; R\}$$

Plotted against the reference term across late samples, each guide gives a straight line with gradient 1 and intercept $-\log(1 + K_i R)$. The ingredients are:

• Dividing RNA by DNA removes cell number and population growth.
• Dividing by each guide's own $t = 0$ removes per-guide capture efficiency and mRNA stability.
• Dividing by the pooled non-targeting guides removes the per-sample depth factor $\frac{k_s}{\phi_s}$.
• Late timepoints are on the transcript plateau, and the growth term is 1.

What remains is $-\log(1 + K_i R)$. Averaging the left side minus the reference term over the late samples gives this value per guide, $\hat{b}_i$. The most negative $\hat{b}_i$ is the strongest guide, so $-\hat{b}i = \log(1 + K_i R)$ ranks the guides. The transient after induction, the mRNA decay rate $d$, and $E{max}$ are not needed.