多任务学习用于多模态生物数据分析

目前的生物技术可以同时测量来自同一细胞的多种模态数据（例如RNA、DNA可及性和蛋白质）。这需要结合不同的分析任务（如多模态整合和跨模态分析）来全面理解这些数据，推断基因调控如何驱动生物多样性。然而，目前的分析方法被设计为执行单个任务，并且大部分仅提供多模态数据的部分图谱（比如单模态和双模态）。因此，作者提出了UnitedNet，这是一种可解释的多任务深度神经网络，能够整合不同的任务来分析单细胞多模态数据。应用于各种多模态数据集（如patch-seq、multiome ATAC+gene expression和空间转录组），与最先进的方法相比，UnitedNet在多模态整合和跨模态预测方面表现出类似或更好的准确性。此外，通过使用可解释的机器学习算法对经过训练的UnitedNet进行剖析，可以直接量化基因表达与其他具有细胞类型特异性的模式之间的关系。UnitedNet是一个全面的端到端框架，可广泛应用于单细胞多模态生物学。该框架有可能促进跨转录组学和其他模态的细胞类型特异性调节发现。

来自：Explainable multi-task learning for multimodality biological data analysis

背景概述

单细胞生物技术的最新进展使同时测量同一细胞的基因表达和其他模态成为可能。例如，patch-seq技术同时测量细胞基因表达和细胞内电活性（intracellular electrical activity），而multiome ATAC+gene expression技术联合测量细胞基因表达和DNA可及性。这种多模态组学数据同时提供了细胞的全面状态。然而，为分析单模态生物数据而开发的方法不能直接应用于多模态数据。与单模态分析相比，最近的研究确定了更多的多模态分析任务，例如：

• 从不同模态中识别具有生物学意义的簇，使人们能够更深入地了解不同生物系统的细胞特性和功能。
• 不同模态之间的跨模态预测，推断出不能容易或同时测量的模态数据。此外，为同一类型的细胞生成的多模态数据提供了发现基因表达和其他模态之间的细胞类型特异性关系的机会，有助于揭示感兴趣的生物学调控机制。

目前已经开发了几种多模态分析方法。对于joint group识别任务，已经开发了多模态数据整合方法，以将不同的模态融合到联合表示中，然后将其用于无监督或有监督分类，以识别细胞类型。对于跨模态预测任务，已经开发了基于自编码器的神经网络来预测不同模态。对于跨模态的相关性发现，Schema代表了最先进的多模态整合方法，可以识别用户定义的主模态中对其他模态重要的特征。

与上述方法相比，一种可以在统一的框架内处理所有任务、量化细胞类型特异性，跨模态相关性，并在没有先验知识的情况下，这样的方法可以简化数据分析，潜在提高每个任务的性能，并有助于从单细胞多模态数据中获得生物学见解。

尽管如此，由于以下两个原因，将多个任务组合到一个框架中可能具有挑战性：

• 首先，每个模态测量都具有独特的统计特征，需要不同的统计假设。虽然已经为不同的模态开发了几个统计模型，但仍然缺乏能够适应多模态未知分布的方法。
• 其次，joint group识别和跨模态预测通常代表单独的目标。具体而言，joint group识别的目的是惩罚错误的细胞群分配，而跨模态预测的目的是最小化重建数据和GT之间的差距。
• 此外，在没有先验知识的情况下，在特定细胞类型中寻找基因表达和其他模态之间的相关关系仍然是一个重大挑战。如果简单地迭代特征集的所有可能组合，对于高维数据来说，在计算上将是困难的。需要一种有效的方法来首先从对感兴趣的特定生物条件（例如细胞类型）中识别一组特征，然后量化这些特征之间的关系。

在这里，作者提出了一种可解释的多任务深度神经网络，以解决多模态数据分析的挑战。该网络具有编码器-解码器-判别器结构，并通过在两个任务之间交替进行训练：joint group识别和跨模态预测。具体而言，这种编码器-解码器-判别器结构并不假设数据分布是已知的，而是隐含地近似每个模态的统计特征。

在joint group识别和跨模态预测之间的交替训练保持并提高了这两项任务的性能。此外，应用可解释的机器学习来剖析网络，可以量化细胞类型特异性、跨模态特征与特征的相关性。作者已经将该网络应用于各种多模态数据集（图1a），包括：

• 带有GT标签的模拟多模态数据；
• 同时测量转录组和细胞内电活性（multi-sensing数据）；
• 同时分析转录组学和DNA可访问性（multi-omics数据）；
• 空间分辨转录组与蛋白质组（多模态spatial-omics数据）；

• 图1：利用联合网络对多模态生物数据进行多任务学习。
• a：具有代表性的多模态生物学数据示意图。
• b：共享隐空间的示意图。编码器可以将来自同一细胞的不同模态测量作为隐编码投影到共享空间。隐编码可以通过解码器投影回模态特定空间。在隐空间中，表示来自同一细胞的不同模态的隐编码可以被整合为用于joint group识别的单模态编码。
• c：显示了基于b中设计的共享潜在空间的联合网络结构示意图。
• d和e显示了联合网络在联合组识别和跨模态预测之间的多任务学习，用于分析多模态数据。
• f和g显示了可解释的应用，以剖析经过训练的UnitedNet，用于识别组与特征的相关性（f）和跨模态特征与特征的关联（g）。

方法

联合网络应用于联合测量数据。联合网络UnitedNet具有跨模态预测和joint group识别的特点，这两个任务是基于学习到的不同模态的联合低维表示完成的，其包含了细胞的潜在信息。假设有$V$种不同的模态。令$n$表示细胞数，$p^{(v)}$表示第$v$模态的特征维数，$X^{(v)}\in R^{n\times p^{(v)}}$为第$v$模态的数据，$x_i^{(v)}$为数据中的第$i$个细胞，对于从$v_1$到$v_2$的模态预测，联合网络从$x_i^{(v_1)}$的隐空间预测$x_i^{(v_2)}$。对于group识别，联合网络首先混合$x_i^{(1)},...,x_i^{(V)}$，然后，基于共享的隐空间，返回group类别 k i ∈ { 1 , . . . , K } k_{i}\in\left\{1,...,K\right\} ki​∈{1,...,K}。

UnitedNet由编码器、解码器、判别器和组识别模块组成。它是基于模态内预测损失、跨模态预测损失、生成器损失、鉴别器损失、对比损失、聚类损失（用于无监督的组识别）和分类损失（用于有监督的组识别）来训练的。关于联合网络的组成部分和损失的详细信息如下。

对于每个模态$v=1,...,V$，UnitedNet有一个编码器$En{c}^{(v)}(\cdot )$用于将细胞$i$对应的$x_i^{(v)}$映射到模态特定的隐空间编码$z_i^{(v)}$：$$z_i^{(v)}=En{c}^{(v)}(x_i^{(v)})$$其中，来自不同模态的低维表示需要具有相同数量的特征数。

解码器$De{c}^{(v)}(\cdot )$以模态特定的隐编码作为输入，并映射到模态$v$的特征上，以下表示模态$v_1$到$v_2$的预测：$${\tilde{x}}_i^{(v_1,v_2)}=De{c}^{(v_2)}(z_i^{(v_1)})$$对于同一个模态$v$下的预测，则为：$${\tilde{x}}_i^{(v)}=De{c}^{(v)}(z_i^{(v)})$$

判别器有助于训练编码器和解码器。模态$v$的判别器$Di{s}^{(v)}(\cdot )$以模态内预测${\tilde{x}}_i^v$或原始输入$x_i^v$为输入，输出二元分类结果，旨在区分${\tilde{x}}_i^v$和$x_i^v$。

组的数量为$K$，组识别模块从所有模态中获取模态特定的编码$(z_i^{(1)},...,z_i^{(V)})$，并将其分配给$K$个组中的一个，首先融合数据：$$z_i=\sum _{v=1}^V{\eta }_v{z}_i^{(v)}$$其中，${\eta }_1,..,{\eta }_V$为非负的可训练权重，并且满足${\sum }_{v=1}^V{\eta }_v=1$。然后，令$z_i$通过全连接层：$h_i=laye{r}_1(z_i)$。然后再进行softmax分配：$$\begin{gathered} {\alpha }_i=laye{r}_2(h_i)=(M_1(h_i),...,M_K(h_i){)}^T \\ {M}_k(h)=\frac{exp(W_{k}h)}{{\sum }_{t=1}^{K}exp(W_{t}h)} \end{gathered}$$其中，$W_k(k=1,...,K)$是系数向量。组识别模块将索引$k_i=argma{x}_{k=1,...,K}{M}_k(h_i)$为细胞$i$的簇类型。

聚类损失

聚类损失由三个部分组成。前两个组件来自基于深度散度聚类（DDC）。它们保证了所得到的组是可分离的和紧致的。第三个组成部分是基于自熵。

对于第一部分：组件1将不同组的聚类概率分配（软分配）之间的二乘二的相关性减小。它增加了组的可分性，因为这些相关性与不同组之间的相似性呈正相关。定义矩阵$S\in R^{n\times n}$，元素$s_{i,j}=exp(-||h_i-h_j|{|}_2^2/(2{\sigma }^2))$，其中$\sigma$是一个超参数。该矩阵测量了不同细胞之间的相似性。用${\tilde{\alpha }}_k$记$(M_k(h_1),...,M_k(h_n){)}^T$，这是细胞从$1$到$n$并分配为组$k$的概率。组件1损失为：$$L_{c1}={\left(\begin{array}{c}K\\ 2\end{array}\right)}^{-1}\sum _{k=1}^{K-1}\sum _{l>k}^K\frac{{\tilde{\alpha }}_k^{T}S{\tilde{\alpha }}_l}{\sqrt{{\tilde{\alpha }}_k^{T}S{\tilde{\alpha }}_k{\tilde{\alpha }}_l^{T}S{\tilde{\alpha }}_l}}$$其中，${\tilde{\alpha }}_k^{T}S{\tilde{\alpha }}_l$被最小化代表$M_k(h)$和$M_l(h)$趋近正交。

对于第二部分：组件2将不同组的软分配值推送到$R^K$中的单纯形的不同角，这也增加了组的可分性。令$e_k\in R^K$表示其第$k$个元素为1，其它元素为0的向量。因此，$e_k$是单纯形的第$k$个角。由于${\alpha }_i$是组识别模块中的第二层输出，令$m_k\in R^n$，其第$i$个元素为$exp(-||{\alpha }_i-e_k|{|}_2^2)$，其测量了软分配${\alpha }_i$和$e_k$之间的距离。组件2损失为：$$L_{c2}={\left(\begin{array}{c}K\\ 2\end{array}\right)}^{-1}\sum _{k=1}^{K-1}\sum _{l>k}^K\frac{m_k^{T}S{m}_l}{\sqrt{m_k^{T}S{m}_k{m}_l^{T}S{m}_l}}$$该损失加强了$exp(-||{\alpha }_i-e_k|{|}_2^2)$和$exp(-||{\alpha }_i-e_l|{|}_2^2)$的正交，因此，软分配输出$(M_1(h),...,M_K(h){)}^T$会倾向于一个单纯形角，而不是同时接近多个单纯形角。因此，同一组中的低维表示将是紧凑的，而来自不同组的低维表示将被分离。

对于第三部分：组件3旨在避免大多数细胞被分配给所有组中的一小部分。令$\overline{\alpha }\in R^K$表示${\alpha }_i$（$i=1,...,n$）求平均，其第$k$个元素为${\overline{\alpha }}_k$。组件3损失为：$$L_{c3}=\sum _{k=1}^K{\overline{\alpha }}_{k}log{\overline{\alpha }}_k$$其目的是以相等的概率分配每个组索引，比如${\overline{\alpha }}_1=,...,{\overline{\alpha }}_K=1/K$。

预测损失

模态内预测损失定义为：$$L_{Wpredict}=\frac{1}{nV}\sum _{i=1}^n\sum _{v=1}^V||{\tilde{x}}_i^{(v)}-x_i^{(v)}|{|}_2$$跨模态预测损失为：$$L_{Cpredict}=\frac{1}{n\left(\begin{array}{c}V\\ 2\end{array}\right)}\sum _{i=1}^n\sum _{v_1<v_2}||{\tilde{x}}_i^{(v_1,v_2)}-x_i^{(v_2)}|{|}_2$$其测量模态$v_1$到$v_2$的跨模态预测与模态$v_2$的原始值之间的距离。

生成损失和判别损失

生成器和判别器通过最小二乘损失训练，将模态内预测特征${\tilde{x}}_i^{(v)}$分配给标签1，将原始特征$x_i^{(v)}$分配给标签0。生成损失为：$$L_{Gen}=\frac{1}{nV}\sum _{i=1}^n\sum _{v=1}^V||Di{s}^{(v)}({\tilde{x}}_i^{(v)})-1|{|}_2^2$$判别损失为：$$L_{Dis}=\frac{1}{nV}\sum _{i=1}^n\sum _{v=1}^V||Di{s}^{(v)}({\tilde{x}}_i^{(v)})|{|}_2^2+\frac{1}{nV}\sum _{i=1}^n\sum _{v=1}^V||Di{s}^{(v)}(x_i^{(v)})-1|{|}_2^2$$其目的是使鉴别器将模态内预测特征分类为0，将原始特征分类为1。这组最小二乘损失提高了训练生成器的质量，因为它符合生成器的基本目标，即生成具有与原始特征相似分布的特征数据。

对比损失

应用对比损失来对齐来自不同模态的潜在编码，定义余弦相似度为：$$s_{i,j}^{(v_1,v_2)}=\frac{(z_i^{(v_1)}{)}^T{z}_j^{(v_2)}}{||z_i^{(v_1)}|{|}_2\cdot ||z_j^{(v_2)}|{|}_2}$$令：$$l_i^{(v_1,v_2)}=-log\frac{exp(\frac{s_{i,j}^{(v_1,v_2)}}{\tau })}{{\sum }_{s^{\prime }\in Neg(z_i^{(v_1)},z_i^{(v_2)})}exp(\frac{s^{\prime }}{\tau })}$$其中，$Neg(z_i^{(v_1)},z_i^{(v_2)})$通过从集合$N_i=\left\{s_{ij}^{(v_1,v_2)}:j=1,...,n,j\ne i,v_1,v_2=1,...,V,argmax{\alpha }_i\ne argmax{\alpha }_j\right\}$采样固定数量的元素获得，$\tau$为超参数。对比损失为：$$L_{Con}=\delta \cdot \frac{1}{n\left(\begin{array}{c}V\\ 2\end{array}\right)}\sum _{i=1}^n\sum _{v_1<v_2}{l}_i^{(v_1,v_2)}$$其中，$\delta$为超参数。

分类损失

定义向量$b_i$，其第$k$个元素为1，其余为0，其中，$k$是细胞$i$的观测类别标注。令$g_i=n_{k(i)}/n$，其中，$n_k$为类别为$k$的细胞数量。通过交叉熵评估分类准确率：$$L_{entropy}=-\frac{1}{n}\sum _{i=1}^n{g}_i\cdot (log(M_1(h_i)),...,log(M_K(h_i)){)}^T{b}_i$$

训练步骤

作者首先提出了在没有细胞标签的情况下训练UnitedNet用于跨模态预测和无监督组识别的过程。它由两个步骤迭代训练：组识别更新步骤和预测更新步骤。在组标识更新步骤中，编码器输出模态特定编码并将其送到组识别模块。然后，组识别模块将模态特定编码融合为共享潜在代码，并获得 $K$-dimensional 软聚类分配。解码器基于模态特定编码输出模态内预测特征${\tilde{x}}_i^{(v)}$。接下来，编码器和聚类模块通过聚类损失进行更新：$$L_{group}=L_{c1}+L_{c2}+L_{c3}$$在预测更新步骤中，编码器从不同模态输出低维表示，将它们送到解码器，并获得跨模态预测特征${\tilde{x}}_i^{(v_1,v_2)}$和模态内预测特征${\tilde{x}}_i^{(v)}$。然后，将预测的特征输入到鉴别器中。鉴别器由鉴别器损失$L_{Dis}$更新。接下来，通过模态内预测损失、跨模态预测损失、生成器损失和对比损失之和来更新编码器和解码器：$$L_{PGC}=L_{Wpredict}+L_{Cpredict}+L_{Gen}+L_{Con}$$以上两个训练步骤分别总结在算法1和2中。

• 算法1

• 算法2

对于有监督学习，为了训练UnitedNet进行监督分类，通过以下方式修改组识别损失：$$L_{group}=L_{entropy}$$

使用SHAP进行特征相关性分析

为了深入了解不同特征的重要性，UnitedNet应用了常用于解释机器学习模型的SHAP（SHapley Additive exPlanations）。这种方法的思想是在固定其他特征的同时，通过线性函数近似特征对输出的影响，并且该函数的系数对应于Shapley值。SHAP的优点包括：

• 基于理论的可解释性
• 广泛的应用范围
• 不需要扰动模型或数据的计算过程，这是许多其他特征重要性方法所要求的。

下面将首先介绍Shapley值的计算过程，然后解释它如何应用于UnitedNet。

假设想评估特征$x_j$对函数$f(x)$（$x=(x_1,...,x_Q{)}^T$， j ∈ { 1 , . . . , Q } j\in\left\{1,...,Q\right\} j∈{1,...,Q}）的重要性。令$F$为$x$的特征集，$S$为$F$的子集。$|F|$和$|S|$为集合元素数。Shapley值计算为： ϕ j ( x ) = ∑ S ⊆ F − { j } ∣ S ∣ ! ( ∣ F ∣ − ∣ S ∣ − 1 ) ! ∣ F ∣ ! [ f S ∪ { j } ( x S ∪ { j } ) − f S ( x S ) ] \phi_{j}(x)=\sum_{S\subseteq F-\left\{j\right\}}\frac{|S|!(|F|-|S|-1)!}{|F|!}[f_{S\cup\left\{j\right\}}(x_{S\cup\left\{j\right\}})-f_{S}(x_{S})] ϕj​(x)=S⊆F−{j}∑​∣F∣!∣S∣!(∣F∣−∣S∣−1)!​[fS∪{j}​(xS∪{j}​)−fS​(xS​)]其中， F − { j } F-\left\{j\right\} F−{j}表示从$F$中删除特征$j$。$f_S(x_S)$通过对$f(x_S,x_{F-S}^{(i)})$求样本均值得到，$x_{F-S}^{(i)}$为不在$S$中的特征的第$i$个观测值。Shapley值通过计算去除第$j$个特征后$f(x)$变化的加权平均值来衡量该特征的重要性。

对于神经网络$f(x)=f^{(1)}\odot f^{(2)}\cdot \cdot \cdot f^{(L)}$，令第$l$层的输出维度为$L^{(l)}$，第$l$层的第$q$个输入特征为$e_q^{(l)}$，其样本均值为${\overline{e}}_q^{(l)}$。可以用以下递归公式通过Deep SHAP以高效的方式估计模型的Shapley值：$${\varphi }_{q,r}(f^{(l-1)}\odot f^{(l)},e^{(l-1)})=(e_q^{(l-1)}-{\overline{e}}_q^{(l-1)})\cdot \sum _{r^{(l-1)}=1}^{L^{(l-1)}}{m}^{(l-1)}(r^{(l-1)},q)\cdot m^{(l)}(r,r^{(l-1)})$$其中，${\varphi }_{q,r}(f^{(l-1)}\odot f^{(l)},e^{(l-1)})$衡量了$e_q^{(l-1)}$相对于$f^{(l-1)}\odot f^{(l)}$的输出的第$r$个元素的重要性。$$m^{(l)}(r,r^{(l-1)})={\varphi }_{r^{(l-1)},r}(f^{(l)},e^{(l)})/(e_{r^{(l-1)}}^{(l)}-{\overline{e}}_{r^{(l-1)}}^{(l)})$$$$m^{(l-1)}(r^{(l-1)},q)={\varphi }_{q,r^{(l-1)}}(f^{(l-1)},e^{(l-1)})/(e_q^{(l-1)}-{\overline{e}}_q^{(l-1)})$$为了识别与特定组具有高度相关性的特征，对于每个细胞，需要计算每个输入特征相对于该组的软分配的Shapley值。然后，将来自被分类为该组的细胞的所有Sharpley值作为绝对值，并用于计算平均值。具有最高平均值的前$n$个特征被解释为与该组具有高度相关性（Patch-seq GABAergic数据集的$n＝7$，multiome ATAC+Gene expression BMMCs数据集的$n＝20$）。

为了量化组内跨模态特征与特征的相关性，作者考虑了在前一步中选择的该组的高相关性特征。接下来，计算每个特征相对于来自另一模态的每个其他特征的Shapley值。然后，通过平均来聚集来自不同细胞的Sharpley值的绝对值。值相对较大的特征被认为是重要的。

评价指标

无监督和有监督组识别用ARI（adjusted rand index），预测性能用 $R^2$（coefficient of determination）和AUC（area under the ROC curve），两个任务之间的关系用Pearson’s correlation评估。

对于ARI：
将从模型中获得的簇与从细胞类型标签中获得的簇进行比较。令$a_{k_1}$（$k_1=1,...,K$）表示来自模型的第$k_1$个簇中的细胞数，$b_{k_2}$（$k_1=1,...,K$）表示来自细胞类型标签的第$k_2$个簇的细胞数。$n_{k_1,k_2}$表示来自模型的第$k_1$个簇和来自细胞类型标签的第$k_2$个簇中的观测次数。当来自模型的聚类结果接近于来自观察到的细胞类型标签的聚类结果时，ARI接近于1，并且对于随机猜测接近于0。

对于$R^2$：
令$y$和$\tilde{y}$表示观测数据和预测数据，$\overline{y}$被定义为与$y$长度相同的向量，并且每个元素都是$y$的样本均值。$R^2$为：$$R^2=1-\frac{||\tilde{y}-y|{|}_2^2}{||\overline{y}-y|{|}_2^2}$$它将来自模型的预测$\tilde{y}$的均方误差（MSE）与以常数值$\overline{y}$作为预测的基线的MSE进行比较。值在负无穷到1之间，当预测等于观测时，$R^2=1$。

对于Pearson’s correlation：
令$\hat{y}$表示与$\tilde{y}$等长度的向量，其每个元素为$\tilde{y}$的样本均值。因此，$y$与$\tilde{y}$之间的Pearson’s correlation为：$$r=\frac{<\tilde{y}-\hat{y},y-\overline{y}>}{||\tilde{y}-\hat{y}|{|}_2\cdot ||y-\overline{y}|{|}_2}$$当预测与观测数据分别具有正或负线性关系时，它取[-1,1]中的值，比如等于1或负1。当预测和观测数据没有线性关系时，它等于零。

对于AUC：
为了对multiome ATAC+Gene expression数据集进行建模，作者对DNA可及性数据进行了二值化。因此，为了评估这些数据的预测，我们采用了ROC曲线下的面积。设$n_0$和$n_1$分别表示来自观测数据的0和1的数目。$p_{0,i}$（$i=1,..,n_0$）和$p_{1,j}$（$j=1,..,n_1$）分别表示两组观测值的模型预测。AUC针对于二分类问题，AUC取0到1之间的值，1对应于完美的预测。

数据集

多任务数据集

Dyngen：
作者用Dyngen来模拟4-modality数据集。具体而言，生成500个具有模拟DNA、pre-mRNA、mRNA和蛋白质模态的细胞，每个模态包含100个维度的特征。同时，Groundtruth细胞类型注释与数据集一起生成。对于Dyngen模拟器的参数，使用Dyngen教程中线性主干模型的默认设置，函数包括backline_linear、initialize_mode和generate_dataset。

MUSE：
作者将MUSE中的模拟器应用于模拟2-modality输入，以评估具有一个低质量模态的UnitedNet的稳健性。作者模拟了11个具有1000个细胞和10种细胞类型的双模态数据集。每个模态包含500个模态特定特征。对于11个数据集中的每一个，用可控的衰减系数模拟其中一个。当与其他方法进行基准测试时，使用0.01、0.1、0.2、0.3、0.4、0.5、0.6、0.7、0.8、0.9和1作为不同的衰减系数。

Patch-seq GAVAergic neuron 数据集：
作者使用Patch-seq数据集，该数据集同时表征了从小鼠视觉皮层中GABAergic interneurons获得的形态学（M）、电生理学（E）和转录组学（T）特征。在进行质量控制后，使用了相同的数据集，其中3395个神经元用于E-T分析，448个神经元用于M-E-T分析。作者对每个模态的输入矩阵进行标准化，使每个细胞中所有特征的平均值和标准差分别为0和1。

Multiome ATAC + gene expression BMMCs数据集：
作者使用了一个多组学ATAC+基因表达数据集，该数据集同时结合了从10个受试者和4个组织位点在BMMC（骨髓单核细胞）组织中获得的基因表达和全基因组DNA可及性。除了先前研究中的质量控制外，还使用了多组ATAC+基因表达BMMC数据集的标准预处理程序。对于基因表达模态的预处理，作者使用中值标准化和log1p变换和标准化，并通过Scanpy选择前4000个高变基因。对于DNA可及性模态的预处理，通过将所有非零值替换为1来对数据进行二值化，并通过Scanpy选择前13634个高变的DNA可及性特征。作者使用ChIPseeker和scanpy.var_names_make_unique来注释DNA可及性峰值。

空间组学数据

生成niche表达模态：
使用测量的细胞或spots的RNA表达，作者整合了每个细胞或spots的空间信息，并生成RNA的加权平均表达。在二维空间坐标$(s_i^1,s_i^2)$和模态$v$的第$i$行$x_i^{(v)}$对应于细胞或spots $i$ 的情况下，可以计算模态$v$的niche模态，用$x^{(vniche)}$（$v=1,..,V$）。对于细胞$i$，$x_i^{(vniche)}$为：$$x_i^{(vniche)}=\sum _{j=1}^J{x}_j^{(v)}\cdot w_{ij}$$其中， j ∈ { 1 , . . . , J } j\in\left\{1,...,J\right\} j∈{1,...,J}表示属于cell或spot $i$ 的$J$-nearest的cells或spots，$w_{ij}$为：$$w_{ij}=\frac{1/distance\left\{(s_i^1,s_i^2),(s_j^1,s_j^2)\right\}}{{\sum }_{j=1}^{J}1/distance\left\{(s_i^1,s_i^2),(s_j^1,s_j^2)\right\}}$$其中$distance$表示两个向量之间的欧几里得距离。

DBiT-seq embryo数据集：
作者使用DBiT-seq embryo（胚胎）数据集，其中采用了DBiT-seq 936个spots的以下三种模态：mRNA表达、蛋白质表达和niche mRNA表达。

• 对于mRNA表达模态，使用scanpy的函数scanpy.pp.normalize_total对原始计数矩阵进行归一化，并选择前568个差异表达基因。
• 对于蛋白质表达，对原始计数矩阵进行了归一化，并使用了22种蛋白质。
• niche模态是基于标准化的mRNA表达产生的。

对于组织区域表征的第一项任务，从原始研究中提取真实组织区域标签，这是基于H&E图像的主要组织区域的解剖学注释。作者将UnitedNet的聚类结果与其他最先进方法的聚类结果进行了比较。通过ARI来验证它们的性能。对于跨模态预测的并行任务，尽管有三种模态被用作UnitedNet模型的输入，但作者专注于第一种和第二种模态之间的预测：mRNA表达和蛋白质表达。由于DBiT-seq公共数据集中只有一个批次，作者将DBiT-seq embryo数据集中的936个spots分为用于预测任务的训练数据集（80%，748个spots）和测试数据集（20%，188个spots）。

DLPFC数据集：
作者使用了12个批次的成人背外侧前额叶皮层（DLPFC）数据集。使用以下三种模式：mRNA表达、从H&E染色图像中提取的形态学特征和 niche mRNA。对于mRNA表达的模态，对原始计数矩阵进行归一化，并选择前2365个差异表达基因。使用预先训练的卷积神经网络从stLearn实现的H&E染色图像中提取形态学特征。50维形态特征被用作每个spots的第二模态。对于有监督组识别任务，使用11个批次及其组织区域注释来训练UnitedNet模型。然后，将训练后的模型应用于剩余批次，以识别组织区域注释，并在H&E图像特征和mRNA表达之间进行跨模态预测。作者将UnitedNet的识别性能与来自原始DLPFC论文的SpatialDE PCA和 pseudobulk PCA进行了比较。在识别任务之后，作者对最近邻居中具有35个spots的SpaGCN之后的聚类结果应用了细化步骤。

结果

特征与特征相关性

作者使用事后可解释学习（SHAP）对训练的UnitedNet进行剖析，以表明Patch-seq GABAergic数据集中的特征相关性。具体而言，使用SHAP将重要性值（称为Shapley值）分配给关于任何给定模型输出的每个输入特征，例如特定识别的细胞组或特定特征的跨模态预测。根据定义，具有高Shapley值的特征是有影响力的。因此，可以根据特征对Shapley值的排名来选择特征。

以Pvalb神经元类型（一种细胞大类）为例，作者定性地验证了SHAP选择的相关性（图2）。对于组与特征的相关性，SHAP成功地选择了Pvalb神经元差异表达的基因、电生理特征和形态学特征的子集（图2a，d–f）。

对于Pvalb神经元特异性跨模态特征与特征相关性（图2b，c），发现在使用长方电流步长的膜片钳电刺激过程中（long square current steps），基因Lrrc38与Pvalb神经细胞平均放电率的电生理特征（average firing rate）表现出更高的相关性。这一结果与之前的研究一致，表明Lrrc38相关蛋白是大钾（BK）通道最关键的调节剂之一，大钾通道对神经元放电动力学和神经递质释放至关重要。这些结果表明，UnitedNet可以潜在地用于促进Patch-seq数据的细胞类型特异性基因功能相关性的鉴定。

• 图2：特征相关性解释（Patch-seq GABAergic数据集）。

然后，作者探索了UnitedNet中可解释学习对ATAC+基因表达数据集的生物学价值。以CD8+T主要细胞类型（CD8+T和CD8+T naive细胞）为例（图3a），结果首先鉴定了CD8+T主细胞类型差异表达的基因亚群（例如，CD8A、A2M、LEF1和NELL2）和DNA可及性位点（例如，D8A、DPP8、KDM2B和KDM6B）（图3b、d、e）。在这些基因和DNA可及位点中，DNA可及性位点PROS1、KDM2B和KDM6B与CD8+T细胞特异性基因表现出更强的相关性，表明它们在CD8+T功能中的关键作用（图3c）。该结果与先前的研究一致，即CD8+T细胞中PROS1表达水平的升高是防止免疫反应过度活跃的关键调节信号。同时，KDM2B表达的缺乏启动了T细胞白血病的发生。值得注意的是，最近的一项研究发现，与KDM2B属于同一基因家族的KDM6B通过诱导效应相关基因中的DNA可及性直接调节CD8+T细胞的产生。结果进一步表明，KDM2B也可能在调节CD8+T细胞的产生中发挥重要作用（这是模型发现的，这一点是重要的）。

• 图3：特征相关性解释（ATAC+Gene表达数据集）。

空间组学上的应用

空间组学是一种重要的模态技术，可以测量完整组织中的空间分辨多组学信息。然而，在分析用于组识别任务的空间组学数据时，空间信息往往没有被充分利用。联合网络可以灵活地整合不同的模态作为输入，包括空间信息。UnitedNet可以利用细胞niche信息（每个细胞的邻域基因表达信息）作为额外的模态来识别具有生物学意义的群体，并增强跨模态预测（图4a）。

作者首先将UnitedNet应用于单批次的DBiT-seq embryo数据集，该数据集同时映射了胚胎组织上的整个转录组和22种蛋白质。具体而言，作者生成了编码空间信息的细胞niche信息（RNA表达的加权平均值），作为第三种分析模态。然后，联合网络将基因表达、蛋白质和niche模态相结合，用于组织区域的无监督联合识别以及基因表达和蛋白质之间的跨模态预测。通过将原始报告中的组织区域的解剖注释视为Ground Truth，对组织区域识别的准确性进行了基准测试。与最先进的方法相比，UnitedNet实现了更高的无监督组识别精度（图4b）。此外，联合网络使几个代表性基因和蛋白质表达之间的空间分辨跨模态预测成为可能（图第4c）。

• 图4：无监督组识别，joint注释转移，空间组学跨模态预测。联合网络的空间组学数据分析管道示意图。空间组学同时测量完整组织网络中的空间分辨多组学数据。UnitedNet提取细胞邻域信息作为附加模态以及用于无监督或有监督群组识别和跨模态预测的其他模态。

接下来，作者将UnitedNet应用于一个带注释的多批次空间组学数据集，用于同时有监督的联合组识别和跨模态预测。作者使用了人类背外侧前额叶皮层（DLPFC）数据集，该数据集在空间上绘制了12批次DLPFC脑切片的基因表达和H&E染色。同样，使用基因表达、基于H&E染色的形态学特征和细胞niche模态作为UnitedNet的输入。UnitedNet可以成功地注释看不见的DLPFC切片，并实现比其他基准测试方法和未使用交替训练方案或细胞小众模式的消融版UnitedNet更高或相当的精度（图4d）。

此外，作者探讨了UnitedNet是否可以减少空间DLPFC数据集中的批次效应，其方式与ATAC+基因表达BMMC数据集的分析类似。结果表明，与其他消融研究相比，UnitedNet在潜在空间中保持了良好的可分性和减少批次效应的能力，使其在组识别任务中具有更高或可比的性能。此外，UnitedNet使几个代表性基因和H&E形态特征之间的跨模态预测成为可能（图4e）。