4+机器学习+单细胞+实验验证，干湿结合经典不容错过！

‍

导语

今天给同学们分享一篇生信文章“Multi-dataset identification of innovative feature genes and molecular mechanisms in keratoconus”，这篇文章发表在J Cell Mol Med期刊上，影响因子为4.3。

结果：

差异基因的鉴定和途径富集

作者从NCBI GEO公共数据库下载了GSE 204791数据集，其中包括16例病例：8例正常组和8例疾病组。采用Limma软件包计算两组差异表达基因，筛选标准为p < 0.05，|log2光纤通道|>1.总共鉴定了1066个差异表达的基因，包括593个上调和473个下调的基因（图2A，B）。随后，作者用通路分析法分析了差异表达的基因。使用Metascape数据库，差异表达的基因主要富集在途径中，例如次级代谢过程、类维生素A代谢过程和单加氧酶活性（图3A）。这些途径的网络显示在图3B中。

关键基因的筛选及诊断模型的构建

作者分别使用GSE 204791和GSE 77938数据集作为训练集和验证集。使用Lasso回归选择在前一步骤中鉴定的差异表达基因用于特征选择。结果显示，通过Lasso回归鉴定了六个基因作为KC的特征基因，即ATOH 7、DBNDD 1、RNF 217-AS 1、ARL 11、MYRF和SNORA 74 B（图4A-4C）。除了ATOH 7，KC组中其他关键基因的表达水平显著上调。将这6个基因作为核心基因进行后续实验，构建KC诊断的预测模型。模型公式为：风险评分= ATOH 7 ×（− 0.063870942255986）+ DBNDD1 × 0.00690986717583838 + RNF217-AS1 × 0.007515606011362 + ARL11 × 0.00814156838823462 + MYRF × 0.0290731999084615 + SNORA74B × 0.198609773078782。结果显示，由六种基因构建的预测模型具有令人满意的诊断性能，AUC为1（图4D）。GSE77938数据集用作验证集，以从外部验证预测模型。结果显示，该模型具有强稳定性，AUC为0.832，经验证集GSE 77938验证（图4E）。

关键基因的单细胞分析

作者从GSE146276下载了单细胞数据，并使用Seurat软件包进行了单细胞分析。使用tSNE算法对细胞进行聚类，并且使用R包中的SingleR来注释每个聚类。将所有亚型注释为三个细胞类别：上皮细胞、角质形成细胞和成纤维细胞。其中，上皮细胞占最高比例（图5A）。这三种细胞类型中关键基因的表达如图5B、C所示。作者发现ATOH 7、DBNDD 1、ARL 11和MYRF都在上皮细胞中表达，并且表达水平的排序为DBNDD 1、ARL 11、MYRF和ATOH 7（图5B）。ARL11是角质形成细胞中表达最多的基因，并且在成纤维细胞中观察到这些特征基因的表达低于除MYRF之外的其他细胞类型（图5C）。

免疫浸润分析

KC患者的角膜免疫微环境，包括免疫细胞、细胞因子、趋化因子和免疫介质，通常会发生改变。通过分析KC数据集中关键基因与免疫浸润之间的关系，作者进一步探索了关键基因影响KC进展的潜在机制。图6A、B显示了每个患者中免疫细胞的比例和免疫细胞之间的相关性。天然B细胞在所有样品中占最高比例（图6A）。活化的记忆性CD 4+ T细胞和浆细胞之间存在正相关性（Pearsonr = 0.38;图6 B）。此外，结果显示KC组中活化的记忆性CD 4+ T细胞和浆细胞的水平显著高于正常组（p< 0.05;图6C）。六种基因与免疫细胞之间的相关性显示在图6D-I中。浆细胞与ARL 11、DBNDD 1和MYRF具有显著的正相关性（所有p< 0.05，图6D、F、G）。活化的记忆性CD 4+ T细胞与ARL 11和SN 0 RA 74 B具有显著的正相关性（均p< 0.05;图6D，I），但与AT 0 H 7具有负相关性（p< 0.01;图6 E）。

基于TISIDB数据库，作者发现这六个关键基因与不同的免疫因子（包括趋化因子、免疫抑制剂、免疫刺激剂、MHC和细胞受体）之间存在显著相关性（图7A-E）。对于MYRF、RNF 217-AS 1和SNORA 74 B观察到类似的趋势，与趋化因子、免疫抑制剂、受体和MHC相比具有更大的负相关性（图7A、B、D、E），但与免疫刺激剂相比具有更强的正相关性，尤其是与CD 86相比（Pearsonr = 1，p< 0.01;图7 C）。此外，ATOH 7与免疫因子之间的相关性几乎与MYRF相反，与免疫抑制剂（VTCN 1，Pearsonr = 1，p< 0.01;图7 B）、趋化因子和MHC（图7 A、B、D）具有明显的正相关性。这些结果表明，这些关键基因与免疫细胞浸润密切相关，并在免疫微环境中发挥重要作用。

基因集富集分析

接下来，作者研究了富含关键基因的特定信号通路，以探讨潜在的分子机制。选择高度显著的途径用于集中显示。ARL 11富集的途径包括抗坏血酸和乙醛酸代谢、甘氨酸丝氨酸和苏氨酸代谢以及鞘糖脂合成神经节系列（图8A）;富含ATOH 7的途径包括脂肪细胞因子信号传导途径、黑素生成和氮代谢。（图8 B）;富含DBNDD 1的途径包括DNA复制、戊糖和葡萄糖醛酸相互转化以及趋化因子信号传导途径（图8 C）;富含MYRF的途径包括MAPK信号传导途径、戊糖磷酸途径和味觉转导途径。（图8D）;富含RNF217-AS 1的途径包括B细胞受体信号传导途径、黑素生成和嗅觉转导（图8 E）;并且富含SNORA 74 B的途径包括钙粘附蛋白连接、RNA降解和磷酸肌醇代谢（图8 F）。

预言转录因子

作者使用关键基因作为基因集，以进一步探索KC中涉及的转录因子调控网络。使用Cistrome DB在线数据库预测相关转录因子。其中，DBNDD 1预测了99个转录因子，ARL 11推断了81个转录因子，MYRF预测了48个转录因子，ATOH 7预测了76个转录因子（图9）。本研究鉴定了几种与免疫应答相关的转录因子。例如，关于IRF家族，在ATOH7、DBNDD 1、MYRF和ARL 11中发现了IRF 4，并且预测MYRF为IRF 5（图9）。在STAT家族中，在MYRF和DBNDD 1中发现了STAT 4，在ARL 11中发现了STAT 3（图9）。作者使用Cytoscape构建了KC中关键基因的综合转录调控网络以进行可视化（图9）。

关键基因的疾病调控网络

作者从GeneCards数据库（https://www.genecards.org/）中获得了1802个KC相关基因。对疾病基因的表达差异的分析显示，基因如KC 6、LCA 5、SPATA 7和ZNF 469的表达在正常组和疾病组之间不同（图10A）。接下来，作者进行了关键基因和KC相关基因之间的相关性分析。作者发现关键基因的表达水平与KC相关基因的表达水平显著相关（图10 B），MYRF和SPATA 7显示出显著的负相关性（Pearsonr =-0.788，p< 0.01），MYRF和KC 6显示出显著的正相关性（Pearsonr = 0.754，p< 0.01）。ATOH 7与SPATA 7呈显著正相关（p< 0. 05）。01;图10 B）和与KC 6的显著负相关性（p< 0.05;图10 B）。

上游miRNA网络分析与潜在治疗药物预测

作者利用Mircode数据库对关键基因进行了反向预测，获得了49个miRNA和74个mRNA-miRNA关系对。使用Cytoscape可视化结果（图11）。数据显示，miR-150与DBNDD 1、SNORA 74 B和ARL 11具有强相互作用（图11）。作者将差异表达的mRNA分为上调组和下调组，并使用CMap数据库预测差异表达基因的药物靶点。结果显示，受药物干扰的表达谱，如恩替司他（图12 A）、赭曲霉毒素-a（图12 B）、苯环酮（图12 C）和GSK-3-β-葡聚糖-II（图12 D），与疾病干扰的表达谱显著负相关，并且药物可能减轻或甚至逆转疾病状态。

免疫组织化学染色

作者通过免疫组化染色验证了正常和KC角膜样本中MYRF和ATOH 7的表达水平。按照方法章节所述计算代表基因表达水平的阳性细胞面积。与上述结果一致，KC样品中MYRF的表达（2.13%）显著高于对照样品（1.62%）（p< 0.01;图13A-C）。相比之下，KC样品中ATOH 7的表达（0.904%）显著低于对照样品（1.22%）（p< 0.05;图13 D-F）。ATOH 7和MYRF两者主要在上皮中表达，但在KC组的基质中观察到更多的MYRF染色（图13 E），其中成纤维细胞占更高的比例，与图5中所示的scRNA-seq分析的结果一致。