虚拟细胞的前夜：扰动测序与转录调控通路的因果机制推断（一）

1．扰动测序，复合扰动效应与基因互作推断

2016年，Aviv Regev课题组在Cell上发表了Perturb-seq扰动测序技术¹。基于腺病毒转染sgRNA所介导的CRISPR-Cas9敲除，借助合理控制感染复数（MOI）以及额外的荧光/抗生素筛选、测序后UMI统计，足以确保每个进入分析的单细胞转录组都只受到了单个敲除，或者说，基因扰动的影响。在这样的前提下，作者得以对细胞转录组变化受扰动的影响进行如下极为简单的线性建模 $Y=X\beta$$Y=X\beta$

其中，Y是观察到的细胞表达矩阵（的对数变换），X是实验的设计矩阵，对于每一个细胞对应的行，不仅包含其受到何种sgRNA扰动的信息，也包含有其他参量的信息（例如测序深度、细胞周期/其他细胞状态、实验批次等等）。这是因为细胞表达矩阵不止受到扰动的影响，也受到其他参量（尤其是测序深度）的影响。对于这些协变量的建模，能够分离它们对于基因表达的影响因素，增强模型的准确性和可解释性。最终，这个简单的线性模型通过ElasticNet回归的方法，得以求解调控矩阵β，从而得到每个设计因素（例如每个sgRNA或者其他技术上的协变量）对转录组的影响。

β是很可能没有解析解的。固然可以列出如下解析形式： ${{(X}^{T}X)}^{-1}{X}^{T}Y={{(X}^{T}X)}^{-1}{(X}^{T}X)\beta$${{(X}^TX)}^{-1}X^TY\ =\ {{(X}^TX)}^{-1}{(X}^TX)\beta$，但是，我们并不能确保$X^{T}X$$X^TX$真的可逆；即便真的可逆，由于这是一个巨大的稀疏矩阵（尤其是sgRNA设计部分的高度稀疏性)，意图求逆在操作上也是不现实的。

这个线性模型无疑十分朴素，但也足够强大，尤其是在单细胞测序刚刚出现的时代，过往基于测序数据对基因调控网络的分析还停留在WGCNA这种分析高表达相关性模块的思路，这篇文章就已经通过结合基因扰动测序和线性分解，从单细胞转录组数据，对转录因子扰动所产生的转录效应进行聚类，得到了转录因子的调控模块，以及这些模块所调控的基因表达程序，进而得到了转录因子模块和基因表达程序的有向调控网络，这简直可说是划时代的技术。然而，这一构建出的有向调控网络的可信性验证，将是一个很大的问题。在验证的TF的规模有限——例如文章中的情况时，尚可以通过生物实验（例如CHIP-seq数据）来验证特定语义下这一调控网络是否成立；但是对于更大规模的数据，我们固然可以构建出这一网络，但它必定是杂乱、多义、缺乏可靠性的。

Perturb-seq技术发表的同时，与Regev课题组合作的Jonathan Weissman课题组也发表了类似的大规模扰动测序研究²，与前述研究不同的是，他们使用了CRISPRi的技术路线——这一技术路线可能有助于降低CRISPR敲除所致的DNA损伤对转录组的整体影响，而代价则是设计矩阵中sgRNA的有无不再直观地带来对应基因表达的“全或无”，因为CRISPRi本身往往并不会导致基因表达的完全消失。六年后，他们沿用了这一技术并将其进一步地延拓至整个人类细胞基因组上，获得了K562细胞中全表达基因（接近10000个）、以及K562与RPE1细胞中超过2000个必需基因的扰动图谱³。

对于扰动后单细胞转录组的分析，这篇文章并没有更多新奇之处：扰动与扰动聚类得到调控模块，转录组和转录组聚类得到表达程序。对于每一个调控模块或表达程序，我们可以根据既有数据库对其进行本体论注释，使其映射到特定的生物结构、过程或者功能上；对于每一对调控模块与表达程序，我们都可以得到一个平均调控强度：意味着这个调控模块对于这个表达程序是起到了增强的作用，抑或是抑制的作用。相比于Regev组的文章，这篇6年后的文章不仅没有提出更新颖的算法，甚至不再执着于构建那个转录因子到转录因子的因果调控网络了，而其背后的原因，或许正如笔者在前文所述：当数据规模不再限于数十个基因的扰动，而是扩大到全部表达基因时，这一构建出的因果调控网络究竟有多可靠，将会是高度存疑的。

但是从后来的视角看，这篇文章中有一些有趣的结果：首先是在同一种细胞中，尽管可以看到，对同一种蛋白复合体中的基因的扰动所造成的转录组变化呈现出高度的一致性，但其实任意两种扰动之间的相关性，也体现出向正相关方向的偏倚，而并非以零点为中心。这意味着有可能扰动将会造成转录组向特定方向的整体偏移——例如扰动本身对于细胞状态的不利影响。其次，同一种扰动在同一种细胞中的转录组效应是高度相关的，但是在不同类型的细胞之间，其转录组效应的相关性则低得多。这也提示扰动的转录组效应本身可能是语义依赖（context-dependent）的。假如后来虚拟细胞技术的开发者能够早些注意到这些结果，那么对于相关工具的效力检验的设计，也许就会审慎得多。当然，这些是无法强求的。

至此，单基因扰动的分析流程几乎可说是发展到了极致。接下来需要做的不过是将这一技术应用于各种各样的样本和生物学过程而已。然而仍有一些问题并没有得到回答。第一，这些实验几乎都是单基因扰动的结果，而在实际问题中，扰动往往涉及到多个基因，并且即便其中各个基因的单个扰动效应已经被解明，复合扰动的结果也不完全是——尽管在多数情况下确实近似成立——其中各个基因效应的线性加和。我们可以将这种复合扰动给出如下粗略建模 $E=A{\beta }_A+B{\beta }_B+A\odot B{\beta }_{AB}$$E=A\beta_A+B\beta_B+A\odot B\beta_{AB}$，其中A代表扰动基因A的设计矩阵列，B代表扰动基因B的设计矩阵列， $A\odot B$$A \odot B$则为两列的哈达玛乘积。视${\beta }_A$$\beta_A$、${\beta }_B$$\beta_B$与${\beta }_{AB}$$\beta_{AB}$的相对关系，双基因复合扰动的效应可以粗略分为缓冲、协同、主导、相加四种。缓冲意味着两种基因扰动的效应相互拮抗，最终复合效应小于二者中任一者；协同意味着两者扰动的效应相互促进，复合效应至少大于其中一者；主导意味着其中一个扰动的效应完全掩盖另一者的效应；相加意味着两者互不干预（ ${\beta }_{AB}\approx 0$$\beta_{AB}\approx0$），复合效应约等于两者效应的简单加和⁴。

在Jonathan Weissman课题组的系列文章中，双基因的扰动互作给出了更加明确的定义。首先，他们于2018年的cell文章中⁵，基于双基因合成致死的概念，提出了GI score (Genetic Interaction score)。这一标量分数被定义为： $GI=\frac{{\rho }_{abobserved}-{\rho }_{abexpected}}{\sigma }$$GI=\ \frac{\rho_{ab\ observed}-\ \rho_{ab\ expected}}{\sigma}$，其中ρ被称作适应度，对于单个基因， ${\rho }_a={log}_2(\frac{{Reads}_{a,T_{end}}}{{Reads}_{a,T_0}})$$\rho_a=\ {log}_2(\frac{{Reads}_{a,\ T_{end}}}{{Reads}_{a,\ T_0}})$，其中${Reads}_{a,T_0}$${Reads}_{a,\ T_0}$和${Reads}_{a,T_{end}}$${Reads}_{a,\ T_{end}}$分别是测定文库中对于特定基因a扰动的sgRNA在导入时和给定时间段后的测序丰度，反映了该扰动对于细胞的致死性。考虑到扰动致死表型本身所具有的饱和限度，在基因a扰动上叠加基因b扰动后，适应度的0假设（简单加性效应）期待拟合为一个二次多项式关系${fit}_Q$${fit}_Q$，故：${{\rho }_{abexpected}=fit}_Q({\rho }_a)$${\rho_{ab\ expected}=\ fit}_Q(\rho_a)$。

另一方面，σ则为预期适应度窗口内，基因扰动$a$$a$+非靶向对照扰动$ntc$$ntc$的预期值与实际观测值的残差的分布离散度，作为对于实验噪声的建模。最终，GI的正负号决定了双基因扰动的互作是协同性的（GI<0），抑或是缓冲性的（GI>0）。

在这一标量GI的基础上，在Weissman课题组2019年发表于Science的文章中⁶，进一步地考虑了单细胞测序的表达谱向量水平上的双基因扰动互作定义。在这篇文章中，作者首先界定了112个与细胞生长相关的基因，然后介导了这些基因中所有两两配对或单基因组合（通过在所有设计的sgRNA中进行两次不放回采样来实现）的CRISPR激活。因此在指定细胞中和指定基因范围内，对任意基因扰动a，b，我们能够得到其单个存在同时存在时，表达组相比于基线水平产生的变化${\delta }_a$$\delta_a$、${\delta }_b$$\delta_b$与${\delta }_{ab}$$\delta_{ab}$。故而我们得以拟合线性组合： ${\delta }_{abexpected}=C_a{\delta }_a+C_b{\delta }_b={\delta }_{abobserved}-\epsilon$$\delta_{ab\ expected}=\ C_a\delta_a+\ C_b\delta_b=\ \delta_{ab\ observed}-\ \varepsilon$

其中$C_a$$C_a$、$C_b$$C_b$即为两种基因的扰动效应在另一者存在时的放缩尺度。结合之前的双扰动复合效应的朴素建模结果，不难得到 ${\beta }_{AB}=(C_a-1){\delta }_a+(C_b-1){\delta }_b$$\beta_{AB}=(C_a-1)\delta_a+(C_b-1)\delta_b$。在高维空间上计算${\delta }_a$$\delta_a$、${\delta }_b$$\delta_b$的共线性，结合系数$C_a$$C_a$与$C_b$$C_b$的取值，于是得以确定两个基因的关联模式：共线性足够高，即意味着两者为协同作用，若共线性低，则检查线性拟合系数：若两者比率接近1，则说明存在缓冲效应；若两者比率偏离1，则意味着其中一个基因具有主导作用（即上位基因)。

与此同时，我们注意到该线性拟合相比实际观测值的残差$\epsilon$$\varepsilon$：这项残差解释了基因互作的非线性效应，在这篇文章中被称作“新颖表型”。从文章中来看，绝大多数的扰动对的拟合$R^2$$R^2$系数是较高的，这也就意味着仅有一小部分双扰动组合存在非线性效应。这对于我们实际的建模来说，可能会是一个提示思路，即应当尽可能地对非线性机制施加稀疏约束。在文中，通过衡量${\delta }_{abexpected}$$\delta_{ab\ expected}$与${\delta }_{abobserved}$$\delta_{ab\ observed}$的方向相关性d（介乎于0~1之间)，也可以得到对于非线性效应强度的度量。

基于这样的扰动互作定义，即可通过分析基因的互作类型和上下位关系，来确定其基因调控通路上的顺序关系，文中的例子是DUSP9-MAPK1-ETS2的调控通路，展示了通过双扰动分析得到上下位关系和互作效应类别，最终所建立的线性顺序与生物学先验是一致的。然而，文章中没有、事实上也难以依据双扰动实验的数据，来自动化地建立一系列基因调控线性顺序（例如以具有边权重的有向无环图来体现），即便能够建立，这些顺序结果在多大程度上是可靠的、以及如何对其进行验证，将是较大的挑战。

此外，基因两两组合扰动实验的工作量，随着靶向基因的规模扩大将呈现平方律增加。在本研究中，单是研究112个生长相关的待测基因的任意两两组合，作者就测试了28680对独特的sgRNA对。考虑到对于任何一对sgRNA，都需要足够多的细胞、以确保测序结果对全转录组的充分覆盖、并且排除掉细胞状态波动对于扰动结果的影响，再扩大研究基因的范围所带来的工作量和成本将十分巨大——这还没有考虑扰动效应对于细胞种类和细胞状态是高度依赖的。当然从另一方面讲，这也使得现有的双基因扰动数据集格外有价值。但是，如果我们将高通量实验穷举的空间增加到3基因甚至更多的扰动组合，实验方法就变得不可行了。

我们回到第二个问题。单细胞测序技术本身作为一种终点性测量，并不能同时得到单细胞在扰动前和扰动后的转录组状态。在因果推断理论中，这被称为“反事实”配对数据。前述研究集中于扰动后状态的研究，仅仅将扰动前细胞的转录组作为衡量扰动效果的“基底”，即 ${\delta }_a=(T_a-T_0)=f\left(A\right)+\epsilon$$\delta_a=\left(T_a-T_0\right)=f\left(A\right)+\varepsilon$ （单基因扰动情况下，A为扰动是否存在/扰动的剂量，$\epsilon$$\varepsilon$为测量和状态噪声因素）。然而正如之前提到，细胞对于扰动的响应是状态依赖的，即是说更可能的情况是 ${\delta }_a=f(A,T_0)+\epsilon$$\delta_a=\ f\left(A,\ T_0\right)+\varepsilon$，其中$T_0$$T_0$作为细胞初始状态的一个表征。在这种情况下，我们该如何在缺少反事实配对数据的情况下，对特定细胞状态下的${\delta }_a$$\delta_a$进行更加合理的推断？更加激进的推演是，是否存在一个具有高度解释性的映射$f$$f$，使得 ${\delta }_{ab}=f(A+B,T_0)+\epsilon$$\delta_{ab}=\ f\left(A+B,\ T_0\right)+\varepsilon$？在拟合这个映射的参数同时，我们也就得到了A和B的互作模式？至此，基于变分推断原理的深度生成模型进入了人们的视野。