吉安哪家网站建设公司好,网站开发php程序员,缪斯设计,重庆欧勒精细陶瓷有限公司网站策划书目前的生物技术可以同时测量来自同一细胞的多种模态数据#xff08;例如RNA、DNA可及性和蛋白质#xff09;。这需要结合不同的分析任务#xff08;如多模态整合和跨模态分析#xff09;来全面理解这些数据#xff0c;推断基因调控如何驱动生物多样性。然而#xff0c;目…目前的生物技术可以同时测量来自同一细胞的多种模态数据例如RNA、DNA可及性和蛋白质。这需要结合不同的分析任务如多模态整合和跨模态分析来全面理解这些数据推断基因调控如何驱动生物多样性。然而目前的分析方法被设计为执行单个任务并且大部分仅提供多模态数据的部分图谱比如单模态和双模态。因此作者提出了UnitedNet这是一种可解释的多任务深度神经网络能够整合不同的任务来分析单细胞多模态数据。应用于各种多模态数据集如patch-seq、multiome ATACgene expression和空间转录组与最先进的方法相比UnitedNet在多模态整合和跨模态预测方面表现出类似或更好的准确性。此外通过使用可解释的机器学习算法对经过训练的UnitedNet进行剖析可以直接量化基因表达与其他具有细胞类型特异性的模式之间的关系。UnitedNet是一个全面的端到端框架可广泛应用于单细胞多模态生物学。该框架有可能促进跨转录组学和其他模态的细胞类型特异性调节发现。
来自Explainable multi-task learning for multimodality biological data analysis 目录 背景概述方法聚类损失预测损失生成损失和判别损失对比损失分类损失训练步骤使用SHAP进行特征相关性分析 评价指标数据集多任务数据集空间组学数据 结果特征与特征相关性空间组学上的应用 背景概述
单细胞生物技术的最新进展使同时测量同一细胞的基因表达和其他模态成为可能。例如patch-seq技术同时测量细胞基因表达和细胞内电活性intracellular electrical activity而multiome ATACgene 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 V V种不同的模态。令 n n n表示细胞数 p ( v ) p^{(v)} p(v)表示第 v v v模态的特征维数 X ( v ) ∈ R n × p ( v ) X^{(v)}\in R^{n\times p^{(v)}} X(v)∈Rn×p(v)为第 v v v模态的数据 x i ( v ) x_{i}^{(v)} xi(v)为数据中的第 i i i个细胞对于从 v 1 v_1 v1到 v 2 v_2 v2的模态预测联合网络从 x i ( v 1 ) x_{i}^{(v_1)} xi(v1)的隐空间预测 x i ( v 2 ) x_{i}^{(v_2)} xi(v2)。对于group识别联合网络首先混合 x i ( 1 ) , . . . , x i ( V ) x_{i}^{(1)},...,x_{i}^{(V)} xi(1),...,xi(V)然后基于共享的隐空间返回group类别 k i ∈ { 1 , . . . , K } k_{i}\in\left\{1,...,K\right\} ki∈{1,...,K}。
UnitedNet由编码器、解码器、判别器和组识别模块组成。它是基于模态内预测损失、跨模态预测损失、生成器损失、鉴别器损失、对比损失、聚类损失用于无监督的组识别和分类损失用于有监督的组识别来训练的。关于联合网络的组成部分和损失的详细信息如下。
对于每个模态 v 1 , . . . , V v1,...,V v1,...,VUnitedNet有一个编码器 E n c ( v ) ( ⋅ ) Enc^{(v)}(\cdot) Enc(v)(⋅)用于将细胞 i i i对应的 x i ( v ) x_{i}^{(v)} xi(v)映射到模态特定的隐空间编码 z i ( v ) z_{i}^{(v)} zi(v) z i ( v ) E n c ( v ) ( x i ( v ) ) z_{i}^{(v)}Enc^{(v)}(x_{i}^{(v)}) zi(v)Enc(v)(xi(v))其中来自不同模态的低维表示需要具有相同数量的特征数。
解码器 D e c ( v ) ( ⋅ ) Dec^{(v)}(\cdot) Dec(v)(⋅)以模态特定的隐编码作为输入并映射到模态 v v v的特征上以下表示模态 v 1 v_{1} v1到 v 2 v_{2} v2的预测 x ~ i ( v 1 , v 2 ) D e c ( v 2 ) ( z i ( v 1 ) ) \widetilde{x}_{i}^{(v_1,v_2)}Dec^{(v_2)}(z_{i}^{(v_1)}) x i(v1,v2)Dec(v2)(zi(v1))对于同一个模态 v v v下的预测则为 x ~ i ( v ) D e c ( v ) ( z i ( v ) ) \widetilde{x}_{i}^{(v)}Dec^{(v)}(z_{i}^{(v)}) x i(v)Dec(v)(zi(v))
判别器有助于训练编码器和解码器。模态 v v v的判别器 D i s ( v ) ( ⋅ ) Dis^{(v)}(\cdot) Dis(v)(⋅)以模态内预测 x ~ i v \widetilde{x}_{i}^{v} x iv或原始输入 x i v x_{i}^{v} xiv为输入输出二元分类结果旨在区分 x ~ i v \widetilde{x}_{i}^{v} x iv和 x i v x_{i}^{v} xiv。
组的数量为 K K K组识别模块从所有模态中获取模态特定的编码 ( z i ( 1 ) , . . . , z i ( V ) ) (z_{i}^{(1)},...,z_{i}^{(V)}) (zi(1),...,zi(V))并将其分配给 K K K个组中的一个首先融合数据 z i ∑ v 1 V η v z i ( v ) z_{i}\sum_{v1}^{V}\eta_{v}z_{i}^{(v)} ziv1∑Vηvzi(v)其中 η 1 , . . , η V \eta_{1},..,\eta_{V} η1,..,ηV为非负的可训练权重并且满足 ∑ v 1 V η v 1 \sum_{v1}^{V}\eta_{v}1 ∑v1Vηv1。然后令 z i z_{i} zi通过全连接层 h i l a y e r 1 ( z i ) h_{i}layer_{1}(z_{i}) hilayer1(zi)。然后再进行softmax分配 α i l a y e r 2 ( h i ) ( M 1 ( h i ) , . . . , M K ( h i ) ) T M k ( h ) e x p ( W k h ) ∑ t 1 K e x p ( W t h ) \alpha_{i}layer_{2}(h_{i})(M_{1}(h_{i}),...,M_{K}(h_{i}))^{T}\\ M_{k}(h)\frac{exp(W_{k}h)}{\sum_{t1}^{K}exp(W_{t}h)} αilayer2(hi)(M1(hi),...,MK(hi))TMk(h)∑t1Kexp(Wth)exp(Wkh)其中 W k ( k 1 , . . . , K ) W_{k}(k1,...,K) Wk(k1,...,K)是系数向量。组识别模块将索引 k i a r g m a x k 1 , . . . , K M k ( h i ) k_{i}argmax_{k1,...,K}M_{k}(h_{i}) kiargmaxk1,...,KMk(hi)为细胞 i i i的簇类型。
聚类损失
聚类损失由三个部分组成。前两个组件来自基于深度散度聚类DDC。它们保证了所得到的组是可分离的和紧致的。第三个组成部分是基于自熵。
对于第一部分组件1将不同组的聚类概率分配软分配之间的二乘二的相关性减小。它增加了组的可分性因为这些相关性与不同组之间的相似性呈正相关。定义矩阵 S ∈ R n × n S\in R^{n\times n} S∈Rn×n元素 s i , j e x p ( − ∣ ∣ h i − h j ∣ ∣ 2 2 / ( 2 σ 2 ) ) s_{i,j}exp(-||h_{i}-h_{j}||_{2}^{2}/(2\sigma^{2})) si,jexp(−∣∣hi−hj∣∣22/(2σ2))其中 σ \sigma σ是一个超参数。该矩阵测量了不同细胞之间的相似性。用 α ~ k \widetilde{\alpha}_{k} α k记 ( M k ( h 1 ) , . . . , M k ( h n ) ) T (M_{k}(h_{1}),...,M_{k}(h_{n}))^{T} (Mk(h1),...,Mk(hn))T这是细胞从 1 1 1到 n n n并分配为组 k k k的概率。组件1损失为 L c 1 ( K 2 ) − 1 ∑ k 1 K − 1 ∑ l k K α ~ k T S α ~ l α ~ k T S α ~ k α ~ l T S α ~ l L_{c1}\begin{pmatrix} K \\ 2 \end{pmatrix}^{-1}\sum_{k1}^{K-1}\sum_{lk}^{K}\frac{\widetilde{\alpha}_{k}^{T}S\widetilde{\alpha}_{l}}{\sqrt{\widetilde{\alpha}_{k}^{T}S\widetilde{\alpha}_{k}\widetilde{\alpha}_{l}^{T}S\widetilde{\alpha}_{l}}} Lc1(K2)−1k1∑K−1lk∑Kα kTSα kα lTSα l α kTSα l其中 α ~ k T S α ~ l \widetilde{\alpha}_{k}^{T}S\widetilde{\alpha}_{l} α kTSα l被最小化代表 M k ( h ) M_{k}(h) Mk(h)和 M l ( h ) M_{l}(h) Ml(h)趋近正交。
对于第二部分组件2将不同组的软分配值推送到 R K R^{K} RK中的单纯形的不同角这也增加了组的可分性。令 e k ∈ R K e_{k}\in R^{K} ek∈RK表示其第 k k k个元素为1其它元素为0的向量。因此 e k e_k ek是单纯形的第 k k k个角。由于 α i \alpha_{i} αi是组识别模块中的第二层输出令 m k ∈ R n m_k\in R^{n} mk∈Rn其第 i i i个元素为 e x p ( − ∣ ∣ α i − e k ∣ ∣ 2 2 ) exp(-||\alpha_{i}-e_{k}||_{2}^{2}) exp(−∣∣αi−ek∣∣22)其测量了软分配 α i \alpha_i αi和 e k e_k ek之间的距离。组件2损失为 L c 2 ( K 2 ) − 1 ∑ k 1 K − 1 ∑ l k K m k T S m l m k T S m k m l T S m l L_{c2}\begin{pmatrix} K \\ 2 \end{pmatrix}^{-1}\sum_{k1}^{K-1}\sum_{lk}^{K}\frac{m_{k}^{T}Sm_{l}}{\sqrt{m_{k}^{T}Sm_{k}m_{l}^{T}Sm_{l}}} Lc2(K2)−1k1∑K−1lk∑KmkTSmkmlTSml mkTSml该损失加强了 e x p ( − ∣ ∣ α i − e k ∣ ∣ 2 2 ) exp(-||\alpha_{i}-e_{k}||_{2}^{2}) exp(−∣∣αi−ek∣∣22)和 e x p ( − ∣ ∣ α i − e l ∣ ∣ 2 2 ) exp(-||\alpha_{i}-e_{l}||_{2}^{2}) exp(−∣∣αi−el∣∣22)的正交因此软分配输出 ( M 1 ( h ) , . . . , M K ( h ) ) T (M_{1}(h),...,M_{K}(h))^{T} (M1(h),...,MK(h))T会倾向于一个单纯形角而不是同时接近多个单纯形角。因此同一组中的低维表示将是紧凑的而来自不同组的低维表示将被分离。
对于第三部分组件3旨在避免大多数细胞被分配给所有组中的一小部分。令 α ‾ ∈ R K \overline{\alpha}\in R^{K} α∈RK表示 α i \alpha_{i} αi i 1 , . . . , n i1,...,n i1,...,n求平均其第 k k k个元素为 α ‾ k \overline{\alpha}_{k} αk。组件3损失为 L c 3 ∑ k 1 K α ‾ k l o g α ‾ k L_{c3}\sum_{k1}^{K}\overline{\alpha}_{k}log \overline{\alpha}_{k} Lc3k1∑Kαklogαk其目的是以相等的概率分配每个组索引比如 α ‾ 1 , . . . , α ‾ K 1 / K \overline{\alpha}_{1},...,\overline{\alpha}_{K}1/K α1,...,αK1/K。
预测损失
模态内预测损失定义为 L W p r e d i c t 1 n V ∑ i 1 n ∑ v 1 V ∣ ∣ x ~ i ( v ) − x i ( v ) ∣ ∣ 2 L_{Wpredict}\frac{1}{nV}\sum_{i1}^{n}\sum_{v1}^{V}||\widetilde{x}_{i}^{(v)}-x_{i}^{(v)}||_2 LWpredictnV1i1∑nv1∑V∣∣x i(v)−xi(v)∣∣2跨模态预测损失为 L C p r e d i c t 1 n ( V 2 ) ∑ i 1 n ∑ v 1 v 2 ∣ ∣ x ~ i ( v 1 , v 2 ) − x i ( v 2 ) ∣ ∣ 2 L_{Cpredict}\frac{1}{n\begin{pmatrix} V \\ 2 \end{pmatrix}}\sum_{i1}^{n}\sum_{v_1v_2}||\widetilde{x}_{i}^{(v_1,v_2)}-x_{i}^{(v_2)}||_2 LCpredictn(V2)1i1∑nv1v2∑∣∣x i(v1,v2)−xi(v2)∣∣2其测量模态 v 1 v_1 v1到 v 2 v_2 v2的跨模态预测与模态 v 2 v_2 v2的原始值之间的距离。
生成损失和判别损失
生成器和判别器通过最小二乘损失训练将模态内预测特征 x ~ i ( v ) \widetilde{x}_{i}^{(v)} x i(v)分配给标签1将原始特征 x i ( v ) x_{i}^{(v)} xi(v)分配给标签0。生成损失为 L G e n 1 n V ∑ i 1 n ∑ v 1 V ∣ ∣ D i s ( v ) ( x ~ i ( v ) ) − 1 ∣ ∣ 2 2 L_{Gen}\frac{1}{nV}\sum_{i1}^{n}\sum_{v1}^{V}||Dis^{(v)}(\widetilde{x}_{i}^{(v)})-1||_{2}^{2} LGennV1i1∑nv1∑V∣∣Dis(v)(x i(v))−1∣∣22判别损失为 L D i s 1 n V ∑ i 1 n ∑ v 1 V ∣ ∣ D i s ( v ) ( x ~ i ( v ) ) ∣ ∣ 2 2 1 n V ∑ i 1 n ∑ v 1 V ∣ ∣ D i s ( v ) ( x i ( v ) ) − 1 ∣ ∣ 2 2 L_{Dis}\frac{1}{nV}\sum_{i1}^{n}\sum_{v1}^{V}||Dis^{(v)}(\widetilde{x}_{i}^{(v)})||_{2}^{2}\frac{1}{nV}\sum_{i1}^{n}\sum_{v1}^{V}||Dis^{(v)}(x_{i}^{(v)})-1||_{2}^{2} LDisnV1i1∑nv1∑V∣∣Dis(v)(x i(v))∣∣22nV1i1∑nv1∑V∣∣Dis(v)(xi(v))−1∣∣22其目的是使鉴别器将模态内预测特征分类为0将原始特征分类为1。这组最小二乘损失提高了训练生成器的质量因为它符合生成器的基本目标即生成具有与原始特征相似分布的特征数据。
对比损失
应用对比损失来对齐来自不同模态的潜在编码定义余弦相似度为 s i , j ( v 1 , v 2 ) ( z i ( v 1 ) ) T z j ( v 2 ) ∣ ∣ z i ( v 1 ) ∣ ∣ 2 ⋅ ∣ ∣ z j ( v 2 ) ∣ ∣ 2 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} si,j(v1,v2)∣∣zi(v1)∣∣2⋅∣∣zj(v2)∣∣2(zi(v1))Tzj(v2)令 l i ( v 1 , v 2 ) − l o g e x p ( s i , j ( v 1 , v 2 ) τ ) ∑ s ′ ∈ N e g ( z i ( v 1 ) , z i ( v 2 ) ) e x p ( s ′ τ ) l_{i}^{(v_1,v_2)}-log\frac{exp(\frac{s_{i,j}^{(v_1,v_2)}}{\tau})}{\sum_{s\in Neg(z_{i}^{(v_1)},z_{i}^{(v_2)})}exp(\frac{s}{\tau})} li(v1,v2)−log∑s′∈Neg(zi(v1),zi(v2))exp(τs′)exp(τsi,j(v1,v2))其中 N e g ( z i ( v 1 ) , z i ( v 2 ) ) Neg(z_{i}^{(v_1)},z_{i}^{(v_2)}) Neg(zi(v1),zi(v2))通过从集合 N i { s i j ( v 1 , v 2 ) : j 1 , . . . , n , j ≠ i , v 1 , v 2 1 , . . . , V , a r g m a x α i ≠ a r g m a x α j } N_{i}\left\{s_{ij}^{(v_1,v_2)}:j1,...,n,j\neq i,v_1,v_21,...,V,argmax\alpha_i\neq argmax\alpha_j\right\} Ni{sij(v1,v2):j1,...,n,ji,v1,v21,...,V,argmaxαiargmaxαj}采样固定数量的元素获得 τ \tau τ为超参数。对比损失为 L C o n δ ⋅ 1 n ( V 2 ) ∑ i 1 n ∑ v 1 v 2 l i ( v 1 , v 2 ) L_{Con}\delta\cdot\frac{1}{n\begin{pmatrix} V \\ 2 \end{pmatrix}}\sum_{i1}^{n}\sum_{v_1v_2}l_{i}^{(v_1,v_2)} LConδ⋅n(V2)1i1∑nv1v2∑li(v1,v2)其中 δ \delta δ为超参数。
分类损失
定义向量 b i b_{i} bi其第 k k k个元素为1其余为0其中 k k k是细胞 i i i的观测类别标注。令 g i n k ( i ) / n g_in_{k(i)}/n gink(i)/n其中 n k n_k nk为类别为 k k k的细胞数量。通过交叉熵评估分类准确率 L e n t r o p y − 1 n ∑ i 1 n g i ⋅ ( l o g ( M 1 ( h i ) ) , . . . , l o g ( M K ( h i ) ) ) T b i L_{entropy}-\frac{1}{n}\sum_{i1}^{n}g_{i}\cdot (log(M_1(h_i)),...,log(M_K(h_i)))^{T}b_i Lentropy−n1i1∑ngi⋅(log(M1(hi)),...,log(MK(hi)))Tbi
训练步骤
作者首先提出了在没有细胞标签的情况下训练UnitedNet用于跨模态预测和无监督组识别的过程。它由两个步骤迭代训练组识别更新步骤和预测更新步骤。在组标识更新步骤中编码器输出模态特定编码并将其送到组识别模块。然后组识别模块将模态特定编码融合为共享潜在代码并获得 K K K-dimensional 软聚类分配。解码器基于模态特定编码输出模态内预测特征 x ~ i ( v ) \widetilde{x}_{i}^{(v)} x i(v)。接下来编码器和聚类模块通过聚类损失进行更新 L g r o u p L c 1 L c 2 L c 3 L_{group}L_{c1}L_{c2}L_{c3} LgroupLc1Lc2Lc3在预测更新步骤中编码器从不同模态输出低维表示将它们送到解码器并获得跨模态预测特征 x ~ i ( v 1 , v 2 ) \widetilde{x}_{i}^{(v_1,v_2)} x i(v1,v2)和模态内预测特征 x ~ i ( v ) \widetilde{x}_{i}^{(v)} x i(v)。然后将预测的特征输入到鉴别器中。鉴别器由鉴别器损失 L D i s L_{Dis} LDis更新。接下来通过模态内预测损失、跨模态预测损失、生成器损失和对比损失之和来更新编码器和解码器 L P G C L W p r e d i c t L C p r e d i c t L G e n L C o n L_{PGC}L_{Wpredict}L_{Cpredict}L_{Gen}L_{Con} LPGCLWpredictLCpredictLGenLCon以上两个训练步骤分别总结在算法1和2中。
算法1 算法2
对于有监督学习为了训练UnitedNet进行监督分类通过以下方式修改组识别损失 L g r o u p L e n t r o p y L_{group}L_{entropy} LgroupLentropy
使用SHAP进行特征相关性分析
为了深入了解不同特征的重要性UnitedNet应用了常用于解释机器学习模型的SHAPSHapley Additive exPlanations。这种方法的思想是在固定其他特征的同时通过线性函数近似特征对输出的影响并且该函数的系数对应于Shapley值。SHAP的优点包括
基于理论的可解释性广泛的应用范围不需要扰动模型或数据的计算过程这是许多其他特征重要性方法所要求的。
下面将首先介绍Shapley值的计算过程然后解释它如何应用于UnitedNet。
假设想评估特征 x j x_j xj对函数 f ( x ) f(x) f(x) x ( x 1 , . . . , x Q ) T x(x_1,...,x_Q)^{T} x(x1,...,xQ)T j ∈ { 1 , . . . , Q } j\in\left\{1,...,Q\right\} j∈{1,...,Q}的重要性。令 F F F为 x x x的特征集 S S S为 F F F的子集。 ∣ F ∣ |F| ∣F∣和 ∣ S ∣ |S| ∣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 F F中删除特征 j j j。 f S ( x S ) f_{S}(x_S) fS(xS)通过对 f ( x S , x F − S ( i ) ) f(x_{S},x^{(i)}_{F-S}) f(xS,xF−S(i))求样本均值得到 x F − S ( i ) x^{(i)}_{F-S} xF−S(i)为不在 S S S中的特征的第 i i i个观测值。Shapley值通过计算去除第 j j j个特征后 f ( x ) f(x) f(x)变化的加权平均值来衡量该特征的重要性。
对于神经网络 f ( x ) f ( 1 ) ⊙ f ( 2 ) ⋅ ⋅ ⋅ f ( L ) f(x)f^{(1)}\odot f^{(2)}\cdot\cdot\cdot f^{(L)} f(x)f(1)⊙f(2)⋅⋅⋅f(L)令第 l l l层的输出维度为 L ( l ) L^{(l)} L(l)第 l l l层的第 q q q个输入特征为 e q ( l ) e_{q}^{(l)} eq(l)其样本均值为 e ‾ q ( l ) \overline{e}_{q}^{(l)} eq(l)。可以用以下递归公式通过Deep SHAP以高效的方式估计模型的Shapley值 ϕ q , r ( f ( l − 1 ) ⊙ f ( l ) , e ( l − 1 ) ) ( e q ( l − 1 ) − e ‾ q ( l − 1 ) ) ⋅ ∑ r ( l − 1 ) 1 L ( l − 1 ) m ( l − 1 ) ( r ( l − 1 ) , q ) ⋅ m ( l ) ( r , r ( l − 1 ) ) \phi_{q,r}(f^{(l-1)}\odot f^{(l)},e^{(l-1)})(e^{(l-1)}_{q}-\overline{e}^{(l-1)}_{q})\cdot \sum_{r^{(l-1)}1}^{L^{(l-1)}}m^{(l-1)}(r^{(l-1)},q)\cdot m^{(l)}(r,r^{(l-1)}) ϕq,r(f(l−1)⊙f(l),e(l−1))(eq(l−1)−eq(l−1))⋅r(l−1)1∑L(l−1)m(l−1)(r(l−1),q)⋅m(l)(r,r(l−1))其中 ϕ q , r ( f ( l − 1 ) ⊙ f ( l ) , e ( l − 1 ) ) \phi_{q,r}(f^{(l-1)}\odot f^{(l)},e^{(l-1)}) ϕq,r(f(l−1)⊙f(l),e(l−1))衡量了 e q ( l − 1 ) e_{q}^{(l-1)} eq(l−1)相对于 f ( l − 1 ) ⊙ f ( l ) f^{(l-1)}\odot f^{(l)} f(l−1)⊙f(l)的输出的第 r r r个元素的重要性。 m ( l ) ( r , r ( l − 1 ) ) ϕ r ( l − 1 ) , r ( f ( l ) , e ( l ) ) / ( e r ( l − 1 ) ( l ) − e ‾ r ( l − 1 ) ( l ) ) m^{(l)}(r,r^{(l-1)})\phi_{r^{(l-1)},r}(f^{(l)},e^{(l)})/(e^{(l)}_{r^{(l-1)}}-\overline{e}^{(l)}_{r^{(l-1)}}) m(l)(r,r(l−1))ϕr(l−1),r(f(l),e(l))/(er(l−1)(l)−er(l−1)(l)) m ( l − 1 ) ( r ( l − 1 ) , q ) ϕ q , r ( l − 1 ) ( f ( l − 1 ) , e ( l − 1 ) ) / ( e q ( l − 1 ) − e ‾ q ( l − 1 ) ) m^{(l-1)}(r^{(l-1)},q)\phi_{q,r^{(l-1)}}(f^{(l-1)},e^{(l-1)})/(e^{(l-1)}_{q}-\overline{e}^{(l-1)}_{q}) m(l−1)(r(l−1),q)ϕq,r(l−1)(f(l−1),e(l−1))/(eq(l−1)−eq(l−1))为了识别与特定组具有高度相关性的特征对于每个细胞需要计算每个输入特征相对于该组的软分配的Shapley值。然后将来自被分类为该组的细胞的所有Sharpley值作为绝对值并用于计算平均值。具有最高平均值的前 n n n个特征被解释为与该组具有高度相关性Patch-seq GABAergic数据集的 n 7 n7 n7multiome ATACGene expression BMMCs数据集的 n 20 n20 n20。
为了量化组内跨模态特征与特征的相关性作者考虑了在前一步中选择的该组的高相关性特征。接下来计算每个特征相对于来自另一模态的每个其他特征的Shapley值。然后通过平均来聚集来自不同细胞的Sharpley值的绝对值。值相对较大的特征被认为是重要的。
评价指标
无监督和有监督组识别用ARIadjusted rand index预测性能用 R 2 R^2 R2coefficient of determination和AUCarea under the ROC curve两个任务之间的关系用Pearson’s correlation评估。
对于ARI 将从模型中获得的簇与从细胞类型标签中获得的簇进行比较。令 a k 1 a_{k_1} ak1 k 1 1 , . . . , K k_11,...,K k11,...,K表示来自模型的第 k 1 k_1 k1个簇中的细胞数 b k 2 b_{k_2} bk2 k 1 1 , . . . , K k_11,...,K k11,...,K表示来自细胞类型标签的第 k 2 k_2 k2个簇的细胞数。 n k 1 , k 2 n_{k_1,k_2} nk1,k2表示来自模型的第 k 1 k_1 k1个簇和来自细胞类型标签的第 k 2 k_2 k2个簇中的观测次数。当来自模型的聚类结果接近于来自观察到的细胞类型标签的聚类结果时ARI接近于1并且对于随机猜测接近于0。
对于 R 2 R^2 R2 令 y y y和 y ~ \widetilde{y} y 表示观测数据和预测数据 y ‾ \overline{y} y被定义为与 y y y长度相同的向量并且每个元素都是 y y y的样本均值。 R 2 R^{2} R2为 R 2 1 − ∣ ∣ y ~ − y ∣ ∣ 2 2 ∣ ∣ y ‾ − y ∣ ∣ 2 2 R^{2}1-\frac{||\widetilde{y}-y||_{2}^{2}}{||\overline{y}-y||_{2}^{2}} R21−∣∣y−y∣∣22∣∣y −y∣∣22它将来自模型的预测 y ~ \widetilde{y} y 的均方误差MSE与以常数值 y ‾ \overline{y} y作为预测的基线的MSE进行比较。值在负无穷到1之间当预测等于观测时 R 2 1 R^21 R21。
对于Pearson’s correlation 令 y ^ \widehat{y} y 表示与 y ~ \widetilde{y} y 等长度的向量其每个元素为 y ~ \widetilde{y} y 的样本均值。因此 y y y与 y ~ \widetilde{y} y 之间的Pearson’s correlation为 r y ~ − y ^ , y − y ‾ ∣ ∣ y ~ − y ^ ∣ ∣ 2 ⋅ ∣ ∣ y − y ‾ ∣ ∣ 2 r\frac{\widetilde{y}-\widehat{y},y-\overline{y}}{||\widetilde{y}-\widehat{y}||_{2}\cdot||y-\overline{y}||_{2}} r∣∣y −y ∣∣2⋅∣∣y−y∣∣2y −y ,y−y当预测与观测数据分别具有正或负线性关系时它取[-1,1]中的值比如等于1或负1。当预测和观测数据没有线性关系时它等于零。
对于AUC 为了对multiome ATACGene expression数据集进行建模作者对DNA可及性数据进行了二值化。因此为了评估这些数据的预测我们采用了ROC曲线下的面积。设 n 0 n_0 n0和 n 1 n_1 n1分别表示来自观测数据的0和1的数目。 p 0 , i p_{0,i} p0,i i 1 , . . , n 0 i1,..,n_0 i1,..,n0和 p 1 , j p_{1,j} p1,j j 1 , . . , n 1 j1,..,n_1 j1,..,n1分别表示两组观测值的模型预测。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 ) (s_{i}^{1},s_{i}^{2}) (si1,si2)和模态 v v v的第 i i i行 x i ( v ) x_{i}^{(v)} xi(v)对应于细胞或spots i i i 的情况下可以计算模态 v v v的niche模态用 x ( v n i c h e ) x^{(v\thinspace niche)} x(vniche) v 1 , . . , V v1,..,V v1,..,V。对于细胞 i i i x i ( v n i c h e ) x^{(v\thinspace niche)}_{i} xi(vniche)为 x i ( v n i c h e ) ∑ j 1 J x j ( v ) ⋅ w i j x^{(v\thinspace niche)}_{i}\sum_{j1}^{J}x^{(v)}_{j}\cdot w_{ij} xi(vniche)j1∑Jxj(v)⋅wij其中 j ∈ { 1 , . . . , J } j\in\left\{1,...,J\right\} j∈{1,...,J}表示属于cell或spot i i i 的 J J J-nearest的cells或spots w i j w_{ij} wij为 w i j 1 / d i s t a n c e { ( s i 1 , s i 2 ) , ( s j 1 , s j 2 ) } ∑ j 1 J 1 / d i s t a n c e { ( s i 1 , s i 2 ) , ( s j 1 , s j 2 ) } w_{ij}\frac{1/distance\left\{(s_{i}^{1},s_{i}^{2}),(s_{j}^{1},s_{j}^{2})\right\}}{\sum_{j1}^{J}1/distance\left\{(s_{i}^{1},s_{i}^{2}),(s_{j}^{1},s_{j}^{2})\right\}} wij∑j1J1/distance{(si1,si2),(sj1,sj2)}1/distance{(si1,si2),(sj1,sj2)}其中 d i s t a n c e distance distance表示两个向量之间的欧几里得距离。
DBiT-seq embryo数据集 作者使用DBiT-seq embryo胚胎数据集其中采用了DBiT-seq 936个spots的以下三种模态mRNA表达、蛋白质表达和niche mRNA表达。
对于mRNA表达模态使用scanpy的函数scanpy.pp.normalize_total对原始计数矩阵进行归一化并选择前568个差异表达基因。对于蛋白质表达对原始计数矩阵进行了归一化并使用了22种蛋白质。niche模态是基于标准化的mRNA表达产生的。
对于组织区域表征的第一项任务从原始研究中提取真实组织区域标签这是基于HE图像的主要组织区域的解剖学注释。作者将UnitedNet的聚类结果与其他最先进方法的聚类结果进行了比较。通过ARI来验证它们的性能。对于跨模态预测的并行任务尽管有三种模态被用作UnitedNet模型的输入但作者专注于第一种和第二种模态之间的预测mRNA表达和蛋白质表达。由于DBiT-seq公共数据集中只有一个批次作者将DBiT-seq embryo数据集中的936个spots分为用于预测任务的训练数据集80%748个spots和测试数据集20%188个spots。
DLPFC数据集 作者使用了12个批次的成人背外侧前额叶皮层DLPFC数据集。使用以下三种模式mRNA表达、从HE染色图像中提取的形态学特征和 niche mRNA。对于mRNA表达的模态对原始计数矩阵进行归一化并选择前2365个差异表达基因。使用预先训练的卷积神经网络从stLearn实现的HE染色图像中提取形态学特征。50维形态特征被用作每个spots的第二模态。对于有监督组识别任务使用11个批次及其组织区域注释来训练UnitedNet模型。然后将训练后的模型应用于剩余批次以识别组织区域注释并在HE图像特征和mRNA表达之间进行跨模态预测。作者将UnitedNet的识别性能与来自原始DLPFC论文的SpatialDE PCA和 pseudobulk PCA进行了比较。在识别任务之后作者对最近邻居中具有35个spots的SpaGCN之后的聚类结果应用了细化步骤。
结果
特征与特征相关性
作者使用事后可解释学习SHAP对训练的UnitedNet进行剖析以表明Patch-seq GABAergic数据集中的特征相关性。具体而言使用SHAP将重要性值称为Shapley值分配给关于任何给定模型输出的每个输入特征例如特定识别的细胞组或特定特征的跨模态预测。根据定义具有高Shapley值的特征是有影响力的。因此可以根据特征对Shapley值的排名来选择特征。
以Pvalb神经元类型一种细胞大类为例作者定性地验证了SHAP选择的相关性图2。对于组与特征的相关性SHAP成功地选择了Pvalb神经元差异表达的基因、电生理特征和形态学特征的子集图2ad–f。
对于Pvalb神经元特异性跨模态特征与特征相关性图2bc发现在使用长方电流步长的膜片钳电刺激过程中long square current steps基因Lrrc38与Pvalb神经细胞平均放电率的电生理特征average firing rate表现出更高的相关性。这一结果与之前的研究一致表明Lrrc38相关蛋白是大钾BK通道最关键的调节剂之一大钾通道对神经元放电动力学和神经递质释放至关重要。这些结果表明UnitedNet可以潜在地用于促进Patch-seq数据的细胞类型特异性基因功能相关性的鉴定。
图2特征相关性解释Patch-seq GABAergic数据集。
然后作者探索了UnitedNet中可解释学习对ATAC基因表达数据集的生物学价值。以CD8T主要细胞类型CD8T和CD8T naive细胞为例图3a结果首先鉴定了CD8T主细胞类型差异表达的基因亚群例如CD8A、A2M、LEF1和NELL2和DNA可及性位点例如D8A、DPP8、KDM2B和KDM6B图3b、d、e。在这些基因和DNA可及位点中DNA可及性位点PROS1、KDM2B和KDM6B与CD8T细胞特异性基因表现出更强的相关性表明它们在CD8T功能中的关键作用图3c。该结果与先前的研究一致即CD8T细胞中PROS1表达水平的升高是防止免疫反应过度活跃的关键调节信号。同时KDM2B表达的缺乏启动了T细胞白血病的发生。值得注意的是最近的一项研究发现与KDM2B属于同一基因家族的KDM6B通过诱导效应相关基因中的DNA可及性直接调节CD8T细胞的产生。结果进一步表明KDM2B也可能在调节CD8T细胞的产生中发挥重要作用这是模型发现的这一点是重要的。 图3特征相关性解释ATACGene表达数据集。
空间组学上的应用
空间组学是一种重要的模态技术可以测量完整组织中的空间分辨多组学信息。然而在分析用于组识别任务的空间组学数据时空间信息往往没有被充分利用。联合网络可以灵活地整合不同的模态作为输入包括空间信息。UnitedNet可以利用细胞niche信息每个细胞的邻域基因表达信息作为额外的模态来识别具有生物学意义的群体并增强跨模态预测图4a。
作者首先将UnitedNet应用于单批次的DBiT-seq embryo数据集该数据集同时映射了胚胎组织上的整个转录组和22种蛋白质。具体而言作者生成了编码空间信息的细胞niche信息RNA表达的加权平均值作为第三种分析模态。然后联合网络将基因表达、蛋白质和niche模态相结合用于组织区域的无监督联合识别以及基因表达和蛋白质之间的跨模态预测。通过将原始报告中的组织区域的解剖注释视为Ground Truth对组织区域识别的准确性进行了基准测试。与最先进的方法相比UnitedNet实现了更高的无监督组识别精度图4b。此外联合网络使几个代表性基因和蛋白质表达之间的空间分辨跨模态预测成为可能图第4c。 图4无监督组识别joint注释转移空间组学跨模态预测。联合网络的空间组学数据分析管道示意图。空间组学同时测量完整组织网络中的空间分辨多组学数据。UnitedNet提取细胞邻域信息作为附加模态以及用于无监督或有监督群组识别和跨模态预测的其他模态。
接下来作者将UnitedNet应用于一个带注释的多批次空间组学数据集用于同时有监督的联合组识别和跨模态预测。作者使用了人类背外侧前额叶皮层DLPFC数据集该数据集在空间上绘制了12批次DLPFC脑切片的基因表达和HE染色。同样使用基因表达、基于HE染色的形态学特征和细胞niche模态作为UnitedNet的输入。UnitedNet可以成功地注释看不见的DLPFC切片并实现比其他基准测试方法和未使用交替训练方案或细胞小众模式的消融版UnitedNet更高或相当的精度图4d。
此外作者探讨了UnitedNet是否可以减少空间DLPFC数据集中的批次效应其方式与ATAC基因表达BMMC数据集的分析类似。结果表明与其他消融研究相比UnitedNet在潜在空间中保持了良好的可分性和减少批次效应的能力使其在组识别任务中具有更高或可比的性能。此外UnitedNet使几个代表性基因和HE形态特征之间的跨模态预测成为可能图4e。
