4+单细胞hdWGCNA+机器学习+空转+实验，卷王必备技能，学会直接毕业！！

‍

导语

今天给同学们分享一篇生信文章“Deciphering the Regulatory Networks of the Migrasome‐Associated Cell Subpopulation in Heterotopic Ossification via Multi‐Omics Analysis ”，这篇文章发表在FASEB J 期刊上，影响因子为4.2。

结果：

3.1. 单细胞数据的质控与标准化

本研究使用的单细胞数据来源于GSE126060 数据集。根据特定参数进行数据过滤，并去除双细胞，最终保留 23,922 个细胞（图 S1A,B ）。作者鉴定了 2000 个高变基因（图 S1C ），随后进行数据标准化、均匀性调整、PCA 和 Harmony 分析（图 S1D–F ）。

3.2. 细胞注释

使用 UMAP 进行降维后，识别出 17 个亚群（图 1A ）。随后对每个亚群进行注释，并将所有群组分配到 11 种细胞类型：巨噬细胞、间质谱系细胞（MLin）、内皮细胞、平滑肌细胞（SMC）/周细胞（Peri）、红细胞、T 细胞、卫星细胞、施万细胞、多形核白细胞（PMNs）、浆细胞样树突状细胞（pDC）和 B 细胞（图 1B ）。使用气泡图（图 1C ）可视化这 11 种细胞类型的经典标记物表达水平，并在条形图（图 1D ）中展示每种细胞类型的比例。随后，使用 AUCell 对单个细胞的代谢通路基因得分进行量化。使用 pheatmap 包生成热图，以可视化不同细胞类型代谢通路的活动矩阵（图 1E ）。结果表明，巨噬细胞、MLin、内皮细胞和 SMC/Peri 在大多数通路中表现出相对较高的活性。

细胞注释。(A) 根据 PCA 中识别出的显著成分，利用 UMAP 算法将细胞分为 17 个簇。(B) 17 个细胞簇的注释。(C) 11 种细胞类型和细胞标志物的 Doplot 气泡图。(D) 两组样本中 11 种细胞类型百分比的差异。(E) 展示不同细胞类型中不同代谢通路活性的热图。

3.3. 迁移体评分定量、细胞通讯分析和二次聚类

接下来，作者量化了 migrasome 评分，结果揭示 MLin 亚群不仅在所有细胞类型中表现出最高的 migrasome 评分，而且与对照组和疾病组之间的 migrasome 评分差异最大，在疾病条件下有明显的上调（图 2A ）。这些发现表明 MLin 与疾病特异性 migrasome 活性之间可能存在潜在关联。因此，作者选择 MLin 进行进一步分析。根据|avg_log2FC| > 0 和调整后 p 值（p_val_adj）< 0.05 的标准，在该细胞类型中鉴定出 1384 个差异表达基因（DEGs），如图形火山图（图 2B ）所示。这些 DEGs 被用作后续分析的基因集。然后，作者根据中位数 migrasome 评分将 MLin 分为 High_MLin 和 Low_MLin 两组。使用 CellChat 软件包分析细胞间通讯。在多种细胞群体之间观察到广泛的配体-受体对网络，其特征是高连通性和各种相互作用强度，表明复杂的细胞间相互作用（图 2C ）。值得注意的是，High_MLin 与其他细胞类型表现出特别活跃的信号行为（图 2D ）。通过 PCA、Harmony、ElbowPlot 和 FindClusters 对 MLin 进行进一步亚群分析，识别出 6 个亚型（图 2E–G ）。亚型 C1 是仅在疾病样本中发现的特定亚群，被确认为疾病特异性亚群（图 2H ）。简而言之，这些发现表明，通过迁移体评分识别的 MLin 亚群可能在异位骨化的进展中发挥关键作用。

迁移体评分的量化、细胞通讯分析及二次聚类。(A) 迁移体评分的量化及两组间 11 种细胞类型的差异。(B) MLin 中差异基因的火山图；蓝色代表差异下调基因，红色代表差异上调基因。(C) 12 种阳性细胞类型间的通讯相互作用网络；边框宽度代表细胞间通讯的概率和强度。(D) 12 种阳性细胞类型间通讯网络中相互作用计数的比较，从左到右依次减少，最高为 High_MLin。(E) 每个主成分的方差排序图。(F) PCA 主成分的可视化及分布；颜色代表样本，点代表细胞。(G) 根据 PCA 中可用的主成分，通过 UMAP 算法将细胞分为 6 个簇。(H) 两组样本间 6 个簇比例的差异。

3.4. 基于高维 WGCNA 和机器学习的 HO 关键基因识别

为识别 MLin 亚型内的基因共表达网络，作者进行了 hdWGCNA 分析。在此分析中，设置 group.by 参数为聚类亚型，并将自定义亚型用于构建共表达网络，旨在探索参与异位骨化发展的生物标志物。软阈值通过“TestSoftPowers”函数确定，并设置为 6（图 3A ）。该分析检测到 5 个基因模块：青色、蓝色、棕色、黄色和绿色模块（图 3B ）。模块与模块之间以及模块与模块效应基因（MEs）的进一步分析（图 3C,D ）显示，青色模块的 MEs 水平在疾病亚组中升高（图 3E ），因此选择了青色模块中特征基因连接度（kMEs）最高的前 100 个基因作为进一步分析的候选特征基因。随后，作者使用 GO（图 S2A ）和 KEGG（图 S2B ）富集分析对 hdWGCNA 获得的基因进行功能预测。接下来，作者将 MLin 中的差异表达基因（DEGs）与 hdWGCNA 的基因进行交集分析，得到 100 个交集基因（图 3F ）。在这些基因中，仅保留在批量 RNA 测序数据集 GSE126118 中可检测表达的基因作为 RF 算法的输入特征。根据特征重要性排序，确定了 Ppib、Tagln、Myl9、Dad1 和 Fkbp11 这 5 个基因作为特征基因（图 3G ）。作者使用临床样本数据集 GSE106253 进行外部验证，并对这 5 个基因进行了差异分析。结果显示疾病组中 Ppib（肽酰脯氨酰异构酶 B）和 Tagln（转凝胶蛋白）显著下调（图 3H ），这两个分子被选为 HO 关键基因进行进一步分析。

通过 hdWGCNA 和 ML 筛选出 HO 的关键基因。（A）软阈值选择。（B）5 个基因模块的树状图。（C）不同模块中特征基因的 kMEs。（D）带有 ME 染色的 MLin 的 UMAP 图。（E）5 个 MLin 簇的模块活性。（F）显示 hdWGCNA 差异基因和模块基因交集的维恩图。（G）RF 选定的特征基因。（H）两组样本间基因差异分析。

3.5. 免疫浸润分析

微环境由免疫细胞、细胞外基质（ECM）、生长因子和炎症因子组成，对疾病诊断、生存结果和临床治疗敏感性有显著影响。作者进行了免疫浸润分析，展示了免疫细胞浸润水平和免疫细胞类型之间的相关性（图 S3A,B ）。分析显示，M1 巨噬细胞、肥大细胞、CD8 记忆 T 细胞和 CD8 未成熟 T 细胞在疾病组和对照组之间存在显著差异（图 S3C ）。作者进一步研究了关键基因与免疫细胞之间的关系，发现 Ppib 与 M1 巨噬细胞、M2 巨噬细胞和 CD8 未成熟 T 细胞显著正相关，但与肥大细胞和 CD8 记忆 T 细胞显著负相关（图 S3D ）。

3.6. 关键基因的 miRNA 网络及关键基因与疾病调控基因在 HO 中的相关性分析

作者使用 miRDB 数据库对关键基因进行反向预测，并鉴定出 56 个 miRNA，形成了 56 个 mRNA−miRNA 对，这些对通过 Cytoscape（图 4A ）可视化为非编码 RNA 网络。作者从 GeneCards 中检索了与 HO 相关的疾病调控基因。作者分析了在 GSE126118 批量 RNA 测序数据中也表达的、按相关性得分排名前 20 的基因的表达水平，发现 Sox9、Spi1、Mmp9、Bmp2、Smad6、Adrb2 和 Cdh11 在两组间表现出显著的表达差异（图 4B ）。作者还对关键基因与疾病调控基因进行了相关性分析。在气泡图中显示的结果表明，Tagln 与 Adrb2 之间存在显著正相关（r = 0.968），而 Ppib 与 Smad6 之间存在显著负相关（r = −0.977）（图 4C ）。作者从 GeneCards 中获取了高相关性得分的 HO 疾病调控基因，并分析了它们与关键基因的相关性。共表达网络揭示了关键基因和疾病相关基因在单细胞中的表达及其相关性。值得注意的是，疾病基因 Spi1、Bmp2、Mmp9、Smad6 和 Sox9 与 Tagln 表现出显著相关性。其中，Bmp2 和 Smad6 与 Tagln 呈正相关，而 Spi1、Mmp9 和 Sox9 与 Tagln 呈负相关（图 5A–E ）。类似地，疾病基因 Spi1、Bmp2、Mmp9、Smad6 和 Sox9 均与 Ppib 呈负相关（图 S4A–E ）。

关键基因的 miRNA 网络以及 HO 中关键基因与疾病调控基因之间的相关性。(A) 关键基因的 miRNA 网络，红色表示 mRNA，蓝色表示 miRNA。(B) 疾病调控基因的表达差异。(C) 关键基因与疾病调控基因的 Pearson 相关性分析，蓝色表示负相关，红色表示正相关。

分析疾病调控基因与 Tagln 在单细胞中的共表达和相关性。(A–E) 分析每个疾病调控基因与 Tagln 的共表达和相关性，包括 UMAP 可视化、基因表达相关性热图和散点回归分析。

3.7. SCENIC 分析

考虑到 MLin 在疾病组和对照组之间的迁移体评分差异最大，作者选择了 MLin 亚组进行 SCENIC 分析，并输出了该亚组内的所有调控子（Step2_regulonTargetsInfo.tsv）。根据中位数迁移体评分，将 MLin 分为高评分（Hsco）和低评分（Lsco）两组。生成了热图以可视化每个细胞的调控子活性评分（图 6A ）。作者计算了每个调控子的调控子特异性评分（RSS），以评估其在 Hsco 和 Lsco 组之间的特异性。结果使用气泡图进行可视化（图 6B ）。随后，根据关键基因 Ppib 和 Tagln 的表达水平，将 MLin 分为高表达（Hexp）和低表达（Lexp）两组。使用散点图可视化了 Hexp 组中前 5 个特异性最高的调控子。在高 Ppib 表达的 MLin 中，Jund、Creb3、Atf4、Jun 和 Junb 表现出最高的调控子特异性（图 6C ）。类似地，在高 Tagln 表达的 MLin 中，Creb3、Sox4、Jund、Atf4 和 Jun 显示出最大的特异性（图 6D ）。

SCENIC 分析。(A)热图显示每个细胞的调控子活性评分。(B)气泡图显示调控子在 Hsco 和 Lsco 中的特异性。(C)散点图显示 Ppib 高表达组(Hexp)中的前 5 个调控子。(D)散点图显示 Tagln 高表达组(Hexp)中的前 5 个调控子。

3.8. 细胞发育轨迹和免疫代谢通路分析

通过首先计算细胞间的相似性并构建细胞分化轨迹，进行了伪时间分析。通过可视化这条轨迹，作者可以追踪细胞的发育过程，并研究不同时间点的基因表达模式。伪时间是基于细胞基因表达信息由 Monocle 计算出的概率，反映了时间顺序。使用伪时间值和细胞类型对细胞进行着色（图 7A,B ）。在分支的不同位置观察到差异基因表达，这在整体热图中没有体现。作者识别了在分支点前后表达差异较大的基因，并在分支热图中进行可视化（图 7C ）。HO 关键基因 Ppib 和 Tagln 的表达在细胞分化过程中总体呈上升趋势（图 7D ）。作者使用 AUCell 对单个细胞的免疫和代谢通路基因得分进行定量，并生成了气泡图以展示 HO 关键基因在这些通路中的活性差异。结果表明 Ppib 和 Tagln 在上皮间质转化和血管生成等通路中具有高活性（图 7E ）。

细胞发育轨迹。(A, B) 细胞伪时间分析和发育轨迹。(C) 每个色素细胞分支的基因表达动态。(D) 关键基因表达与细胞发育轨迹之间的关系。(E) 显示 Tagln 和 Ppib 在不同通路中富集结果的气泡图。

3.9. 解卷积分析及空间转录组学中关键基因的表达

解卷积结果显示了样本中每个位置的细胞类型比例。在GSM8458769 中，大多数位置被 MLin 和卫星细胞占据。在 GSM8458770 中，大多数位置被 MLin 和 PMNs 占据（图 8A,B ）。为验证解卷积的准确性，作者使用 FindAllMarkers 函数为每种细胞类型鉴定差异表达基因（DEGs）（图 8C ）。使用的过滤参数为：logfc.threshold = 0，min.pct = 0.1，only.pos = F。然后作者分析了空间转录组学中关键基因的表达情况，通过散点图（图 8D ）和气泡图（图 8E ）可视化 Ppib 和 Tagln 的表达。散点图显示了关键基因在空间转录组学数据中的特定分布。

解卷积分析和空间转录组中 HO 关键基因的表达。(A, B) 每个点的细胞类型和比例。(C) 每种细胞类型的差异基因。(D, E) 空间转录组样本中 HO 关键基因的表达。

3.10. 实验验证 HO 中关键基因的表达

为验证关键基因的表达水平，进行了 qRT–PCR 实验。结果表明，在 HO 模型中，Ppib 和 Tagln 基因的表达在受伤后 1 周（图9A ）和 4 周（图 9B ）显著上调，随后在 7 周（图 9C ）和 10 周（图 9D ）显著下调。

HO 中关键基因表达的实验验证。qRT–PCR 分析小鼠 HO 模型在受伤后 1（A）、4（B）、7（C）和 10（D）周的 Ppib 和 Tagln 表达。

总结

本研究通过单细胞转录组学的迁移体评分分析，确定了 MLin 是异位骨化进展过程中的关键细胞亚型。随后，利用 hdWGCNA、RF 算法及临床样本数据集验证，Ppib 和 Tagln 被鉴定为参与异位骨化的关键分子。整合 bulk transcriptome 进行免疫浸润、miRNA 网络、共表达网络、转录因子调控网络及信号通路分析，从而阐明异位骨化发生发展背后的潜在调控机制。此外，空间转录组学进一步揭示了细胞亚型分布的空间模式及关键基因表达。最终，实验验证表明 Ppib 和 Tagln 表现出与异位骨化进展相关的动态表达模式。这些发现为未来异位骨化的临床诊断和治疗策略提供了新的启示和有前景的方向。对这篇文章感兴趣的老师，欢迎扫码咨询！

往期推荐

纯生信选刊

非肿瘤生信

预后模型

单基因生信

单细胞系列