过去我们在“看”细胞，这一次研究人员开始逐个“动”基因
单细胞RNA测序（single-cell RNA sequencing，scRNA-seq）已经让我们能够观察一个组织中成千上万种细胞状态，但“看到两个基因同时变化”并不意味着知道谁在调控谁。
这正是基因调控网络（gene regulatory network，GRN）研究长期面对的问题：相关性容易观察，因果关系却很难建立。
这项研究采用的Perturb-seq，可以理解为把基因扰动与单细胞转录组测量结合起来。研究人员首先利用CRISPR干扰（CRISPR interference，CRISPRi）降低特定基因的表达，再通过单细胞RNA测序观察：当这个基因被压低之后，整个细胞的转录组发生了什么变化。
与完全敲除不同，CRISPRi更接近一种“降低表达量”的功能扰动。关键在于，它可以被大规模并行化：不同细胞携带针对不同基因的引导RNA（guide RNA，gRNA），最后再根据每个细胞携带的gRNA，把“扰动了谁”和“转录组发生了什么变化”对应起来。
研究建立的文库覆盖12,779个基因。更重要的是，研究没有只观察一种细胞状态，而是设计了三个时间窗口：静息状态（Rest）、再次刺激后8小时（Stim8hr）以及刺激后48小时（Stim48hr）。也就是说，研究人员不仅想知道“一个基因有什么功能”，还想知道：这个功能会不会随着T细胞所处的状态发生改变？
3340万个细胞过筛，最终留下近2200万个“可解释的扰动”
这种研究最容易被忽略的问题，是规模大并不等于数据可靠。
经过质量控制后，研究共获得3340万个细胞。其中，2199万个细胞，占65.8%，能够明确分配到单个靶向或非靶向gRNA；469万个细胞没有可靠检测到gRNA，另外670万个细胞同时被分配到多个gRNA，后者中相当一部分可能属于双细胞（doublets），因此没有进入主要分析。
在26,504条gRNA中，23,573条，也就是88.9%，至少在50个细胞中被检测到。平均而言，在每一种实验条件下，每个被扰动基因都有大约575个细胞提供信息。
随后研究人员检查CRISPRi到底有没有把目标基因压下去。在24,174条能够检测的gRNA中，73%产生了显著的靶基因下调，中位表达变化为log fold change=-2.33。
经过这些质量控制，研究最终能够估计11,527个基因被扰动后的全转录组效应。其中7,807个基因，占67%，在至少一种状态下显著影响了3个或更多下游基因。把这些关系全部展开，一共得到：
|
2,035,311条 显著的“调控基因—下游基因”反式调控关系（trans-regulatory effects） |
这已经不再是一张简单的“基因列表”，而是一个可以被系统分析的功能网络。
最重要的发现之一：基因的功能没有脱离“场景”独立存在
如果只记住这项研究中的一个观点，或许应该是：
|
基因调控具有很强的情境依赖性（context specificity）。 |
研究人员从1,860个调控基因产生的3,341个强扰动状态中，识别出了111个功能调控模块。这些模块涉及细胞信号、胞内运输、代谢、基因组稳定性以及转录调控等多个过程。
但值得注意的是，大约三分之一的调控模块只在三个实验状态中的一个或两个状态表现出较强的协同效应。这种差异并不能简单用基因表达量或者细胞数量解释。
甚至同一组调控因子，在三个状态下都可能保持活跃，但它们控制的下游基因却发生改变。
T细胞激活就是很好的例子。
刺激8小时后发挥作用的早期调控因子中，包括ATP2A2、ATP2A3等与钙稳态（calcium homeostasis）相关的基因；到了48小时之后，ERCC4、RAD51B、RNF8等与细胞周期和基因组稳定性有关的基因开始变得更加突出。
两者控制的下游程序也不同：8小时阶段更多影响IL2RA、LAG3等早期激活标志物，而48小时阶段则明显牵涉氧化磷酸化（oxidative phosphorylation）等代谢过程。
研究人员进一步挑选6个预测的早期或晚期调控因子进行独立实验，其中5个的时间特异性效应获得验证。
所以，“这个基因是做什么的？”也许本身就是一个不完整的问题。
更准确的问题可能应该是：
|
它在什么细胞、什么状态、什么时间点、什么外界刺激下做什么？ |
一个细胞因子背后，可能站着上千个调控基因
CD4+ T细胞发挥功能的重要方式之一，是分泌细胞因子（cytokine）。传统实验往往围绕IL-2、IFN-γ或者某一个特定细胞因子逐一寻找上游调控基因。
而这项研究一次分析了30种经典细胞因子。
结果发现，有1,556个基因扰动能够在至少一种状态下显著影响至少一种细胞因子的表达，统计阈值达到假发现率（false discovery rate，FDR）<1%。这些调控基因不仅来自熟悉的T细胞受体（T cell receptor，TCR）信号通路，还涉及Mediator、SAGA转录共激活复合物，以及线粒体功能、氨基酸感知、RNA加工和糖酵解等过程。
其中，IL10和IL21尤其值得关注。
研究发现，NFKB2和KDM1A竟表现出方向相反的双重作用：它们是IL10的正向调控因子，却同时是IL21的负向调控因子。 CYB5R4则被预测为IL21的正向调控因子；ATP2A2、ORAI1以及ELP2、ELP3等则表现为IL21负向调控因子。
为了判断这些结果是不是大规模筛选造成的统计假象，研究人员又独立验证了9个已知或候选调控因子，每个靶点使用两条独立gRNA，并同时进行了RNA测序和细胞内蛋白检测。
结果是：12组被检验的调控关系在方向上全部与Perturb-seq一致，其中10组达到统计学显著。
但这里也出现了一个很有意思的“不一致”：NFKB2扰动导致的IL10转录变化，并没有在蛋白水平完全重现。
这提醒我们，mRNA变化并不必然等价于最终蛋白功能变化。转录后调控（post-transcriptional regulation）依然可能重新改变结果。
能不能从一张转录组“照片”，反推出是谁推动了细胞状态？
接下来，研究走出了单纯做基因功能筛选的一步。
人群单细胞研究已经积累了海量“观察性数据”：我们知道Th1细胞和Th2细胞表达哪些不同基因，也知道老年人与年轻人的T细胞转录组有什么差别。
问题是：谁造成了这些差别？
研究人员尝试把一个细胞状态的表达特征，表示成多个基因扰动效应的组合。换句话说，如果一个自然存在的Th2细胞呈现某种转录特征，他们就寻找：哪些基因扰动的组合最能够“重建”这种状态？
在Th1/Th2极化（Th1/Th2 polarization）分析中，这套模型在发现队列中的平均交叉验证相关系数达到R=0.39，在独立验证队列中仍达到R=0.27，而且明显优于使用K562细胞系扰动数据建立的模型。
模型找回了许多熟悉的调控因子：IFNGR1、IFNGR2和JAK2指向Th1程序，IL4R、STAT6和GATA3则指向Th2程序。
更有价值的是，它还提出了一批此前缺乏明确功能证据的候选调控因子。进一步验证的13个高质量候选调控因子中，12个都让Th1/Th2特征基因向预测方向发生了显著移动，其中包括4个此前功能不明确的Th1调控因子和2个Th2调控因子。
这说明Perturb-seq的价值不只是发现某个基因影响某个分子，它还可能成为连接实验性因果扰动与人群细胞图谱的一座桥梁。
4名供者的扰动实验，为什么能解释数百人的T细胞衰老特征？
研究中另一个颇具启发性的分析来自衰老。
研究人员利用OneK1K人群队列中的CD4+ T细胞数据分析年龄相关变化，总共纳入782名供者，并另外保留199名供者进行验证。
随后，他们尝试用前面仅来自4名供者的Perturb-seq实验解释这些人群中的年龄相关表达变化。
结果显示，静息状态T细胞的扰动数据对年龄相关转录特征具有预测能力，交叉验证R达到0.366，并且优于刺激状态下的T细胞以及K562细胞数据。
模型首先指向了一个并不意外的通路——mTORC1信号通路（mechanistic target of rapamycin complex 1），包括TSC1、TSC2、DEPDC5、NPRL3和RPTOR；与此同时，蛋白SUMO化（SUMOylation）相关酶也被识别出来。
但数据并没有给出一个简单的“mTOR越高，衰老越明显”的故事。
部分mTORC1调控基因预测出的方向，反而与直觉预期相反。论文对此非常谨慎，提出负反馈、代偿机制以及高龄人群中的生存者偏倚（survivor bias）等多种可能解释，而没有把相关结果直接解释为简单的因果方向。
|
这个细节非常重要。大规模功能组学最有价值的地方，有时并不是给出一个简单结论，而是迫使我们意识到：一个通路的长期状态，与短期敲低其中某个基因产生的响应，并不是同一件事。 |
从GWAS的“相关位点”，走向疾病发生的功能通路
基因组关联研究（genome-wide association study，GWAS）已经找到大量自身免疫病相关位点，但一个风险位点并不会直接告诉我们：它在哪种细胞中发挥作用？影响什么基因？最终改变哪条免疫程序？
研究人员因此进一步将Perturb-seq与人群遗传学数据结合。
在UK Biobank淋巴细胞计数分析中，他们利用罕见编码功能缺失变异（loss-of-function variant，LoF）的效应与Perturb-seq中的基因调控效应进行比较。有趣的是，刺激后的原代CD4+ T细胞能够观察到明显的遗传信号富集，而静息T细胞、K562细胞等并没有表现出相同模式。
到了自身免疫病分析，研究人员整理了14种自身免疫疾病相关基因。在77个可检测的调控模块中，33个模块的调控因子或状态特异性下游基因出现显著的自身免疫病基因富集；这种富集主要出现在刺激后的T细胞模块，而不是只在静息状态起作用的模块。
研究还进一步提出LRRC25和FAM20B可能参与T细胞激活和Th1极化。特别是FAM20B，其扰动产生的全转录组效应在两条独立gRNA之间相关系数达到R=0.81，跨供者平均相关达到R=0.50。
需要强调的是，这些数据更适合被理解为功能候选机制的优先级证据，而不是已经证明这些基因就是疾病的直接致病靶点。
这一区别非常关键。
这项研究最大的价值，也许不是“发现了多少新基因”
2200万个细胞、200多万条调控关系、上千个细胞因子调控因子，这些数字当然足够醒目。
但如果只关注数量，反而容易错过这项工作的核心。
过去二十多年，人类遗传学解决了大量“哪些遗传变异与疾病有关”的问题；单细胞组学则告诉我们“疾病状态下哪些细胞、哪些基因发生了变化”。
两类数据都很丰富，却缺少中间的一层：
|
如果改变一个基因，细胞究竟会发生什么？ |
Perturb-seq正在尝试补上这一层，从而建立一条更完整的逻辑链：
|
遗传变异 → 调控基因 → 基因程序 → 细胞状态 → 人体性状与疾病 |
这也是该研究提出其数据可能服务于AI虚拟细胞（AI virtual cell）模型的重要原因。如果未来的模型希望预测“扰动一个基因后细胞会怎样变化”，它需要的就不只是海量观察数据，而是大量经过实验验证的扰动-响应关系。
但这项研究也远没有把问题解决完。
实验只有4名供者；每个基因主要使用2条gRNA；CRISPRi只是部分抑制而不是完全敲除；基于探针的单细胞RNA测序无法完整覆盖某些非编码RNA、剪接异构体和UTR层面的调控；伪批量分析（pseudobulk analysis）提高了统计稳定性，却可能掩盖细胞间的响应异质性。更重要的是，这些T细胞来自外周血，并在体外培养，真实组织微环境中的调控网络很可能更加复杂。
因此，这项工作的意义并不是宣布我们已经获得了一张完整的“T细胞线路图”，而是提供了一套可以继续扩展的实验框架。
而它留下的一个问题，可能比任何单个新发现都更值得思考：
|
当一个基因在静息、激活8小时和激活48小时后的作用都可能不同，我们还能不能脱离细胞类型、时间和微环境，孤立地讨论“这个基因的功能”？ |
也许未来功能基因组学（functional genomics）的基本单位，不再只是“基因”，而是：
基因 × 细胞类型 × 状态 × 时间
这可能也是理解复杂疾病遗传机制时，需要重新建立的一套坐标系。
参考文献
