学者追踪 | Cell | Seurat作者新作，多重单细胞替代

[图片上传失败...(image-d9e4ae-1727280015647)]

Basic Information

• 英文标题： Multiplexed single-cell characterization of alternative polyadenylation regulators
• 中文标题：多重单细胞替代聚腺苷酸化调控因子的特征分析
• 发表日期：8 August 2024
• 文章类型：Resource
• 所属期刊：Cell
• 文章作者：Madeline H. Kowalski | Rahul Satija
• 文章链接：https://www.sciencedirect.com/science/article/pii/S0092867424006457

Highlights

Para_01

• CPA 调节因子的遗传扰动揭示了聚腺苷酸化位点使用的变化是协调一致的
• PASTA 可量化单细胞 RNA 测序数据中多聚腺苷酸化位点的使用异质性
• RNA生命周期的不同组成部分影响内含子多A位点的使用
• 深度学习识别决定扰动响应的序列特征

Summary

Para_01

1. 大多数哺乳动物基因具有多个多聚腺苷酸化位点，这些位点代表了由剪切和多聚腺苷酸化（CPA）机制调控的转录多样性的重要来源。
2. 为了更好地理解这些蛋白质如何控制多聚腺苷酸化位点的选择，我们介绍了CPA-Perturb-seq，这是一个包含42个CPA调节因子的多重扰动筛选数据集，它使用3'单细胞RNA测序读出，能够对全转录组的聚腺苷酸化位点使用进行推断。
3. 我们开发了一个框架，用于检测依赖于扰动的聚腺苷酸化变化，并描述了受共调节的多聚腺苷酸化位点的模块。
4. 我们发现了一组内含子多聚腺苷酸化位点，它们受到核RNA生命周期不同组成部分的调控，包括延伸、剪接、终止和监控。
5. 我们训练并验证了一个深度神经网络（APARENT-Perturb），用于串联多聚腺苷酸化位点的使用，勾画了一个顺式调控代码，该代码可预测扰动响应并揭示调控复合物之间的相互作用。
6. 我们的工作强调了多重单细胞扰动筛选进一步理解转录后调控的潜力。

Graphical abstract

[图片上传失败...(image-574d9f-1727280015646)]

Keywords

• cleavage and polyadenylation; alternative polyadenylation; post-transcriptional regulation; Perturb-seq; RNA processing

Introduction

Para_01

1. RNA切割和多腺苷酸化（CPA）是转录后调控机制，对于真核生物前mRNA的成熟是必需的。
2. 大多数哺乳动物基因含有多个多腺苷酸化位点，使得单个基因能够通过替代多腺苷酸化编码多个mRNA转录本。"," Sentence_03 ":" 由这一过程产生的不同的3'端可以影响RNA生命周期的多个不同阶段。
3. 例如，3'非翻译区（UTR）的缩短可能会影响转录本的稳定性和定位。
4. 而 内含子位点的替代多腺苷酸化可能导致截短的编码或非编码转录本的生成。
5. 更一般地说，多腺苷酸化的广泛变化已在许多生物学背景下得到证实，包括细胞增殖、肿瘤发生、胚胎发育和分泌细胞分化。
6. 生化和分子研究表明，一组核心和辅助蛋白质负责调控多腺苷酸化位点的选择。
7. 例如，CPA特异性因子（CPSF）复合物催化切割，切割因子Im（CFIm）和切割因子IIm（CFIIm）复合物结合辅助识别序列，而多腺苷酸化聚合酶负责添加多腺苷酸化尾。

Para_02

1. 全基因组3'端转录组技术可以用来分析多腺苷酸化位点使用情况的变化，包括CPA调节因子的遗传扰动。
2. 尽管有些研究进行了个别或小规模的扰动，其他研究则利用小干扰RNA筛选方法产生了更多的资源。
3. 这些研究已经表征了单个调节因子在全球范围内促进3'非翻译区内的近端或远端多腺苷酸化位点使用的倾向，并调节一种称为串联替代多腺苷酸化（图1A）的过程。
4. 尽管不同的扰动影响不同数量的位点，但目前尚不清楚这种变异是否反映功能性的共调节：即，是否有多腺苷酸化位点的组对不同的调节因子或亚复合物的扰动特别敏感。
5. 如果是这样，识别这些共调节位点的模块及其使用的分子特征代表了增进我们对CPA调节理解的关键目标。

[图片上传失败...(image-ff806b-1727280015646)]

• 图 1. CPA-Perturb-seq 概述（A）（顶部）用于生成 CPA-Perturb-seq 数据集的实验工作流程示意图。（底部）与串联或内含子多腺苷酸化相关的扰动依赖性变化的示意图。"," Sentence_02": "（B）图示组成并参与切割和多腺苷酸化机制的的核心调控复合物。"," Sentence_03": "（C）读覆盖图描绘了 CBX3 位点上替代多腺苷酸化位点的差异使用。每条轨迹代表一个伪批量平均细胞，按其扰动分组。下图显示了在检测到的多腺苷酸化位点之前的 ENSEMBL 基因模型和峰（量化区域）。"," Sentence_04": "（D）通过 CPA-Perturb-seq 对 HEK293FT 细胞进行 UMAP 可视化。细胞根据目标基因的同一性进行着色，使用与（C）中相同的颜色。该可视化基于转录组范围内多腺苷酸化位点计数的线性判别分析（LDA）计算得到。另见图 S1-S3。

Para_02

1. 替代性聚腺苷酸化也可以发生在内含子聚A位点（内含子替代性聚腺苷酸化；图1A），这与3' UTR的变化不同，会导致编码序列的改变。
2. 内含子聚A位点的使用与多种调控蛋白相关，这些蛋白控制RNA转录本的合成或核内处理。
3. 尽管在内含子聚腺苷酸化改变已在疾病状态下得到识别，但目前尚不清楚内含子聚A位点是否普遍对转录动力学变化敏感，或者是否特定的调控因子子集决定了不同内含子位点的使用。

Para_03

1. 多重单细胞技术，如Perturb-seq，利用单细胞RNA测序（scRNA-seq）对分子扰动的全转录组进行高通量表征，具有解决这些问题的激动人心的潜力。
2. 虽然scRNA-seq通常用于分析基因表达水平的异质性，但这些数据也可以用来描述转录结构的变化。
3. 大多数scRNA-seq协议旨在捕获聚腺苷酸化mRNA转录本的3'端。
4. 因此，这些方法非常适合在单细胞分辨率下量化全转录组的聚A位点使用情况以及基因丰度，揭示细胞分化和疾病过程中聚腺苷酸化的动态变化。

Para_04

1. 在这里，我们介绍CPA-Perturb-seq，一个资源库，其中我们在多路复用的3'端单细胞RNA测序屏幕中扰动已知的CPA调节因子，并在单细胞分辨率上量化每个扰动对多聚A位点使用的影响。
2. 我们引入了新的统计方法来量化稀疏单细胞数据集中多聚A位点使用的变化，并识别了共调节多聚A位点的不同模块。
3. 我们发现，内含子多聚A位点使用的共调节是由核RNA生命周期不同元素对扰动的差异敏感性驱动的，而对于串联位点，共调节是由顺式调控密码驱动的，其中单个序列元素调节反应性。
4. 我们通过将开创性的深度学习模型扩展到多个遗传背景来学习这个密码，这些模型用于替代性多腺苷酸化，并通过使用大规模平行报告子分析（MPRA）来验证我们的发现。
5. 最后，我们展示了我们的计算工具如何应用于任何3'端单细胞RNA测序数据集，并使用基因组规模的Perturb-seq资源表征了参与RNA处理的数百个基因的调控效果。

Results

Multiplexed Perturb-seq screens of 3′ polyA site usage

多重 Perturb-seq 筛选 3' 多聚腺苷酸位点使用情况

Para_01

1. 我们试图了解基因在CPA中的系统扰动如何在单细胞分辨率上影响选择性多腺苷酸化（图1A）。
2. 我们设计了一个包含162个单导向RNA（sgRNA）的库，针对42个基因和10个非目标（NT）对照（表S1）。
3. 我们的目标集合包括18个已知的核心CPA复合物成员基因，包括CFIm、CFIIm、CPSF和裂解刺激因子（CSTF）复合物（图1B）。
4. 我们还包括了23个先前被报道影响相对多腺苷酸化位点使用的基因（表S1）。

Para_02

1. 我们在 HEK293FT 细胞中进行了合并的 CRISPR 干扰 (CRISPRi) 筛选，并使用了 Perturb-seq 实验工作流程（STAR 方法）来同时捕获每个细胞接收到的指导物的身份以及 3' scRNA-seq 读数（图 1A 和 S1）。
2. 我们的主要分析集中在了深度剖析的 HEK293FT 数据集上（每次扰动中位数 2,168 个细胞），但也在 K562 细胞中重复了实验（每次扰动中位数 960 个细胞）。
3. 在每种细胞系的两个生物学重复实验（独立的病毒转导）中，我们总共获得了 140,415 个单细胞（表 S2），在这些细胞中我们成功分配了一个 sgRNA。
4. 我们应用了我们之前开发的计算管道 Mixscape，48 来解决在合并的单细胞 CRISPR 筛选中已经描述过的混淆变量源的问题。33,34,48
5. 对于 42 个调节因子中的 6 个，Mixscape 将所有细胞分类为"未扰动"，这表明尽管目标基因发生了敲低（KD），但对转录组的影响很小（图 S2A）。
6. 对于其余的 36 个基因，Mixscape 将 69% 的细胞分类为扰动。

[图片上传失败...(image-2b4cd2-1727280015646)]

• 图S1。使用3 'scRNA-seq分析多腺苷酸化位点使用情况，与图1相关
• （A）直方图显示CPA-Perturb-seq数据集中每个基因鉴定的多腺苷酸化位点数量（与polyA_DB v3.2相交后）。
• （B）端粒和内含子中检测到的多腺苷酸化位点相对于典型切割动机AATAAA和ATTAAA的位置偏好。
• （C）用于在单核苷酸分辨率下鉴定切割位点的基因特异性3 'cDNA末端（3 'RACE）放大的示意图。
• （D）从3 'RACE和Illumina测序获得的含有多腺苷酸化序列（切割位点半径的证据）的读取频率（y轴）的视觉化。蓝色条表示预测的转录本，红色箭头指示从CPA-Perturb-seq推断的多腺苷酸化位点。
• （E和F）CPA-Perturb-seq来自JKAMP（左）和CIAPIN1（右）位点的读取覆盖图。每个轨迹代表一个伪批量平均细胞，按目标基因扰动分组。JKAMP在NUDT21扰动下表现出近端位点使用的强烈增加，而CIAPIN1在NUDT21扰动下不表现出多腺苷酸化位点的变化。
• （G）将CPA-Perturb-seq中的NT与NUDT21敲低进行比较时，近端位点使用的变化（x轴）与7个基因的3 'RACE（y轴）。（H）在使用不同技术在HEK293FT细胞中基线分析多腺苷酸化位点使用情况时，显示重复之间可重复性的散点图。每个点代表重复1（x轴）与重复2（y轴）中的百分比等位基因使用。点按核密度估计（截断于2.5）着色，并计算重复之间的皮尔逊相关系数。（I）皮尔逊相关系数矩阵显示在分析多腺苷酸化位点使用情况时，三种不同技术之间在整个转录组范围内的定量一致性。

[图片上传失败...(image-30352a-1727280015646)]

• 补充图 S2. CPA-Perturb-seq 数据集中的扰动响应，与图 1 相关
• （A）Perturb-seq 数据中观察到的目标基因的对数倍变化（左）与该扰动的差异表达基因数量（右）。在未检测到或仅检测到少数差异表达（DE）基因的情况下，Mixscape 将所有细胞归类为"未扰动"。这可能发生在目标扰动不成功的情况下（例如，CLP1），也可能发生在目标基因被下调但未观察到全局效应的情况下（例如，PABC4）。
• （B）CBX3 位点的读段覆盖图。每条轨迹表示细胞的伪批量平均值，按其扰动分组，并按其单独的 sgRNA 分开。
• （C）聚腺苷酸位点/调节因子计数矩阵的层次聚类。每列表示在减去 NT 对照细胞伪批量平均值后，受到相同目标基因扰动的单细胞伪批量平均值。
• （D）CBX3 位点单个细胞内远端聚腺苷酸位点使用比例的分布，按每个细胞内 CBX3 的总读段数划分，包括 NT 细胞（顶部）和 NUDT21 扰动的细胞（底部）。该位点覆盖度低的细胞（1-5 reads，左）在使用近端和远端位点之间显示出主要的双峰分布，但具有足够测序深度的细胞（中，右）表明即使是个别细胞，聚腺苷酸位点的使用也是异质性的。
• （E）3 个 NT 细胞（灰色）和 3 个 NUDT21 扰动细胞（绿色）在 CBX3 位点的读段覆盖情况，说明了单个细胞内聚腺苷酸位点使用的异质性。

Para_02

1. 我们利用单细胞RNA测序数据来量化每个细胞的基因表达和转录组范围内的聚腺苷酸位点使用情况（方法详见STAR Methods）。
2. 我们首先使用polyApipe来识别一组可能的聚腺苷酸位点，然后量化它们在单个细胞中的使用情况，从而为下游分析生成一个聚腺苷酸位点/细胞计数矩阵。
3. 我们的分析仅限于距离polyA_DB v3.2中确定的聚腺苷酸位点50个核苷酸以内的聚腺苷酸位点，polyA_DB v3.2是一个从多个人类细胞系生成的聚腺苷酸位点数据库。
4. 我们还只包括位于基因的内含子或最后一个外显子中的位点（方法详见STAR Methods）。

Para_03

1. 我们共鉴定出35,882个多腺苷酸化位点，分布在12,617个检测到的基因中。
2. 我们发现8,558个基因在我们的数据集中表现出使用两个或更多的多腺苷酸化位点（其中5,661个基因表现出使用三个或更多）（图S1A）。
3. 大多数多腺苷酸化位点包含在剪切位点上游的典型AATAAA/ATTAAA剪切基序，正如预期的那样（图S1B）。
4. 我们使用3′RACE结合Illumina测序在7个位点验证了我们预测的多腺苷酸化位点，发现我们预测的多腺苷酸化位点和3′转录本末端敏感映射高度一致（STAR方法；表S3；图S1C和S1D）。
5. 我们在NUDT21扰动细胞中重复了3′RACE实验，以评估我们检测多腺苷酸化位点使用变化的能力，当比较NUDT21扰动在3′RACE和CPA-Perturb-seq数据之间的效果时，观察到高度一致（R=0.86）（图S1E-S1G）。
6. 我们还发现我们估计的多腺苷酸化位点使用比例与使用批量黄金标准检测PAPERCLIP和A-seq获得的转录组范围定量之间有高度相关性（STAR方法；图S1H和S1I）。

Para_04

1. 我们观察到了不同调控因子在多聚A位点使用上的多样化影响（图1C和S3），并且这些变化在生物学重复实验和多克隆sgRNAs（每个基因3-4个；图S2B）中均可重现。
2. 我们将受干扰细胞的polyA位点/细胞计数矩阵与4336个NT对照一起作为线性判别分析（LDA）、UMAP可视化（图1D）和无监督聚类分析polyA位点矩阵（图S2C）的输入。
3. 这些分析揭示细胞不仅根据扰动聚集，而且还聚集到更广泛的复合物中。
4. 例如，NUDT21和CPSF6（CFIm复合物的两个成员）扰动后的细胞轮廓高度相关，CPSF（CPSF1-4和FIP1L1）、CSTF（CSTF1/3）和多聚酶相关因子（PAF）（PAF1、CTR9、LEO1和CDC73）复合物的成员轮廓也是高度相关的。

[图片上传失败...(image-6ccf16-1727280015645)]

• 图S3：CPA-Perturb-Seq数据集中代表性的读覆盖图，与图1相关（A-G）读覆盖图描绘了代表性位点上替代性polyA位点的差异使用情况。每条轨迹代表了一个根据目标基因扰动分组的细胞的伪批量平均值。所选位点基于图1、2、3、4和5中显示的读覆盖图，并展示了36个调节因子和NT对照细胞的数据。在ATP6V1G1位点（A）上，我们主要观察到调节因子之间的总基因丰度的变化，而不是相对polyA位点使用情况的变化。

Para_04

1. 这些结果表明，我们的数据集可以用来发现特定于复制的共同调控的多聚腺苷酸化位点"模块"，每一个都对功能相关的调控因子的扰动有响应。
2. 然而，我们注意到，多聚腺苷酸化位点/细胞计数矩阵的变化可以反映多聚腺苷酸化位点利用的变化以及基因总体丰度的变化。
3. 例如，当敲低CSTF3时，我们发现了基因近端多聚腺苷酸化位点利用的变化仅与总RNA丰度的变化相对应的情况（ATP6V1G1），仅与由于3' UTR缩短导致的转录本长度变化相对应的情况（HNRNPH3），或者与丰度和相对异构体使用的变化都相关的情况（CDK1）。

[图片上传失败...(image-4958b6-1727280015645)]

• 图2：PolyA残留物在单细胞分辨率下量化替代性多腺苷酸化。
• （A）NT细胞（x轴）和CSTF3扰动细胞（y轴）中6019个近端多腺苷酸化位点的平均使用情况。只考虑至少有两个串联多腺苷酸化位点的基因。不同条件之间的变化可能反映相对多腺苷酸化位点使用情况、总基因表达量或两者的变化。
• （B-D）在（A）中突出显示的三个位点的读覆盖图。阴影标记近端多腺苷酸化位点。
• （E）计算多腺苷酸化残留物的程序示意图（完整描述见STAR方法部分）。
• （F-H）小提琴图描绘了NT和CSTF3扰动细胞的单细胞基因表达水平（左）或近端多腺苷酸化位点的单细胞多腺苷酸化残留物。RNA比较中的不显著（NS）表示绝对log2FC<0.25或使用威尔科克森秩和检验的Bonferroni调整p值>0.05。多腺苷酸化残留物比较中的NS表示差异多腺苷酸化分析中描述的百分比变化<0.05或调整后的p值>0.05（详见STAR方法部分）

Quantifying relative polyadenylation levels at single-cell resolution

在单细胞分辨率下量化相对多腺苷酸化水平

Para_01

1. 为了具体描述扰动驱动的替代性多腺苷酸化效应，我们试图设计一种计算方法来解析这两种效应。
2. 虽然计算基因内每个位点的多腺苷酸酸计数值之比通常用于批量分析中研究替代性多腺苷酸化，但在单细胞RNA测序数据中计算这些比例通常是不可行的或噪声较大，这是由于数据的稀疏性。
3. 我们清楚地观察到了同一细胞内多个多腺苷酸酸位点的使用，但在较低的测序深度下，单细胞多腺苷酸位点使用比例存在噪声（见图S2D和S2E）。
4. 相反，对于每个单细胞中的每个多腺苷酸位点，我们的目标是建模与对照细胞相比的过度使用或使用不足的程度。

Para_02

1. 我们注意到，这个问题在概念上与量化单个细胞中基因表达的变化相似，正如我们和其他人使用广义线性模型所解决的问题一样。
2. 我们扩展了这一框架以模拟替代性多腺苷酸化（图2E）。
3. 我们使用了狄利克特-多项分布来模拟NT细胞中多腺苷酸酸位点使用的背景分布，同时控制基因表达。
4. 与标准多项式相比，狄利克特多项式允许过度分散，这种做法与使用负二项分布来模拟基因丰度时的泊松过度分散类似。
5. 这种过度分散考虑了背景群体中自然生物异质性和"内在"噪声。
6. 正如在sctransform中所做的那样，我们首先为每个多腺苷酸位点单独参数化过度分散估计，但随后对这些估计进行正则化处理，以适应类似位点（STAR方法）。
7. 我们的程序的输出是每个多腺苷酸位点的统计模型，描述了其在4,336个NT对照细胞中的背景使用情况。

Para_03

1. 通过比较每个细胞中每个聚腺苷酸位点观察到的计数与Dirichlet多项式模型预期的值和方差，我们计算了每个聚腺苷酸位点的Pearson残差（聚腺苷酸残差）。这个残差的符号和大小描述了每个细胞相对于每个聚腺苷酸位点的背景分布的相对偏差。正残差反映了相对于背景分布，单个细胞中某个聚腺苷酸位点的使用频率更高，反之亦然。我们使用线性模型检测聚腺苷酸残差的变化，使我们能够识别与扰动相关的聚腺苷酸位点使用变化，而不会受到基因丰度变化的影响（STAR方法）。当应用于我们之前的示例（图2F-2H）时，这种方法成功地识别出了我们在其中观察到转录结构、丰度或两者都发生变化的位点。使用基于LDA的聚腺苷酸残差矩阵可视化（图S4A）确认我们的观察到的共调节模式是由聚腺苷酸变化的协调变化驱动的。

[图片上传失败...(image-408e3d-1727280015645)]

• 图 S4：与图 3 和 4 相关的串联和内含子多腺苷酸化扰动驱动变化。
• （A）通过 CPA-Perturb-seq 分析的 HEK293FT（左）和 K562（右）细胞的 UMAP 可视化。
• 细胞根据目标基因的身份进行着色（每个基因的颜色与图 1C 中使用的颜色相同）。
• 可视化是基于串联和内含子多腺苷酸化位点的多腺苷酸残差线性判别分析（LDA）计算得出的。
• （B）与图 3A 相同，但来自 K562 数据集。
• （C）箱线图表示 NUDT21 扰动后的基因表达观察到的 log2 倍变化（HEK293FT 细胞）。
• 基因根据 NUDT21 扰动后观察到的 3' UTR 变化的程度被划分为十分位数。
• （D）内含子多腺苷酸化位点使用发生变化的基因数量，根据 HEK293FT（左）和 K562（右）细胞内含子位点的使用增加或减少进行分类。
• 调节因子 PAF1、CTR9、CDC73、SCAF8、SCAF4、CPSF3L、RABP2 和 PABPN1 主要表现出内含子多腺苷酸化位点的使用增加。
• （E）展示 HEK293FT 细胞内含子多腺苷酸化调节因子之间关系的皮尔逊相关矩阵。
• 相关性是使用差异多腺苷酸化分析期间学习的线性模型系数计算的（STAR 方法）。
• 基因通过层次聚类进行排序。
• （F）与 NT 对照细胞相比，PAF1 扰动（左）和 CPSF3L 扰动（右）细胞的多腺苷酸化位点使用百分比变化（y 轴）。
• 这作为与下一个多腺苷酸位点距离的函数（在 x 轴十分位上分箱）。
• 与下一个多腺苷酸位点的距离与内含子位点对 PAF1 扰动的响应相关，但与 CPSF3L 扰动无关。
• 皮尔逊相关性是在距离（对数刻度）和百分比变化之间计算的。
• （G）与（F）相同，但对于 PAF1 扰动，分别显示位于第一个内含子（左）和不位于第一个内含子（右）的位点。
• （H）与（F）相同，但将内含子多腺苷酸化位点按 GC 含量十分位进行分箱。
• （I）与（F）相同，但将内含子多腺苷酸化位点按内含子宽度十分位进行分箱。
• （J）热图显示 K562 的模块 A 基因和模块 B 基因的远端位点的多腺苷酸残差。
• 热图中显示的多腺苷酸位点身份和顺序与图 4C 相同。
• （K）与图 4F 相同，显示 K562 细胞中模块 A 和 B 远端位点的使用。

Characterizing perturbation-dependent changes in polyadenylation

表征依赖于扰动变化的聚腺苷酸化

Para_01

1. 我们鉴定出7,402个基因在至少36个基因扰动中的至少一个中表现出差异性的选择性多腺苷酸化（至少有一个多A位点的使用差异）。
2. 但是，我们在调控因子之间观察到了实质性的差异。
3. 例如，CFIm复合体的成员NUDT21表现出最强的扰动响应（图3A），包括多A位点使用和总转录本丰度的广泛变化。
4. 相比之下，PABPC1的扰动（该蛋白在核输出后与多A尾结合）主要影响转录本水平的变化，而不是结构（图3A和S4B）。

[图片上传失败...(image-1dbf4a-1727280015645)]

• 图3. 在CPA-Perturb-seq中定性串联和内含子替代多聚腺苷酸化(A)(左) HEK293FT数据集中每个调节因子扰动后多聚腺苷酸化位点使用、基因表达或两者都有变化的基因数量。(中) 扰动驱动内含子或串联多聚腺苷酸化位点使用发生变化的基因数量。(右) 串联多聚腺苷酸化位点使用发生变化的基因数量，按3' UTR缩短或3' UTR延长分类。
• （B和C）读覆盖图显示ZSCAN9（B）和EXOSC4（C）位点内含子位点（方框内）的差异使用。
• （D）热图显示在PAF复合物成员（PAF1/CTR9/CDC73）、反终止子（SCAF4/SCAF8）、PABPN1和CPSF3L/RPAP2扰动后唯一差异使用的内含子位点的多聚腺苷酸残差。每个热图单元格显示按sgRNA身份分组的细胞后的伪批量平均值。
• （E）热图显示预测每种扰动内含子多聚腺苷酸位点使用时不同特征的重要性。颜色代表从预测线性模型（STAR方法）获得的每个协变量的t统计量，也用于层次聚类。
• （F）元基因图显示每个调节因子内含子多聚腺苷酸位点显著变化的规范化位置。
• （G）对于每个调节因子，使用有显著变化的内含子多聚腺苷酸位点的内含子的GC含量。
• （H）对于每个调节因子，包含有多聚腺苷酸位点使用显著变化的内含子的宽度。
• （I）对于每个调节因子，按Mukherjee等人的分类，按内含子剪接速度分组，包含多聚腺苷酸位点使用显著变化的内含子的比例。参见图S4。

Para_01

1. 我们接下来根据聚腺苷酸化位点使用情况的变化将其分类为内含子聚腺苷酸化或串联聚腺苷酸化，对于串联位点，我们确定了它们是代表使用增加还是近端（缩短）或远端（延长）位点（STAR方法；图3A）。大多数调控因子影响串联聚腺苷酸化，表现出超过70%的偏向于缩短或延长，这也在我们K562数据集中得到了验证（图S4B；表S4）。3' UTR缩短通常与总基因丰度的增加有关（图S4C），这与3' UTR长度与可能影响RNA稳定性的调控元件存在的关系相一致。
2. 一般情况下，待整理文本中有多少句话，输出的json就应该有多少个item。

Para_02

1. 将我们的结果与之前利用批量3'端测序技术的研究进行比较，突显了Perturb-seq技术在定义扰动特征方面的优势。之前的研究一致发现，NUDT21扰动影响了一部分基因（范围在375到1600个）的多聚腺苷酸位点使用，并导致串联UTR的3' UTR缩短。
2. 在我们的HEK293FT和K562数据集中（分别见图3A和S4B），我们观察到高灵敏度（超过5400个基因在扰动后多聚腺苷酸位点使用发生了显著变化）以及高特异性（超过91%的串联UTR变化导致缩短）。
3. 同样，RBBP6扰动与3' UTR延长有关，但我们的数据集以高特异性识别出更多基因（超过94%的3' UTR缩短，见图3A和S4B）。

Distinct regulatory mechanisms drive intronic polyA site usage

不同的调控机制驱动内含子多A位点的使用

Para_01

1. 我们确定了八个调节因子，其中干扰反应主要与内含子多腺苷酸化增加相关（图3A和S4D）。
2. 这些包括反终止蛋白SCAF8和SCAF4，22 PAF复合物成员66（PAF1、CDC73和CTR9）、PABPN1、RPAP2和CPSF3L，这是整合器复合物的一个成员67。
3. 这些调节因子每一个都直接与RNA聚合酶II或新合成的转录本相互作用，68,69,70,71,72但是这些调节因子的干扰影响了不同的内含子位点。
4. 这表明内含子多腺苷酸化有多种调节模式（图3B-3D）。

Para_02

1. 我们对具有相关扰动响应的调控因子进行了层次聚类分析（图 S4E）。
2. 我们发现 PAF 复合物成员（PAF1、CDC73 和 CTR9）调控类似的内含子多 A 位点（图 3D），并且整合素复合物成员 CPSF3L 和相关蛋白 RPAP2 之间存在类似的扰动响应。
3. 我们进行了线性建模以识别预测内含子多 A 位点使用变化的特征（图 3E），并发现特征预测强度在不同调控因子间差异很大。
4. PABPN1 的扰动与内含子首个多 A 位点的使用增加强烈相关（图 3F）。
5. 与这一发现一致的是，PABPN1 作为一个串联的 3' UTR 调控因子发挥作用（图 3A），但它也通过与识别 5' 帽结构的因子结合，参与核监控和降解短的聚腺苷酸化转录本（参考文献 68）。

Para_03

1. 我们发现，内含子GC含量以及转录本中到下一个切割位点的距离（衡量切割事件间延伸时间的一个指标）都能预测对PAF1扰动的反应性。这种关系在位于第一个内含子的polyA位点最强，但也适用于下游位点（图S4F和S4G）。在我们的数据集中，只有21%的受调控的内含子polyA位点对PAF1扰动有反应，这突显了这一特定位点子集对延伸动力学的广泛变化有反应。
2. 这种关系在位于第一个内含子的polyA位点最强，但也适用于下游位点（图S4F和S4G）。
3. 在我们的数据集中，只有21%的受调控的内含子polyA位点对PAF1扰动有反应，这突显了这一特定位点子集对延伸动力学的广泛变化有反应。

Para_04

1. 整合子亚单位CPSF3L及其相关因子RPAP2在整个基因体上调节内含子多聚A位点（图3F）。这些响应位点主要位于短内含子（中位长度：2,660 bp）和GC含量较高的内含子中（图3G、3H、S4H和S4I）。整合子复合物与多种不同的功能有关，包括参与小核RNA（snRNA）生物发生67以及在识别暂停的启动子近端RNAPII时驱动提前终止72。由于我们没有观察到内含子富集在转录起始位点附近（图3F），因此我们认为可能是snRNA处理缺陷导致了更广泛的剪接异常，影响了内含子多聚腺苷酸化。我们利用之前发表的RNA代谢数据集，计算了我们数据集中2,212个基因的内含子剪接速度（图3I）62。我们发现CPSF3L/RPAP2特征富集于具有"慢"剪接动态的内含子，这与其升高的GC含量77一起表明它们剪接效率低下。

Para_05

1. 最后，我们发现通过干扰两种抗终止蛋白，响应内含子多A位点的身份出现了明确的分支。SCAF4主要调节短内含子（中位长度：7,529 bp）中多A位点的使用，这些内含子具有高GC含量，而SCAF8则调节长内含子（中位长度：33,798 bp）中位点的使用，这些内含子具有低GC含量（图3G和3H）。我们的发现扩展了之前的开创性工作，即SCAF4和SCAF8通过冗余机制阻止内含子多A位点的使用。通过我们检测方法的敏感性，我们能够检测到两种单独干扰后的表型，证明了这些蛋白质在数百个位点上的非冗余作用，并识别了指导这种选择性的决定因素。我们的分析表明，内含子多聚腺苷酸化不是一个全局性调控的现象，不同的位点集合对调控不同RNA核生命周期组分的因素的干扰具有独特的敏感性。

Modules of co-regulated tandem polyA sites exhibit distinct functional properties

共调节串联多A位点的模块表现出不同的功能特性

Para_01

1. 我们对串联多腺苷酸化调节因子的层次聚类分析进行了重复，以确定相关扰动的响应(STAR方法；图4A)。扰动簇反映了核心CPA复合物的成员结构以及共调节的额外证据。
2. 例如，RBBP6、FIP1L1和PCF11不是同一复合物的成员，但它们的扰动导致在重叠位点3' UTR延长。
3. 我们在K562数据集中观察到了高度一致的相关性(图4B)。

[图片上传失败...(image-cf1583-1727280015644)]

• 图4显示了共调节的多腺苷酸化位点模块在功能上的差异。（A）展示在HEK293FT细胞中串联干扰之间的关系的皮尔逊相关系数矩阵。相关性是通过使用差异多腺苷酸化分析期间学习的线性模型系数来计算的（STAR方法）。基因通过层次聚类进行排序。（B）与（A）相同，但相关系数矩阵是由K562多腺苷酸化残差生成的。（C）热图显示模块A基因（CSTF和CPSF与CPSF6/NUDT21的方向相反）和模块B基因（CSTF和CPSF与CPSF6/NUDT21的方向相同）的远端峰位点的多腺苷酸化残差。每个模块显示了排名前100的多腺苷酸化位点（按CSTF干扰排序）。（D）模块A和B所属基因的示意图。（E）读取覆盖率图显示了属于模块A（左，CCT6A）和模块B（右，TMEM106C）的代表性基因的多腺苷酸化位点使用情况。（F）密度图显示了NT对照细胞中属于模块A（左）与模块B（右）的基因的远端位点使用情况。另见图S4。

Para_01

1. 我们发现，相关性结构并非仅仅由全球性偏好于缩短和延长驱动，还受到位点特定差异在扰动响应中的影响。
2. RBBP6的扰动（偏好于3' UTR延长）以及CFIm复合体成员CPSF6和NUDT21的扰动（偏好于3' UTR缩短）表现出强烈的反相关响应，反映出全局性的相反调控。
3. 相比之下，CSTF和CPSF复合体成员（偏好于3' UTR延长）与CFIm成员仅表现出较弱的反相关性，反映出更为复杂的共调控模式。

Para_02

1. 为了进一步探索这个问题，我们考虑了一组在CFIm扰动后转录缩短的基因（STAR方法）。
2. 我们根据CSTF扰动反应将这些基因进行了分类，并观察到了一个预期的模块（图4C-4E，模块A），该模块包含323个多聚A位点（20%），在这些位点上，CSTF扰动导致了相反的延长反应。
3. 然而，我们还确定了一个包含149个基因的模块（模块B）（9%），其中CSTF扰动导致了缩短，这种现象模仿了CFIm扰动。
4. 剩余的71%的位点在CSTF扰动后未表现出利用上的变化。
5. 我们在K562细胞中的相同位点观察到了可重复的模式（图S4J）。

Para_03

1. 此外，我们发现，在CSTF扰动时3' UTR延长（模块A）的基因在NT细胞中强烈倾向于近端位点的使用，而具有相反CSTF扰动表型的基因（模块B）则表现出远端位点偏差（图4F和S4K）。这些基因远端偏差程度的变异很可能是由它们近端和远端位点CSTF活性的相对强度差异驱动的。

APARENT-Perturb reveals an interactive cis-regulatory code

APARENT-Perturb 揭示了一个交互性的顺式调控代码

Para_01

1. 我们发现的差异多聚腺苷酸化可重复模式强调了局部序列在决定多聚A位点对干扰反应中的作用。
2. 我们试图扩展准确预测基线条件下全基因组替代多聚腺苷酸化模式的深度学习模型，以预测干扰反应。
3. 这些模型成功捕捉非线性相互作用的能力，包括基序之间的位置和组合依赖性，突显了它们学习复杂顺式调控决定因素的能力。

Para_02

1. 为了预测未受干扰细胞中的基线聚腺苷酸化位点使用情况，我们使用了APARENT2模型，这是一个最初在HEK293FT细胞中测量的MPRA数据集上训练的残差神经网络。
2. 受到MTSplice模型的启发，我们随后训练了一种新的基于集成学习的多任务干扰网络（APARENT-Perturb），它能预测我们在10个最强的干扰中聚腺苷酸化位点的使用情况，输入的是200个核苷酸序列与核心六聚体的对齐以及APARENT2基线得分（见图5A）。
3. APARENT-Perturb准确预测了在NT条件下保留基因的聚腺苷酸化位点的同种形式比例（RS=0.70），以及在干扰中的比例（0.65≤RS≤0.73），这是通过10倍交叉验证测量的（见图5B和S5A）。

[图片上传失败...(image-ae3e99-1727280015644)]

• 图 5. 一个多任务神经网络预测 RNA 序列的扰动响应
• （A）APARENT-Perturb 的示意图，这是一个基于集成的神经网络架构，用于预测扰动响应。绿色/蓝色/红色输出头对应于模型对 K 个扰动条件的预测。
• （B）在预测远端异构体比例（顶部行）或相对于 NT 条件的远端异构体比例差异（底部行）时的 10 折交叉验证性能。
• （C）KMT5A 基因中 2 个示例扰动的特定序列归因分数。归因分数是在与 NT 细胞计算残差后显示的。
• （D）10 个扰动位置的平均归因分数。对于每个扰动，展示了 3 个顶级 MoDISco 基序（STAR 方法）。
• （E）热图显示每个基因最远端位点每个扰动位置的平均归因分数。
• （F）CSTF1 扰动在模块 A 与模块 B 中对近端（左）和远端（右）位点的平均归因分数。核心六聚体和下游序列元素（DSEs）的位置分别用实线和虚线垂直线标记。图表显示了单碱基对分辨率的平均归因分数（点）以及 loess 平滑趋势（线）。
• （G）针对 GC 丰富（红色）或 AT 丰富（蓝色）背景中的双重 TGTA 基序的表型分析。y 轴反映了在预测 NUDT21 扰动时，同时插入两个基序与一次插入一个基序的效果比较（STAR 方法）。
• （H）基于 RBBP6 扰动的典型六聚体和 GT 丰富基序的表型分析。另见图 S5。

[图片上传失败...(image-d42dc3-1727280015644)]

• 图 S5. APARENT-Perturb 的评估与解释，与图 5 (A) 相关。（左）在三种扰动条件（NUDT21、CSTF3 和 RBBP6）下的预测与实测远端聚腺苷酸位点使用情况。（右）相对于非目标（NT）条件的预测与实测近端或远端聚腺苷酸位点使用差异。仅包含交叉验证过程中保留集的预测。
• （B）与图 5C 相同，但适用于其他位点。
• （C）所有 10 个学习扰动的平均残差归因得分，按近端和远端聚腺苷酸位点分开。每个图下方显示了与近端和远端位点对应的排名前 5 的 MoDISco 基序。这些基序都与正贡献得分相关（即增加扰动幅度），除了 RBBP6 的远端基序和 THOC5 的近端基序；这些基序与负贡献得分相关（减少扰动幅度）。
• （D）按扰动分组的特定类别 MoDISco 基序命中的平均归因得分。红色表示基序类别与扰动幅度的增加相关。蓝色表示基序类别与扰动幅度的减少相关。
• （E）（左）NUDT21 和 CSTF1/3 扰动中具有高归因得分的序列特征共存情况。NUDT21 扰动中远端 pA 信号的 平均归因得分（红色）以及 NUDT21 扰动与 CSTF1（橙色）或 CSTF3（紫色）之间归因得分的平均差异。（右）NUDT21 和 CSTF3 扰动中符合 U/G 富含基序的 MoDISco 基序命中共存分析（用费舍尔精确检验计算的优势比和 p 值）。
• （F）与图 5D 相同，但显示了响应 NUDT21 扰动的所有远端聚腺苷酸位点。聚腺苷酸位点根据模型分配的重要性分数（DSE；用垂直线表示）分为十分位。对于第 1 十分位的位点，DSE 与 NUDT21 的正归因得分相关，而第 10 十分位的位点具有负贡献得分。
• （G）CPA-Perturb-seq 数据集中调节因子在十分位 1 和十分位 10 的聚腺苷酸位点的平均聚腺苷酸残差（按扰动缩放）。这些十分位是由 NUDT21 扰动模型定义的，但也对不同 CSTF 扰动有不同的响应。∗表示 t 检验中十分位 1 和十分位 10 的残差 p 值 <0.05；∗∗表示 p 值 <0.001。
• （H）在远端聚腺苷酸信号中精确两个野生型 UGUA 基序的组合消融，这表明了整体亚加性交互。用指数函数估计的拟表型优势比（y 轴），即在 NUDT21 条件下预测的扰动对数优势比，当同时替换两个基序为随机序列时与一次替换一个基序时对数优势比之和的差异。
• （I）在变化的侧翼核苷酸环境下对双重 UGUA 基序的插入模拟。y 轴表示在 NUDT21 条件下插入两个基序和一次插入一个基序时预测的扰动对数优势比的指数函数差异。
• （J）典型核心六聚体基序和下游 U/G 富含基序的组合消融。将 U/G 富含基序替换为随机序列，而核心六聚体替换为随机选择的一位点突变的六聚体。
• （K）当在 NUDT21 条件下相对于 NT 对远端扰动对数优势比拟合回归模型时，计数特征和组合指标变量的线性回归系数。系数分布是从 1000 次引导重采样生成的。
• （L）当对 RBBP6 扰动对数优势比拟合时的回归系数分析。

Para_02

1. 接下来，我们进行了计算机模拟突变（ISM），得到了一系列核苷酸水平的"归因得分"，这些得分反映了每个碱基对模型预测的贡献。
2. 重要的是，通过减去NT（基线）输出的得分，我们可以孤立每个序列在预测扰动响应中的重要性。
3. 例如，KMT5A基因中远端多A位点的归因得分突出了预测能够驱动对NUDT21扰动响应的上游TGTA基序和驱动对CSTF3扰动响应的独特的下游GT丰富区基序（图5C和S5B）。
4. 对于每种扰动，我们平均了跨位点的ISM得分以识别含有重要序列元件的区域（图5D和5E）。
5. 然后，我们使用了一个基序发现工具TF-MoDISco，将每种扰动的归因得分聚类成一组显著的基序（图5D，S5C和S5D）。
6. 这些结果概括并扩展了之前建立的结合基序和位置。
7. 例如，CSTF1和CSTF3在下游区域展示了T-或GT丰富序列的重要性峰。
8. NUDT21和CPSF6在多A位点的上游区域展示了高平均重要性，但也用A-和T-丰富的翼端扩展了规范的TGTA基序。
9. 有趣的是，NUDT21和CPSF6的扰动响应还受到大约30-50 bp下游的序列元件的驱动（下游元件[DSE]）。
10. 这个DSE与预测对CSTF扰动重要性的区域重叠。
11. 这反映了两者在相同位点上的功能序列共同富集，表明这些因子之间存在调控相互作用。

Para_03

1. 我们之前观察到CSTF和CFIm复合物成员可以共同调节polyA位点，无论是在同一方向还是相反方向（图4C-4F）。
2. 我们发现，在CFIm扰动导致转录缩短而CSTF扰动导致延长的基因（模块A）中，近端polyA位点的DSE由具有高CSTF归属分数的序列元素组成。
3. 然而，在两个复合物扰动都导致转录缩短的基因（模块B）中，近端位点的序列元素显著较弱（图5F，左，p < 2.0 × 10−5，威尔科克森双尾秩和检验）。
4. 相比之下，我们观察到模块B基因在远端位点的归属分数增加（图5F，右，p < 1.6 × 10−4）。
5. 综合这些发现，我们提出了一个模型，在该模型中，近端polyA位点的序列内容对于建立近端偏向和扰动响应尤其重要。

Para_04

1. 最后，我们模拟了单个和成对的主题插入，以识别CPA调节因子之间的上位性相互作用，正如转录因子所做的那样。
2. 例如，CFIm复合物包括一个NUDT21同源二聚体，但目前尚不清楚多个TGTA基序是否以及如何影响结合。
3. APARENT-Perturb给具有A-和T-丰富侧翼的NUDT21基序赋予了更高的重要性得分，因此我们在进行插入时测试了多种可能的侧翼序列。

Para_05

1. 当在短距离内插入相邻的TGTA基序时，我们观察到当两个基序都被GC丰富的序列包围时，NUDT21的扰动有协同效应，而AT丰富的背景与次加性相互作用相关（图5G、S5H和S5I）。
2. 插入典型的核心六聚体和GT丰富的DSE元素，这两种元素都与RBBP6和CSTF的调控有关，展现出一种协同的表型关系，这种关系在20-bp插入距离时达到最大（图5H和S5J）。
3. 我们使用多项式特征回归验证了这些结果（图S5K-S5L）。
4. 我们得出结论，将深度学习模型应用于Perturb-seq数据集可以揭示顺式调控景观，该景观编码了在多个复合物中复杂的共调控模式。

Validating sequence predictions by massively parallel screening with perturbations

通过大量并行筛选扰动来验证序列预测的有效性

Para_01

1. 我们的结果表明APARENT-Perturb能够识别驱动聚腺苷酸化位点选择的序列，将这些序列分配给特定的调节因子，并识别它们之间的相互作用。
2. 在任何基因座验证这些预测都需要两个组成部分：证明突变预测重要性高的序列元素会改变聚腺苷酸化位点的使用，并且这些改变依赖于所分配的调节因子的存在。
3. 我们在373个基因座设计了修改后的序列（总共产生3,802个野生型[WT]和突变序列），并利用之前描述的报道基因构建体，进行了一些小的修改，以可扩展地验证APARENT-Perturb的预测（见STAR方法；图6A和S6A）。
4. 我们在有和没有CRISPRi扰动CSTF3的样本中进行了MPRA，使我们能够探索每个序列突变与基因扰动的影响。

[图片上传失败...(image-c301ca-1727280015644)]

• 图6.通过在多种遗传背景下进行大规模平行报告分析(MPRA)验证APARENT-Perturb。
• （A）用于验证APARENT-Perturb的大规模平行报告分析（MPRA）示意图。
• （B）根据APARENT-Perturb预测的373个WT位点的近端聚腺苷酸位点使用情况（对数几率）（x轴）和通过MPRA测量的结果（y轴）。预测在NT（顶部）和CSTF3扰动（底部）样本中均准确，对于近端和远端位点均准确。
• （C）比较CSTF3和NT样本的近端聚腺苷酸位点使用情况（对数几率比）的变化，针对APARENT-Perturb预测的响应CSTF3扰动的序列（n = 107，右侧）与不响应的序列（n = 109，左侧）。**表示p值<0.0001，威尔科克森检验比较响应与中性序列的对数几率比。
• （D）基于CPA-Perturb-seq推断的切割位点构建的GYG2 3' UTR的基因模型。红色突出显示的区域插入到MPRA构建体中。
• （E）（顶部）APARENT-Perturb CSTF模型在GYG2基因（chrX:2882818）远端位点的归因分数，对于野生型（WT）序列（顶部）以及在进行了序列改变后的结果。这些改变包括重新排列预测的CSTF响应元素以最小化扰动效应（中间）以及设计合成突变以最大化CSTF3响应。彩色核苷酸表示改变的核苷酸。
• （底部）GYG2位点的MPRA数据的可视化。线条表示包含聚腺苷酸序列（聚腺苷酸读数）的读数比例，这表示近端切割，用于NT和CST3F扰动细胞（STAR方法）。
• （F）在107个CSTF3响应位点上执行（E）中描述的序列改变的效果，包括（左）和不包括（右）典型的CSTF结合基序。CSTF3与NT的对数几率比（y轴）对于改善和重新排列的序列改变（x轴）进行了展示。*表示p值<0.0001；表示p值=0.0001-0.05；NS，无统计学意义。
• （G）比较NUDT21与NT样本的近端聚腺苷酸位点使用情况的对数几率比（y轴），对于WT序列以及在（F）中重新排列CSTF序列元素后的结果。**表示p值<0.0001。
• （H）FAM13C位点上聚腺苷酸读数的比例，表示近端切割，在NUDT21扰动和NT细胞中，对于WT序列以及在重新排列CSTF序列元素后的结果。
• （I）WT模块A序列（n = 47，绿色）和将CSTF响应序列元素插入近端聚腺苷酸位点后（右侧，橙色）远端聚腺苷酸位点使用的密度。
• （J）基于CPA-Perturb-seq推断的切割位点构建的PGRM1（模块B基因）3' UTR的基因模型。红色突出显示的区域插入到MPRA构建体中。
• （K）将CSTF响应序列元素插入（J）中描绘的构建体中对聚腺苷酸读数比例的影响（y轴）。
• （L）将NUDT21基序（TGTA）插入中性序列时对聚腺苷酸位点使用的对数几率比（y轴，在NT细胞中插入的对数几率比）的影响，当基序被AT富集（左）和GC富集（右）侧翼包围时，按远端位点。
• （M）在PIP5K1C位点上单独或双重插入TTTGTAAT基序的序列的APARENT-Perturb ISM分数。
• （N）在多个距离（x轴）处插入具有AT富集侧翼的TGTA基序的表型几率比（y轴）对于NT（灰色）和NUDT21（绿色）样本。**表示单侧t检验的p值<0.0001，表型几率比为1。另见图S6。

[图片上传失败...(image-d7e25b-1727280015644)]

• 图S6展示了在多种遗传背景下对数千个APARENT-Perturb预测的分析，与图6相关（A）示意图展示了野生型（WT）和打乱序列的示例构建体，它们被插入到MPRA构建体中，以及来自每个序列的代表性mRNA分子的读取。构建体的前20 bp由一个独特的条形码（红色）组成，确保每个近端位点构建体的可识别性。每个构建体还以一个核心六聚体（绿色）为锚。在所示示例中，我们确定了一个9-bp的下游元件，预测CSTF3扰动响应（橙色），并用随机序列（紫色）替换了该元件。通过直接读取近端位点的使用情况（当我们测序多A尾时，以蓝色突出显示），我们测量了序列修饰如何改变多A位点的使用。
• （B）密度图显示了所有3802个测试序列的远端位点使用率（远端位点使用情况的分数），分为5个远端位点背景。
• （C）散点图表示了在不同生物重复实验中多A位点使用测量（近端位点使用情况的分数）的可重复性，对于NT（左）和CSTF3（中），颜色代表核密度估计。对于独立针对CSTF3的引导序列，我们也观察到了高可重复性（右）。
• （D）与图6F相同，但按远端位点背景分开。
• （E）与图6F相同，但对于APARENT-Perturb预测的CSTF3非响应序列。
• （F）APARENT-Perturb NUDT21模型在FAM13C基因（chr10:59246514）的远端位点的归因得分，对于野生型（WT）序列（顶部）和在预测的CSTF3序列元件（绿色突出显示）打乱后的序列（底部）。上游NUDT21结合位点（TGTA）以橙色显示。
• （G）与图6G相同，但按远端位点分开。
• （H）直方图显示了APARENT-Perturb预测的9-bp区域的最大归因得分距离（核心六聚体的下游），针对CSTF3模型（左）。这些区域内的WT序列中最常见的4-mer包括T-和GT-丰富元件，但不包括NUDT21结合位点TGTA（右）。
• （I）在模块A基因中打乱预测的CSTF响应元件（n = 49，左）或在模块B基因中插入GT-丰富区域（n = 47，右）对近端位点使用产生相反影响（通过比较修改后的序列与WT的log-odds比，y轴），在5个远端位点中一致。
• （J）（左）CSTF3扰动响应（log-odds比，y轴）对于WT模块A基因的近端位点使用和进行序列打乱后的情况。（右）WT模块B位点的扰动响应和进行序列插入后的情况。∗∗表示p值<0.0001。
• （K）插入TGTA基序，带有AT-丰富（左）或GC-丰富（右）侧翼对近端位点使用的影响（插入与WT序列比较的log-odds比，y轴）。∗∗表示p值<0.0001；∗表示p值=0.0001–0.05。

Para_01

1. 我们将3,802个测试位点逐一插入报告基因的近端位点，并使用了带有五个不同强度的远端位点的构建体（方法详见STAR Methods；图S6B；表S5）。
2. 对于每种遗传条件，我们用我们的MPRA库转染细胞，进行了两次生物学重复，并观察到了可重复的多聚A位点定量结果（相关系数R为0.96–0.98；图S7C）。

[图片上传失败...(image-3b6ec3-1727280015644)]

• 图 S7. 从全基因组 Perturb-seq 中鉴定替代性多腺苷酸化调控因子，与图 7 相关（A）与图 7C 相同。坐标轴显示了相关性矩阵中所有基因名称。（B）GWPS 数据集中具有串联多 A 位点使用显著变化的基因数量，根据 3' UTR 缩短或 3' UTR 延长进行分类。（C）GWPS 数据集中具有内含子多 A 位点使用显著变化的基因数量，根据内含子多 A 位点使用增加或减少进行分类。（D–I）覆盖图描绘了每个模块中代表性基因的替代性多 A 位点的差异测序读段覆盖度和使用情况。

Para_01

1. 首先，我们将APARENT-Perturb对CSTF3扰动反应的预测与我们的MPRA进行了比较。
2. 在NT对照细胞（平均R = 0.88）和CSTF3KD细胞（平均R = 0.89）中，我们观察到每个五个远端位点都与预测值明确一致（图6B）。
3. 预测的响应位点在CSTF3扰动后表现出位移（p < 3.4 × 10−37），而预测的非响应位点表现出最小的变化（图6C）。

Para_02

1. 我们接下来测试了APARENT-Perturb是否能够识别特定的序列元素驱动这种关系。
2. 在每个位点上，我们识别并打乱了具有最大ISM分数的9-bp序列元素，用于CSTF3扰动（STAR方法；图6D和6E中的示例）。
3. 这些突变正确地消除了NT和CSTF3扰动细胞之间的差异（图6F和S6D）。
4. 此外，我们还设计了"超响应"位点，类似于基于模型的合成增强子序列对TF结合的设计。86,87
5. 我们对原始位点进行了10个或20个单核苷酸突变，APARENT-Perturb预测这将增加位点的响应性（STAR方法）。
6. 这些序列在CSTF3扰动后的MPRA中显示出对聚腺苷酸位点使用的改变增加（图6F）。

Para_03

1. 我们还发现，APARENT-Perturb 成功地识别出缺乏与 CSTF 结合有关的 GT-或 T-富集基序的 CSTF3 响应元件（参考文献 88,89）（占响应序列的 10%，图 6F）。相反，我们预测并验证为对 CSTF3 扰动无反应的一组序列中包含一个典型的 DSE（图 S6E）。这突显出 APARENT-Perturb 在预测调节子活性方面优于基于通用基序的方法。

Para_04

1. 我们之前的分析（图5D、5E和S5D-S5F）表明，驱动CSTF调控的序列特征对NUDT21的反应性也很重要。
2. 在NUDT21 KD细胞中测试相同的构建体时发现，打乱CSTF基序也会影响对NUDT21干扰的反应性（p < 3 × 10^-11，图6G、6H和S6G）。
3. 这些调控元件位于规范六聚体的下游，并且不包含下游TGTA基序（图S6H）。
4. 与CSTF3干扰相比，这种效应有所减弱，这可能是由剩余的NUDT21调控基序引起的（图S6F）。
5. 这很可能反映了NUDT21和CSTF3之间的由序列驱动的间接调控相互作用。

Para_05

1. 我们还测试了我们的假设，即差异性的CSTF结合决定了我们在模块A与模块B基因中观察到的近端与远端偏差。
2. 为了驱动CSTF3反应，我们在模块B基因的近端位点插入了一个由APARENT-Perturb预测的9个碱基序列。
3. 令人惊讶的是，我们发现这个单个序列的修改逆转了模块之间的差异，将位点从远端偏差（WT序列）转变为近端偏差（突变序列）（图6I-6K和S6I）。
4. 这些更改的效果在CSTF3干扰的细胞中有所减少，突显出CSTF3活性驱动这种行为（图S6J）。
5. 相比之下，打乱模块A近端位点的预测CSTF反应元件导致向远端聚A位点使用的转变（图S6I）。
6. 这证明了APARENT-Perturb除了识别扰动反应外，还成功识别了决定近端与远端偏差的序列元件。

Para_06

1. 最后，我们测试了APARENT-Perturb对序列基序之间相互作用的预测。我们在49个位点上进行了单次和双次的NUDT21 TGTA基序插入，测试了AT富集侧翼和GC富集侧翼。
2. 与APARENT-Perturb的预测一致，与GC富集基序插入相比，AT富集基序的单次插入在聚腺苷酸位点使用上导致了显著更强的变化（p < 3.0 × 10^-54）（图6L）。
3. 这种插入的强度在NUDT21干扰细胞中相比于NT对照组被减弱了（图S6K）。
4. 我们还观察到，当在短距离内进行带有AT富集侧翼的TGTA基序的双插入时，存在次加性效应（图6M和6N），尽管GC侧翼基序的插入显示出最小的影响。

Identification of CPA regulators from genome-wide screening datasets

从全基因组筛选数据集中鉴定CPA调节因子

Para_01

1. 为了了解附加调节因子如何影响聚A位点的使用，我们对最近发表的K562细胞中的全基因组Perturb-seq（GWPS）数据集进行了重新分析。
2. 我们计算了聚A残差，并将其作为差异聚腺苷酸化分析（STAR方法）的输入。
3. 由于GWPS数据集中每个扰动所包含的细胞数量远少于我们的数据集（对于36个重叠的扰动，中位数分别为94个细胞和1142个细胞），因此我们识别出在聚A位点使用方面发生变化的基因数量明显更少（每个重叠扰动中位数为72个基因，而我们的数据中为1389个）。
4. 然而，GWPS数据集仍然能够准确地表征每个调节因子。
5. 例如，我们观察到通过同时扰动串联聚腺苷酸化调节因子，两个数据集都趋向于3' UTR缩短或延长的一致偏差（图7A）。

[图片上传失败...(image-5ec2b1-1727280015643)]

• 图 7. 在基因组规模 Perturb-seq 数据集中相对多聚腺苷酸位点使用异质性的特征化（A）在 CPA-Perturb-seq 数据集中扰动串联调节因子后观察到的 3' UTR 缩短偏好（x 轴）和 GWPS 数据集（y 轴）。
• （B）与图 3B 相同，但对于 GWPS 数据集。
• （C）相关性矩阵描述 GWPS 数据集中的扰动之间的关系，如图 4A。（左侧显示六个相关模块各自的代表基因。所有基因均列于图 S7A。）
• （D）代表多聚腺苷酸位点（n = 100），在扰动外泌体/PAXT 复合物成员后其使用增加。
• （E）CPSF3L 扰动（来自 CPA-Perturb-seq）与 GWPS 数据中的其他整合复合物成员和剪接体因子相关性最高（y 轴）。
• （F）代表读段覆盖图，展示在 CPSF3L、其他整合复合物成员和 SMN 复合物成员扰动后多聚腺苷酸位点使用共享的变化。
• （G）CPA-Perturb-seq 内含子特征在 GWPS 模块响应位点中的富集（y 轴）。点的大小对应于超几何富集检验的 -log10(p 值)，颜色对应于每个特征的 平均多聚腺苷酸残差。
• （H）我们数据集中所有内含子多聚腺苷酸位点的剪接供体位点（5'）MaxEnt 评分与 CPSF3L/RPAP2 扰动中显著增加使用的位点进行比较。∗∗表示 p 值 <0.0001，Wilcoxon 检验。另见图 S7。

Para_01

1. 为了关注直接修饰RNA的调节因子，我们将分析限制在一组1280个RNA结合蛋白上。
2. 我们识别出134个扰动，这些扰动至少在100个基因中引起了多聚A位点使用的改变（图7B；表S4），包括高度相关的扰动（图7C和S7A）。
3. CFIm复合物成员CPSF6和NUDT21以及转录-输出复合物成员THOC3的扰动表现出相关的3' UTR缩短（图S7B和S7D），而由上调移位复合物成员、小核糖体亚基和核糖体成熟因子组成的模块与3' UTR延长相关（模块4；图S7G）。
4. 尽管这些基因是翻译控制和RNA稳定性的良好研究调节因子，但之前未曾有报道与多聚A位点选择调节相关。
5. 我们还识别出一个包含大核糖体亚基组件以及翻译起始因子EIF6的模块（模块2；图S7E），这表明替代性多腺苷酸化和多个RNA调节过程之间存在紧密的交互作用。

Para_02

1. 我们发现多聚腺苷酸尾外泌体靶向（PAXT）复合物（模块3）的成员之间存在相互关联的扰动响应，该复合物还包含核帽结合复合物成员NCBP2和剪接调节因子MBNL1。
2. 对此模块的扰动主要与内含子多腺苷酸化转录本的上调有关（图7B、7D和S7C）。
3. 这种反应可能是由于PAXT复合物在降解过早终止的RNA转录物中的作用，尽管具体识别过早转录本的监控机制尚不清楚。

Para_03

1. 最后，我们旨在验证我们之前的假设，即与整合器扰动相关的内含子多A位点变化是由于剪接动态变化引起的（图3D和3I）。
2. 我们的CPA-Perturb-seq整合器签名与运动神经元生存（SMN）复合体的成员以及其他整合器复合体成员（包括一个额外的模块，模块5；图7E，7F和S7H）相关性最强。
3. 我们确认GWPS中受剪接模块调控的位点集合独特地与我们的整合器扰动响应特征重叠（图7G）。
4. 最后，我们观察到这些位点位于具有较弱的典型5′供体剪接位点得分的内含子中，这与剪接效率降低一致（图7H）。
5. 这为我们的假设提供了正交支持，即整合器扰动改变内含子中多A位点的使用，这些内含子的剪接效率低下，表明它们对剪接动态变化敏感。

Para_04

1. 我们得出结论，3'端scRNA-seq数据可以与定制的计算管线结合，用于探索多聚A位点使用的细胞异质性，并且我们开发了一个开源的R包，即使用相对转录丰度进行多聚A位点分析（PASTA），它实现了本文描述的分析方法。
2. PASTA与我们的分析工具包Seurat完全兼容，软件发布包括了演示如何用户可以探索循环人外周血单核细胞数据集中的替代多聚腺苷酸化细胞异质性的指南。
3. 这些数据和代码资源将促进在多种生物系统中异质替代多聚腺苷酸化的表征，并深入理解调控转录后调控的序列和调控因子。

Discussion

Para_01

1. 在这项研究中，我们证明了Perturb-seq技术，该技术已被广泛用于研究转录调控网络，也可以成功地应用于研究转录后调控。我们引入了一个统计框架，以在单细胞分辨率下量化调节因子之间的相对多聚A位点使用变化，并确定共同调控的多聚A位点模块。
2. 我们引入了一个统计框架，以在单细胞分辨率下量化调节因子之间的相对多聚A位点使用变化，并确定共同调控的多聚A位点模块。

Para_02

1. 我们的CPA-Perturb-seq数据集揭示了扰动反应的显著异质性，包括与每个调节因子相关的变化数量、类型和方向性。
2. 这表明替代性多腺苷酸化并不是统一调节的，即所有的多A位点对核心调节因子的扰动同等敏感。
3. 相反，我们一致观察到了多A位点模块间的不同调节反应的证据。

Para_03

1. 利用我们的深度神经网络APARENT-Perturb，我们发现这种局部调控结构部分由剪切位点周围的序列特异性元素编码。
2. 通过将先前训练的基于序列的模型与我们的扰动数据整合，我们直接学习序列元素与调控因子之间的关联，为顺式调控元件的功能提供了更机械的理解，包括调控因子之间的相互作用。
3. 我们注意到这种策略可以扩展到其他基于序列的深度学习模型。

Para_04

1. 尽管我们的分析旨在关注影响CPA决策的调控机制，但我们反复观察到额外的RNA调控过程改变了替代多腺苷酸化转录本的相对丰度。
2. 我们发现，影响RNA聚合酶延伸、RNA输出、翻译和剪接的蛋白质的扰动导致多腺苷酸化位点的差异使用，突显了RNA调控过程之间的广泛相互依赖性。
3. 未来的工作可能会利用这些相互依赖性，从3'端单细胞RNA测序数据中推断RNA动力学参数。
4. 更广泛地说，我们的统计方法可以扩展到描述转录组多样性的其他来源，例如选择性剪接。

Para_05

1. 我们的方法在测量和解释多腺苷酸化变化方面存在局限性。尽管调控因子之间的相关性扰动响应可能反映了共享功能，但我们不能排除间接效应（即一个调控因子的扰动影响另一个的表达）也可能影响相关性结构。
2. 扰动可能会改变单个转录本的生产和降解速率，而Perturb-seq无法明确区分这两种现象。
3. 此外，大幅度缩短多聚A尾长度的扰动可能会产生我们在本研究中无法解决的转录本捕获偏差。
4. 最后，我们对多腺苷酸化位点的注释基于polyA_DB数据库，我们可能偶尔会将串联多腺苷酸化位点误分类为内含子。
5. 未来的研究如果将CPA-Perturb-seq工作流程与RNA代谢标记101,102或长读测序103相结合，可以更准确地定量异构体，代表了我们工作的激动人心扩展。

Para_06

1. 展望未来，我们相信单细胞RNA测序（scRNA-seq）分析后的转录后调控，来自扰动筛选和原始样本的分析将相互提供信息。
2. 功能性基因组学工具如Perturb-seq非常适合用来识别分子调节器的靶标。
3. 我们设想，从建立因果关系的实验中推断出的分子特征代表着重要的资源，可用于解释那些因果关系不明的分子特征，如疾病状况。
4. 因此，这些数据集的整合代表了一条潜在的途径，可以系统地重建指导RNA生命周期的调控网络。

STAR★Methods

Key resources table

关键资源表

Resource availability

资源可用性

Lead contact

主要联系人

Para_01

1. 进一步的信息和资源、试剂的请求应直接联系负责人 Rahul Satija（rsatija@nygenome.org），他将负责满足这些需求。

Materials availability

材料可用性

Para_01

1. 本研究未生成独特的试剂。

Data and code availability

数据和代码的可获得性

Para_01

• 本次研究生成的CPA-Perturb-seq数据集可在https://zenodo.org/record/7619593#.Y-P7Zi1h2X0下载。所有原始数据均可在GEO数据库中找到，访问号为GEO: GSE269600。
• CPA-Perturb-seq 数据可以通过定制的 UCSC GenomeBrowser 轨道进行探索，该轨道可在 https://satijalab.org/cpa-perturb-seq 获得。
• Seurat 和 PASTA 均可作为开源 R 包在 https://github.com/satijalab/seurat 和 https://github.com/satijalab/PASTA 获得。
• 用于训练和解释APARENT-Perturb模型的代码可在https://github.com/johli/aparent-perturb找到。
• 本文报告的数据需要重新分析时所需的任何额外信息，可应要求从主要联系人处获得。