Chestnut Studying
摘要
Macrophages induce the expression of hundreds of genes in response to immune threats. However, current technology limits our ability to capture single-cell inducible gene expression dynamics. Here, we generated high-resolution time series single-cell RNA sequencing (scRNA-seq) data from mouse macrophages responding to six stimuli, and imputed ensembles of real-time single-cell gene expression trajectories (scGETs). We found that dynamic information contained in scGETs substantially contributes to macrophage stimulus-response specificity (SRS). Dynamic information also identified correlations between immune response genes, indicating biological coordination. Furthermore, we showed that the microenvironmental context of polarizing cytokines profoundly affects scGETs, such that measuring response dynamics offered clearer discrimination of the polarization state of individual macrophage cells than single time-point measurements. Our findings highlight the important contribution of dynamic information contained in the single-cell expression responses of immune genes in characterizing the SRS and functional states of macrophages.
巨噬细胞在免疫威胁下会诱导数百个基因的表达。然而，目前的技术限制了我们捕捉单细胞诱导基因表达动态的能力。在这里，我们从小鼠巨噬细胞中获得了高分辨率的时间序列单细胞RNA测序（scRNA-seq）数据，这些数据对六种刺激做出了反应，并推断出实时单细胞基因表达轨迹（scGETs）的集合。我们发现，scGETs中包含的动态信息极大地促进了巨噬细胞刺激反应特异性（SRS）。动态信息还揭示了免疫反应基因之间的相关性，表明了生物协调性。此外，我们还发现，极化细胞因子的微环境对scGETs有深刻影响，因此，与单时间点测量相比，测量反应动态能够更清晰地区分单个巨噬细胞的极化状态。我们的研究结果强调了免疫基因的单细胞表达反应中包含的动态信息在描述巨噬细胞的SRS和功能状态方面的重要贡献。
实验结果1
模拟scGET的动力学
通过活细胞显微镜对巨噬细胞进行单细胞信号传导研究，揭示了转录因子或激酶活性的可诱导动态变化。这些信号传导蛋白轨迹的特征取决于刺激，但也表现出细胞间的异质性（图1A）。然而，在每个细胞内，信号传导蛋白的动态活性会导致数百个独特基因的进一步下游激活。为了探索单个细胞中基因表达的动态变化，作者根据先前的文献构建并模拟了免疫反应基因的数学模型，每个基因的合成和降解速率由基因特异性参数和先前测量的脂多糖（LPS）刺激后信号效应分子活性来调节。基因被分类 五个先前已描述的基因反应模块（由激活蛋白-1 [AP1]、NF-κB、干扰素调节因子[IRF]、NF-κB|p38或NF-κB|IRF调节的基因），其中每个基因反应模块中的基因在激活和降解控制效应物方面有所不同（图1B）。这六种刺激中的每一种都以特定于刺激的动态方式激活信号传导效应。模拟结果显示，在LPS刺激下，每个细胞中200个独特基因的诱导动态完全相同，而细胞间的异质性则源于单细胞信号传导的异质性（图1C）。
为了量化基因表达的动态模式，作者考虑了一系列表征基因表达轨迹的指标：相对峰值振幅（Peak Amp）、从时间零开始的最高对数变化倍数（Max LFC）、时间序列的总表达量（Integral）和激活速度（Speed）（图1D）。通过计算模拟的真实基因表达轨迹的这些特征，作者用一组可解释的值总结了轨迹的单细胞异质性（图1E）。值得注意的是，由于每个单元格包含数百个基因的轨迹，因此可以确定单个细胞中基因与基因之间的相关性。正如预期的那样，在单个细胞中，受具有共同转录因子的基因调节因子（GRM）调节的基因具有高度相关的轨迹特征，而受不同转录因子调节的基因的相关性较低（图1F）。
实验结果2
一种归因方法可以从时间序列数据中识别出scGETs
为了确定scGETs是否可以从scRNA-seq时间序列数据中推算出来，作者开发了一种方法，该方法包括两个不同的步骤：相邻时间点之间测量细胞的直接连接以及未测量时间框架之间的插值（图2A）。与其他scRNA-seq分析方法类似，作者对所有时间点的数据采用了主成分分析（PCA），目的是在恢复数百个细胞的单个基因轨迹之前，先组装单个细胞的轨迹（图2B）。为了改善单细胞测量的稀疏性，作者通过k均值聚类构建了细胞原型（元细胞）。细胞原型在时间点之间通过随机游走进行连接，随机游走的权重由转换概率矩阵决定，与相邻时间点之间细胞原型对的欧氏距离成反比（图2C）。加权随机游走使得细胞与细胞之间的连接排列成为可能，但降低了轨迹在不太可能路径上的概率（图2D）；转换概率矩阵的指数化使得某些细胞原型的包含更加确定，但降低了噪声的容忍度。为了在测量的时间点之间进行插值并获取单细胞轨迹的时间过程，作者随后在每个主成分上拟合了样条曲线。结果是一组分解为主成分的轨迹曲线，每个主成分的样条曲线描述了细胞转录组轨迹的特征（图2E）。然后，使用原始的PCA载荷，将细胞轨迹投影回原始基因表达空间，以恢复每个基因的数百个细胞的相关转录组尺度GET（即每个细胞数百个单个基因的实时表达动态）（图2F）。然后，通过分解动态特征来描述scGET（图2G）。
实验结果3
根据模型模拟的真实数据评估估算的表达轨迹
为了评估推导出的轨迹，作者逐个步骤将输出结果与模型模拟提供的真实动态进行比较。作者将按转移概率矩阵加权的随机游走与完全随机的单元间连接进行了比较，发现无论使用多少时间点，使用加权随机游走生成的路径都更接近真实路径（图2H）。作者还比较了选择不同时间点分布对单元间连接的影响，发现前向加权的五个时间点比均匀分布的五个时间点效果更好（图2H）。作者改变了用于计算细胞间联系的组件数量，发现增加组件数量并不能显著提高性能，但使用加权随机游走和五个前向时间点在任何组件数量下都能产生更接近真实值的路径。
在模拟的具有不同mRNA降解速率的基因中，估算的轨迹集合的轨迹特征与真实轨迹的特征是一致的（图2I）。作者探索了不同的插值选项，以平滑选定单细胞路径上的轨迹，因为这种选择可能会影响轨迹特征的量化，包括局部估计散点图平滑（LOESS）、具有三个（spline.df3）或五个（spline.df5）自由度的38次立方样条 、使用逐个剔除交叉验证（spline.cv）的立方样条函数以及单调Hermite样条函数（spline.mono）。作者发现，每种方法都能在一系列具有不同合成和降解特性的基因中合理地概括出地面真实值的不同轨迹特征。然而，在各种插值方法中，具有快速降解和高解离常数的基因轨迹的积分值被更准确地概括。作者继续通过逐个剔除交叉验证来拟合三次样条曲线，因为它们既提供了基因特异性特征，又保持了与特征的地面真实分布相比的单细胞异质性。
实验结果4
时间序列单细胞转录组测序揭示了巨噬细胞反应中随时间变化的异质性
根据对代入法的性能分析，作者在刺激（0、0.25、1、3和8小时）后，从五个相对较早的时间点（0、0.25、1、3和8小时）的数千个巨噬细胞中生成了时间序列scRNA-seq数据。作者用一组六种免疫配体刺激了处于三种不同极化状态的巨噬细胞，这些免疫配体包括细胞因子或不同的细菌或病毒成分，它们分别与不同的Toll样受体（TLR）或细胞因子受体结合： LPS（TLR4）、poly(I:C)（PIC）（TLR3）、Pam3CSK4（P3C）（TLR2）、CpG（TLR9）、肿瘤坏死因子（TNF）（TNF受体[TNFR]）和干扰素（IFN）-β（干扰素-α/β受体[IFNAR]）（图3A）。这产生了超过72个scRNA-seq样本，其中包含超过10万个单个巨噬细胞。作者注意到，基因表达异质性在很大程度上取决于时间点，因此没有一个时间点能够准确反映巨噬细胞对刺激的反应的异质性（图3B）。为了降低脱落率并提高细胞因子基因等目标基因的测序深度（图3B），作者对一组500个免疫反应基因进行了测序，这些基因先前已被确定为在巨噬细胞中具有高SRS12。
作者从M0巨噬细胞（0、0.5、3、5和8小时）的重复时间和互补时间点收集了重复数据，以评估数据质量和可重复性。重复时间点产生了一致的数据。单细胞数据的伪块模式也与先前发表的大规模巨噬细胞刺激反应RNA-seq数据一致。
实验结果5
免疫威胁下单个巨噬细胞中scGET的归因
由于时间序列scRNA-seq数据是快照，它们无法捕捉单个细胞免疫反应转录组的一个基本方面——其可诱导动态。作者将伪时间算法和scGET插补方法应用于作者用LPS刺激的巨噬细胞的时间序列scRNA-seq数据。与伪时间轨迹不同，提供了单个细胞沿着有时分支轨迹的顺序，而scGETs则提供了异质基因诱导特征，如倍数变化、峰值幅度、诱导速度或mRNA丰度的累积总量。将几个基因的伪时间轨迹与单细胞实时动态进行比较，可以说明这种对比。在多个基因中，根据表达谱将单细胞按时间顺序排列，并且可以根据基因的伪时间模式进行聚类；在scGETs的插补中，由于数据是在时间序列中收集的，因此顺序是已知的，并且相邻时间点之间的单细胞原型是明确关联的，从而可以根据scGETs对单细胞进行聚类（图S4H）。
作者将scGET插补方法应用于M0巨噬细胞群在6种不同配体刺激下的时间序列scRNA-seq数据，从而恢复了刺激8小时内的连续单细胞表达动态，例如Tnf（图3C）。作者通过比较每个测量时间点scRNA-seq数据中推导轨迹集合的平均值和方差与测量平均值和方差的差异，评估了每个基因的推导单细胞轨迹的准确性。值得注意的是，一些基因的估算结果比其他基因更准确。Tnf是估算结果较差的基因之一，其均值和方差的归一化均方根偏差（RMSD）相对较高。然而，即使Tnf的轨迹仍然与数据相当匹配，并保持了预期的SRS（图3C）。其他基因的平均值和方差在测量时间点与数据更接近，例如Ccl5和Cxcl10。此外，scGET的诱导模式与每个基因的已知GRM一致，同时提供了表达动态中细胞间异质性的估计值。
重要的是，通过scGETs的估算，作者可以利用基因表达的整个轨迹信息对单个巨噬细胞进行聚类，而不仅仅是它们在特定时间点的基因表达。单个细胞因子基因Tnf的scGET显示，仅使用Tnf作为读数，巨噬细胞对至少三个刺激组（IFN-β、细菌PAMP（LPS、P3C、CpG）和TNF/PIC）的反应明显不同（图3D）。相比之下，仅在单一时间点（如刺激后3小时）测量的Tnf表达甚至无法区分两个刺激组，12 这支持了这样一个观点，即对免疫威胁的定制反应依赖于Tnf表达动态，而不仅仅是表达水平。
此外，scGETs的归因恢复了每个细胞中所有测量基因的共同调节表达动态，从而能够根据多基因时间进程的异质性评估刺激特异性反应（图3E）。通过对五种基因（促炎细胞因子Tnf和Il6、淋巴细胞趋化因子Ccl5和Cxcl10以及NF-κB负反馈调节因子Nfkbia）动态信息组合进行分层聚类，进一步区分了响应不同刺激的单个细胞，表明不同巨噬细胞反应的稳健性是通过多基因互补动态表达特征实现的。用TNF和PIC刺激的细胞，仅通过Tnf动态无法区分（图3E），现在可以通过仅对PIC的Cxcl10和Ccl5的晚期诱导来区分（图3E）。测量刺激后1小时等单个早期时间点无法捕捉基因表达的时间差异。绘制表达分布边界进一步表明，刺激对之间的特异性取决于特定基因在时间过程中的异质性调节，例如，在Nfkbia的晚期时间点，LPS与TNF反应分布部分重叠，但在Tnf和Cxcl10上则截然不同（图3F）。作者注意到，从约100个细胞的估算中可以看出分布边界是稳定的，这表明scGET的估算数量（1000个细胞）足以捕捉刺激特异性反应分布的特征（图3G）。
实验结果6
scGETs的动态特性传递刺激信息
为了量化scGET的动力学模式，作者再次考虑了基因表达反应的四个动力学特征：峰值振幅、最大LFC、积分和速度。对于Tnf、Nfkbia和Cxcl10，绘制轨迹特征的单细胞分布图表明，峰值振幅和积分比最大LFC或速度更具刺激特异性，单细胞分布重叠较少（图4A）。为了量化SRS，作者计算了六种刺激输入和表达输出之间的最大互信息（max MI）。max MI越高，表明刺激信息越多，基因的动态特征越具有刺激特异性。作者发现，与以动态特征形式提供的轨迹信息相比，每个基因在单个测量时间点的max MI始终较低（图4B）。此外，我们注意到，与单个时间点相比，使用四个随机链接的时间点可以略微提高最大MI，但通过加权随机游走链接的四个时间点对所有基因的SRS都有显著提高（图4B）。
为了研究scGET的特征对SRS的揭示程度，作者接下来评估了每个基因的每个特征的最大MI。作者发现Integral和Peak Amp的平均最大MI最高，而Max LFC和Speed的信息量较少（图4B）。有趣的是，某些基因更依赖于动态特征的差异来保留刺激特异性信息。例如，NF-κB靶基因Nfkbia在某个时间点的表达量远低于通过scGET插补方法关联的多个时间点或Nfkbia动态特征（图4B）。这一发现与之前关于Nfkbia诱导动力学重要性的研究结果一致。Nfkbia编码一种关键的NF-κB反馈蛋白，在免疫威胁中产生适当的反应。
作者注意到，对于某些基因子集，轨迹特征揭示的刺激信息远远多于单个时间点的表达信息，而其他基因的信息量则较少，例如Cxcl10（图4B）。为了进一步研究，作者根据先前文献中激活免疫反应基因的信号通路对它们进行了分类，分为五类（AP1、NF-κB、IRF、NF-κB“与”p38、NF-κB“或”IRF）。有趣的是， 由单个转录因子（如AP1、IRF或NF-κB）调控的基因，在考虑表达积分而非表达水平时，信息增益约为1比特，但由NF-κB|IRF“顺序或”门调控的基因，信息增益仅为0.5比特左右（图4C）。这种模式表明，在基因启动子处，两个转录因子的组合活性可以通过组合控制来区分刺激，而较少依赖基因表达的时间动态，而由单个转录因子诱导的基因则更有可能通过动态控制来显示SRS（图4C）。
实验结果7
轨迹特征揭示了时间点测量中不明显的基因相关性
在单细胞的转录组测量中，由于每个细胞内相似的染色质或信号环境，不同信号通路的靶基因可能存在关联。为了确定单细胞中动态依赖的基因-基因关联，作者比较了作为相同或不同转录因子靶标的基因对的关联性，无论是单时间点的scRNA-seq数据还是scGET特征。在1、3或8 h的任何单个响应时间点，两个突出的NF-κB调节基因Tnf和Nfkbia均显示弱相关或无相关（图5A）。单个时间点缺乏相关性可能归因于激活时间不同：Tnf表达更持久，而Nfkbia更前向，因此任何单个时间点都无法很好地捕获总表达。然而，考虑到轨迹特征，Integral显示基因与基因之间在各种刺激下以及每种刺激条件下均具有高度相关性（图5B）。检查这两种基因的单细胞Integral值时，还发现与Nfkbia相比，Tnf表达在CpG、P3C或LPS刺激下的细胞中有所增加（图5B）。CpG、P3C和LPS是免疫配体，通过Myd88发出信号激活MAPKp38，导致NF-κB|p38调节的Tnf表达相对于仅由NF-κB调节的Nfkbia表达增加。只有scGET功能可以清楚地揭示Tnf-Nfkbia相关性的微妙变化，这反映了这两个基因的已知调节。
然后，作者研究了编码Cxcr3配体的两个基因，Cxcl9和Cxcl10，它们都是由巨噬细胞产生的趋化因子，可导致T细胞募集（图5C）。在刺激后的单个时间点，Cxcl9和Cxcl10在单个细胞中的相关性并不高，只有Cxcl10显示由LPS、IFN-β或PIC诱导的IRF信号的早期激活（图5C）。然而，Cxcl9和Cxcl10的积分与时间点数据相比具有更大的相关性（图5D）。此外，积分检查表明，在响应相同刺激的每个细胞群中，Cxcl9和Cxcl10在单个细胞中呈正相关，这表明通过共享激活信号通路或相似的响应染色质环境来实现共同调节。从更广泛的角度来看，在3小时这一单一时间点，所有基因的相关性热图显示，受相同GRM调控的基因之间存在微弱的相关性，而根据Integral对基因进行聚类后，NF-κB和IRF调控的基因簇之间的相关性更强，从而更好地区分了它们（图5E）。当基因对Integrals受相同的GRM或具有共同转录因子的GRM调控时，它们之间表现出高度相关性，但一些控制不同GRM的基因对也表现出相关性（图5F）。对这些基因对的检查揭示了它们之间的相关性，例如NF-κB调节的细胞因子（如Il1a、Ccl5）与IRF调节的基因（如病毒诱导的Heatr9）之间的相关性，后者参与细胞因子的前馈调节。因此，单细胞基因诱导的完整动态轨迹对于捕捉基因对之间的许多其他相关关系、识别共享染色质或信号环境的影响至关重要。
实验结果8
极化改变了巨噬细胞特异性免疫反应基因的表达
当巨噬细胞被微环境细胞因子极化时，scGETs的动态特征分布会如何变化？接下来，作者使用来自数千个M1（IFN-γ）和M2（白介素[IL]-4）极化巨噬细胞的时间序列scRNA-seq数据来推断scGETs，这些巨噬细胞同样被六种免疫配体中的每种配体刺激。根据scGETs对细胞进行分层聚类，初步检查显示，多个免疫反应基因的平均GETs发生了预期的变化（图6A）。例如，在细菌PAMP（脂多糖、CpG和P3C）的作用下，NF-κB靶基因Nfkbia在M1（IFN-γ）巨噬细胞中表现出更短暂的活性，3小时后活性减弱，而在M0巨噬细胞中则表现出更持久的活性。NF-κB|IRF靶基因Ccl5和Cxcl10在M1(IFN-γ)巨噬细胞中诱导速度更快，而在M2(IL-4)巨噬细胞中诱导速度则慢得多（图6A）。M1(IFN-γ)和M2(IL-4)巨噬细胞的聚类分析表明，反应仍具有刺激特异性，但特异性可能归因于基因动态的不同方面。
为了识别在极化状态下变化最大的具有scGET动态特征的基因，作者使用来自所有极化状态的细胞对scGET动态特征进行了PCA分析。作者发现，在PC载荷中，某些GRM调节的基因的权重更高。对于最大LFC，NF-κB和NF-κB|p38调节的基因在PC1上的权重最大，而IRF和NF-κB|IRF调节的基因在PC2上的权重最大，其中Mx1和Ifit1（图6B）。绘制所有单细胞的Mx1和Ifit1最大LFC曲线图，可以看出，这些IRF基因的刺激反应分布差异在M1（IFN-γ）巨噬细胞中消失了，最明显的是在CpG和LPS两种细菌刺激之间（图6C）。IFN-γ分化巨噬细胞特异性的丧失是由于刺激之间的平均反应更相似（例如Mx1最大LFC）或反应的异质性增加（例如Ifit1最大LFC）。对速度进行的类似分析表明，M0巨噬细胞中由反应速度（例如Mx1或Ifit1速度）产生的特异性在M2（IL-4）巨噬细胞中消失，最明显的是在IFN-β和所有其他刺激之间（图6D）。综上所述，此处评估的微环境细胞因子（IFN-γ和IL-4）对轨迹特征的异质性有显著影响，例如Max LFC和IRF基因的速度。
由于极化深刻地影响了许多基因的表达动态，作者想知道这些变化对巨噬细胞SRS的影响有多大。作者确定了动态表达特征最能反映刺激特性的基因组合。考虑到M0巨噬细胞中scGETs的积分，仅用三个基因（Mx2、Nfkbiz、Egr3）（IRF基因、NF-κB |p38基因和AP1基因的组合），接近理论最大值2.58比特（区分六个刺激，如2 1 = 2、2 2 = 4、2 2.58 = 6）。在M1（IFN-γ）（Mx2、Tnf、Gna15）和M2（IL-4）巨噬细胞（Cited、Ehd1、Pou2f2）中，三个基因的组合也可以达到接近理论上限的最大MI（图6E）。相比之下，使用单个时间点（3小时）的前三个基因组合最多只能达到2比特的信息（图6E）。绘制所选基因的积分图证实了这样的结论：动态轨迹特征将六种刺激物明显区分开来，尽管由于CpG和P3C分布的重叠，一些刺激物信息仍然丢失（图6F）。因此，虽然Cxcl10等特定基因的SRS在极化巨噬细胞中可能会降低，但巨噬细胞在极化后仍保持SRS，只是依赖不同的基因组合来调节特异性。
实验结果9
scGETs能够有效区分极化状态
巨噬细胞的状态通常通过准稳态下的蛋白质标记物或单细胞转录组学来评估。然而，巨噬细胞的功能是在刺激下发挥的。因此，作者假设刺激反应基因表达动态可能提供更准确的细胞功能状态标记。为了研究稳态测量和scGET在区分极化状态方面的差异，作者重点关注了用于识别M1（IFN-γ）或M2（IL-4）极化的典型标记基因，包括Nos2（M1标记基因）、Cd86（M1标记基因）和Retnla（M2标记基因）。将它们在未受刺激时的稳定状态下的表达值与它们的LPS反应积分值进行比较，结果表明scGET积分值可以更清晰地区分三种巨噬细胞极化状态（图7A）。平均而言，Nos2在M1巨噬细胞中的表达水平更高，但Nos2的单细胞表达分布与极化状态重叠，因此不能仅根据一个时间点的数据就确定所有单细胞是否属于M1型（图7B）。根据这三种标记基因对极化状态进行分类，结果表明，与在单个时间点测量的稳态值或反应相比，轨迹特征（尤其是积分、峰值振幅或最大LFC）极大地改善了极化状态的识别（F1评分更高）（图7C）。
作者想知道稳态下的其他基因是否比经典标记基因更能区分极化状态。根据MI分析，最能区分未受刺激的M0、M1（IFN-γ）和M2（IL-4）巨噬细胞的三个基因分别是Irf1、Fgl2和Tgtp1。虽然稳态表达值无法区分M0和M2(IL-4)巨噬细胞，但同一三个基因的LPS反应积分可以完美区分三种极化状态，这表明编码TF（Irf1）或分泌蛋白（Fgl2）的基因的反应动态对于确定细胞功能状态也很重要。
为了研究解释极化状态差异所需的基因数量，作者接下来使用稳态或刺激反应数据构建了一个多变量最小绝对收缩和选择算子（LASSO）惩罚回归模型。作者发现，对稳态表达值（0小时）进行LASSO正则化回归后，模型中包含约150个基因，这些基因基于单时间点响应值（1、3和8小时）构建 ∼50-120个基因，而那些使用轨迹特征（如积分）构建的模型仅包含∼25-30个基因，这再次表明GET在指定单细胞功能状态时具有更高的信息含量。
