UnpairReg：回归分析用于不配对多组学分析

尽管最近取得了进展，但很难对同一细胞上的所有多组学单细胞数据模式进行分析。因此，产生了大量的单细胞不配对组学数据。作者提出了一种名为UnpairReg的方法，用于对未配对观测值进行回归分析，以整合单细胞多组学数据。在真实和模拟数据上，UnpairReg提供了细胞基因表达的准确估计（在只有染色质可及性数据可用的情况下）。从UnpairReg推断的顺调控网络与eQTL mapping高度一致。此外，UnpairReg通过联合分析单细胞基因表达和染色质可及性数据来提高细胞类型鉴定的准确性。

来自：Integration of single‑cell multi‑omics data by regression analysis on unpaired observations

背景概述

scRNA-seq提供全基因组转录分析，scATAC-seq鉴定可访问的染色质区域，single cell bisulfite sequencing测量DNA甲基化，single cell CUT&Tag在单细胞水平上分析组蛋白修饰或转录因子。目前已经开发了几种单细胞多模态测序技术，如用于基因表达和蛋白质表达联合分析的CITE-seq，用于基因表达和染色质可及性联合分析的MultiOme，以及用于基因表达和组蛋白修饰联合分析的Paired-Tag。然而，在同一时间观察同一单细胞的所有图谱是高成本的。一种替代方法是生成其他模态，或者整合不配对的其他模态。

据报道，多组学分析通过整合不同类型的分子数据，提供了对细胞的全面理解。目前已经产生了大量未配对的scRNA-seq和scATAC-seq。scRNA-seq与scATAC-seq的整合分析可以更准确地识别亚群，并提供更详细的基因调控信息。为了联合分析这两类数据，很多方法都需要顺式调控元件(cis-regulatory elements, REs)和靶基因(target genes, TGs)之间的连接。例如，Coupled NMF和DC3，分别从来自不同细胞背景的外部bulk数据和bulk 3D染色质可及性数据中学习RE-TG连接。

基于距离的连接是有问题的，因为一些REs不调节距离近的基因。适当的外部bulk数据并不总是可用的，从外部bulk数据中获得的关联将忽略那些未包含在数据库中的细胞群特有的RE-TG关系。连接的更一般的形式是使用染色质可及性来预测基因表达。包括Seurat和Signac在内的许多方法计算scATAC-seq细胞的基因活性评分。Cicero定义了一种不同版本的基因活性评分，它是附近REs的加权和，其权重取决于RE与启动子可及性的相关性；MAESTRO将基因活性评分定义为附近REs的加权和，其中权重是RE-TG距离的指数递减函数；SnapATAC通过平滑基因体中的计数来定义基因可及性评分。准确的连接，包括RE-TG连接或从REs预测TG表达，是这些方法是否有效的关键。迫切需要一种不使用外部数据来学习连接的统计方法。

主要贡献：通过引入回归模型来预测非配对细胞上染色质可及性的基因表达，从而填补空白。该模型允许对特征和响应变量进行非配对计算，例如REs的可及性是特征，而靶基因的表达是响应变量。传统的回归模型需要对同一样本上的特征变量和响应变量进行观察(我们称之为配对观测)。

作者将回归问题转化为协方差水平的二次方程，其中响应变量的内积表示为特征变量的内积的二次方程。通过拟合这些方程，我们可以了解每个REs的系数，并预测仅测量染色质可及性的细胞的基因表达。从该模型中了解到的顺式调控有助于理解细胞类型特异性调控机制。模型预测的基因表达是REs附近基因可及性的加权和，其中权重是从统计模型中获得的系数。因此，这种方法应该比以前的方法提供更准确的基因表达估计。对基因表达的准确预测将提高细胞类型鉴定的能力，并有助于联合分析。

结果

未配对观测上的回归分析

作者提出了一种统计方法，用于综合分析来自相同组织/背景但来自不同细胞的未配对单细胞基因表达和染色质可及性数据。方法是基于非配对观测的回归模型，因此命名为UnpairReg。图1显示了UnpairReg。目标是：

• 预测只有染色质可及性数据可用的细胞的基因表达
• 了解调控元件(REs)和靶基因(TGs)之间的顺式调控关系

在图1中，$O$是$n_1\times p_1$的矩阵，表示$n_1$个细胞在$p_1$个REs上的染色质可及性。TGs的表达记录在$E$矩阵中，其中$E$是$p_2$基因和$n_2$细胞的基因表达$n_2\times p_2$矩阵。

如果我们在同一细胞上测量了基因表达和染色质可及性数据(成对数据)，我们可以通过拟合回归模型$E=O\beta +\epsilon$来了解RE-TG的关联，其中TG的表达表示为接近该基因的REs可及性的加权和，$\beta$表示要估计的系数(图1A)。然而，这种类型的回归分析对于未配对的数据是不可行的，其中基因表达和染色质可及性不是在同一细胞上测量的(图1B)。

• 图1A：基于配对染色质可及性和基因表达数据的线性回归模型示意图。在同一细胞中观测基因表达水平和染色质可及性。TG1和TG2来自两个不同的基因组位置。通过数据观察到RE1-TG1调控、RE5-TG2调控、TG1-TG2相关、RE1-RE5相关。
• 图1B：基于未配对数据的UnpairReg模型示意图。基因表达水平仅在细胞1至细胞5中观察到，而染色质可及性在细胞6至细胞10中观察到。在此数据中观察到TG1-TG2相关和RE1-RE5相关，但REs与TGs之间没有其他显著相关。因此，推断RE1-TG1调控和RE5-TG2调控。

通过图1中一个直观的例子来介绍UnpairReg的主要思想。假设我们有两个强相关的基因，TG1和TG2。每一个基因都有三个相邻的REs，分别为{RE1, RE2, RE3}和{RE4, RE5, RE6}。对于染色质可及性数据，如果我们观察到RE1和RE5有很强的相关性，而其他对不相关，那么很容易推断TG1的调控因子更可能是RE1而不是{RE2, RE3}，TG2的调控因子更可能是RE5而不是{RE4, RE6}。

在未配对的数据中，对于基因表达可用的细胞(我们称之为RNA-seq细胞)，观察不到REs的染色质可及性。因此，我们无法通过线性回归模型了解RE-TG的相关性(图1B)。只能利用其余的信息，RE1-RE5和TG1-TG2的相关性，以及基因组定位信息，这些对于确定RE对TG的调控非常重要。

UnpairReg从协方差矩阵中获取RE-RE和TG-TG相关性信息，再推断出顺式调控信息，并预测只有染色质可及性数据可用的细胞(ATAC-seq细胞)的基因表达。在数学上，基因协方差矩阵是RE-RE协方差矩阵的二次。

在一定的假设条件下，将线性回归问题转化为协方差矩阵回归问题：$$E^{T}E={\beta }^T{O}^{T}O\beta +{\epsilon }^T\epsilon$$其中，$E^{T}E$表示gene-gene协方差矩阵，$O^{T}O$表示RE-RE协方差矩阵。这种协方差回归允许我们获得与线性回归相似的系数。估计的回归系数反映了顺式调控，它们也可以用来预测只有染色质可及性数据可用的细胞的基因表达水平。

模拟数据上的表现

为了说明UnpairReg在顺式调控推断和基因表达预测中的作用，作者模拟了不同dropout下的顺式调控系数、基因表达和染色质可及性数据。

首先，作者评估了系数估计的有效性。以未配对数据为输入，通过UnpairReg估计顺式调控系数。由于我们知道该模拟数据中的真实顺调节系数，因此我们通过计算Pearson相关系数(PCC)将估计的系数与GT进行比较，以评估系数估计。图2A显示，在0.6到0.9的dropout下，PCC的范围从0.1到0.4，这表明UnpairReg准确地估计了系数。

为了评估基因表达预测的性能，作者将真实的基因表达数据作为GT，并计算预测与GT之间的PCC。在这里，计算跨基因的每个细胞的PCC，称为细胞水平PCC，以及跨细胞的每个基因的基因水平PCC。以真实的基因表达数据为基础，将UnpairReg与观察到的基因表达数据进行比较。

请注意，真实的基因表达数据与观察到的基因表达数据是不同的。对于未配对的数据，它们都没有被观察到。前者反映了每个单细胞中基因的实际表达水平，并且无法通过测序技术观察到。后者可以通过单细胞多组数据(配对数据)来测量，它提供了基因表达的粗略估计，因为它受到dropout的影响。图2B显示了我们方法的平均细胞水平PCC和观察到的基因表达数据。观察到的基因表达的平均细胞水平PCC随着dropout的增加从0.23下降到0.11，而预测基因水平的PCC从0.89下降到0.78。

• 图2：A.UnpairReg预测了不同dropout下的$\beta$和实际调控系数的相似性。B.预测基因表达水平与实际基因表达水平之间的PCCs。细胞水平的PCC是跨基因计算每个细胞的PCC。B至E中的颜色表示预测基因表达数据的不同类型。橙色表示UnpairReg预测的基因表达，蓝色表示dropout后观察到的基因表达，灰色表示随机数据。C.基因水平PCCs在预测基因表达水平和基本真实值之间的平均值。基因水平PCC是计算每个基因跨细胞的PCC。D.细胞水平PCC, dropout率为0.87。E.基因水平PCC在dropout为0.87的情况下。

MultiOme上的表现

为了评估UnpairReg在实际单细胞数据上预测基因表达的准确性，UnpairReg被应用于来自10X Genomics的外周血单核细胞(PBMC)。将UnpairReg与另一种名为基因活性评分(GAS)的基因表达预测方法进行比较。GAS是Seurat中用于整合scRNA-seq和scATAC-seq的方法。

• 图3：基因表达预测与配对数据一致。C.y轴和x轴表示细胞中预测的基因表达和log10(1+计数)。

首先通过将基因表达计数矩阵作为基础真值来评估每个细胞基因表达预测的准确性。图3A显示了所有细胞的细胞级PCCs，其中y轴代表UnpairReg，x轴分别代表PBMC和HHBT数据的GAS。这些结果表明UnpairReg对几乎所有细胞都有更好的表现。对于PBMC数据，细胞水平的平均PCC为0.55，比GAS的平均PCC 0.11大5倍，对于HHBT数据，UnpairReg和GAS的平均PCCs分别为0.39和0.19。

然后，比较了所有基因的PCC和GAS的基因水平。PBMC的平均基因水平PCC为0.15，是GAS的4倍(图3B)。此外考虑到GT受到dropout的影响，无法反映真实的基因表达水平。为了对抗dropout，作者对基因表达数据进行估算，然后将预测结果与估算的基因表达数据进行比较。结果得到了进一步提高。为了详细研究预测结果，作者选择一个细胞，并显示观察到的与预测的表达模式(图3C)。

除了通过相关性量化每个细胞基因表达预测的准确性外，我们还通过比对误差进一步量化它，比对误差定义为比真正匹配的距离更近的细胞数量（图3D）。

为了评估推断顺式调控的能力，作者计算了顺式调控系数与eQTL的一致性。我们从GTEx和eQTLGen中下载全血中由eQTL定义的变异基因链接，并使用它们验证RE-TG的预测。由于RE和TG之间的距离对预测很重要，作者根据RE-TG对的距离(0-3 kb、3-10 kb、10-20 kb、20-50 kb、50-100 kb和100-150 kb)将它们分为不同的组。在每个距离组中，以eQTL数据为基础真实值，通过滑动顺式调节得分计算AUPR。将UnpairReg与基于距离的方法、增强子和启动子可及性之间的PCC以及随机预测进行了比较（图3E）。