一种特发性肺纤维化浆细胞特征基因预后模型构建方法

文献发布时间：2024-04-18 19:59:31

技术领域

本发明涉及特发肺纤维化浆细胞特征基因技术领域，特别是涉及一种特发性肺纤维化浆细胞特征基因预后模型构建方法。

背景技术

特发性肺纤维化(IPF)是一种罕见的、慢性、侵袭性和纤维化性的间质性肺炎，其特征是呼吸困难和肺功能逐渐恶化，导致死亡。虽然IPF患者的自然病史各不相同，但诊断后的中位生存期仅为2-3年，5年生存率不超过40％，尤其是IPF急性加重的短期死亡率可能高达60％或更高。由于IPF发病机制的复杂性，目前仍缺乏有效的药物来预防IPF的进展或降低死亡率。因此，迫切需要探索与IPF相关的潜在预后特征并开发治疗靶点。

反复上皮细胞损伤、成纤维细胞活化和细胞外基质沉积是IPF最基本的病理机制。越来越多的证据表明，免疫功能失调在IPF的发病机制中起着不可或缺的作用，其中浆细胞的过度适应性免疫反应可能是促进和维持IPF肺中持续炎症的关键机制。几项研究表明，在IPF肺中观察到T、B细胞、树突状细胞、巨噬细胞等细胞群体浸润的异位淋巴结构，CXCL13的异常表达诱导CD4

目前IPF尚缺乏有效的治疗药物。美国FDA批准的两种药物吡非尼酮和尼达尼布虽然有报道能够阻止疾病进展或降低死亡率，但仍有部分患者对药物应答不明显。由于免疫反应被认为是IPF的重要发病机制，专门针对基本免疫机制的新型疗法逐渐出现。其中，靶向人抗CD20抗体的利妥昔单抗被用于治疗IPF恶化患者的临床试验，但其对肺功能的改善有限，这可能是由于利妥昔单抗仅降低了患者对某些抗原的召回反应而总免疫球蛋白水平保持不变。同时，浆细胞缺乏CD20受体的表达也使得无法被利妥昔单抗消除。因此有研究者对博来霉素治疗小鼠使用蛋白酶体抑制剂硼替佐米清除浆细胞，结果发现博来霉素诱导的肺纤维化得到显著抑制。然而，硼替佐米的毒性和靶外效应十分常见，需要进一步评估其用于IPF临床治疗的可行性。另外令人兴奋的是，近年来新兴的CD19 CAR-T细胞疗法能不仅有效地消除重症及难治性系统性红斑狼疮(SLE)患者的B细胞和浆细胞，还让SLE患者实现无治疗缓解和B细胞重建。考虑到自身免疫综合征和IPF之间的许多相似之处，IPF患者可能会从CD19CAR-T细胞治疗中获益，但需要评估细胞因子释放综合征和免疫效应细胞相关神经毒性综合征等常见副作用的危害，以及CD19低表达的长寿命浆细胞的潜在残留与致病能力的关系。总之，在探索与IPF预后密切的浆细胞相关基因基础上，再挖掘IPF的可替代性治疗方式是十分有意义的。

因此，借助人工智能技术设计的一种特发性肺纤维化浆细胞特征基因预后模型，这可能有助于优化IPF的精准治疗，进一步改善IPF患者的临床预后结果。

发明内容

本发明的目的是提供一种特发性肺纤维化浆细胞特征基因预后模型构建方法，建立能够区分低风险组和高风险组患者的预后模型，可应用于特发性肺纤维化患者预后情况的辅助判断，有助于新型药物靶点的研发和新型监测方案的建立。

一种特发性肺纤维化浆细胞特征基因预后模型构建方法，包括如下步骤：

步骤1：进行单细胞转录组数据分析；

步骤2：进行组织水平转录组数据差异表达基因分析；

步骤3：进行免疫细胞浸润分析；

步骤4：进行加权基因共表达网络分析；

步骤5：进行浆细胞相关候选基因的筛选；

步骤6：进行浆细胞相关候选基因的评分；

步骤7：进行浆细胞相关候选基因的功能通路富集分析；

步骤8：进行IPF预后模型的构建及内部评价；

步骤9：进行IPF预后模型的外部评价。

可选的，步骤1中，进行单细胞转录组数据分析，具体为：

在GEO数据库中检索数据集GSE132771，将每个样本的基因-细胞表达矩阵导入R包Seurat v4.3.0.1。首先将少于50个基因或多于5000个基因或线粒体表达比例超过10％的细胞剔除。接着，使用SCTransform函数对每个样本的矩阵进行归一化处理并识别出3000个高变异基因。使用FindIntegrationAnchors和IntegrateData函数进行典型相关性分析以整合所有样本的矩阵并去除批次效应。使用RunPCA函数进行主成分分析。对于细胞类群，在FindNeighbours函数中使用了前50个主成分，并使用FindClusters函数以0.6的分辨率进行聚类。使用RunUMAP函数的UMAP算法对前50个主成分降维至二维空间。使用R包SingleRv2.2.0对细胞类群注释后再进行人工校正。使用FindAllMarkers函数(关键参数：test.use＝wilcox；min.pct＝0.1；logfc.threshold＝0.25)识别每种细胞类型的DEGs。使用FindConservedMarkers函数计算IPF组和对照组之间的各细胞簇特异性标记基因，纳入标准为平均log2FC>0.25且调整后的P值<0.05(初始P值使用Bonferroni算法校正)。使用ClusterGVis R软件包对每种细胞类型的前30个DEGs进行GO分析，并根据P值排名选择前五条注释。

可选的，步骤2中，进行组织水平转录组数据差异表达基因分析，具体为：

在GEO数据库中检索GSE150910和GSE70866数据集。对于GSE150910数据集的RNA测序数据分析，在R中使用DESeq2 v1.40.2包筛选IPF组和对照组之间的DEGs。所有样本的采集时间批次通过DESeqDataSetFromMatrix函数去除。采用Benjamini-Hochberg方法对初始P值进行调整。对于GSE70866数据集的RNA微阵列数据分析，使用Sangerbox网站进行DEG分析，具体而言使用鲁棒多阵列平均算法对原始数据进行归一化处理，并使用ComBat算法消除批次效应。使用limma v3.56.2R包筛选了GSE70866数据集中IPF组和对照组之间的DEGs。初始P值同样使用Benjamini-Hochberg方法进行调整。DEG的筛选标准为调整后的P值<0.05且|log2FC|>1。

可选的，步骤3中，进行免疫细胞浸润分析，具体为：

在GSE150910数据集中，基于LM22文件(一个白细胞特征基因矩阵)使用CIBERSORT算法并设置1000次置换来计算22种免疫细胞的浸润分数，最终从GSE150910数据集中选择每个样本的浆细胞浸润分数纳入后续分析。

可选的，步骤4中，进行加权基因共表达网络分析，具体为：

在Sangerbox网站上使用WGCNAv1.72-1R包生成GSE150910数据集的基因共表达网络。简言之，计算每个基因的中位数绝对偏差，并排除中位数绝对偏差最小的前50％的基因。为构建无向网络计算了合适的软阈值β值。然后，将加权邻接矩阵转化为拓扑重叠矩阵，以降低噪音并估计其互连性。此外，采用平均连锁层次聚类和动态树切割方法对基因模块进行检测。选择与IPF诊断和浆细胞浸润分数显著且最高正相关的基因模块纳入后续分析。在该模块中，当基因的模块成员数>0.7，基因显著性>0.2时，定义为枢纽基因。

可选的，步骤5中，进行浆细胞相关候选基因的筛选，具体为：

通过jvenn工具将浆细胞标记物分别与加权基因共表达网络分析的枢纽基因、GSE150910和GSE70866数据集的DEGs进行交集后取并集的基因作为候选基因。

可选的，步骤6中，进行浆细胞相关候选基因的评分，具体为：

在R中使用AUCell v1.22.0分析候选基因在浆细胞中是否活跃表达，基于基因集富集分析和根据候选基因的曲线下面积值构建每个细胞的基因表达排名。在基因集中表达较高比例候选基因的细胞具有较高的曲线下面积值。使用AUCell_exploreThresholds函数计算基于基因集区分激活细胞的阈值。随后，使用ggplot2 v3.4.3R包生成UMAP嵌入图，绘制每个细胞的曲线下面积值以可视化活跃的细胞簇。

可选的，步骤7中，进行浆细胞相关候选基因的功能通路富集分析，具体为：

在Sangerbox网站上通过clusterProfilerv3.14.3R包对候选基因进行GO分析，包括生物过程、分子功能和细胞成分，以及基于候选基因的京都基因和KEGG富集分析。根据P值排名选择前十个最显著的GOs和通路作为分析对象。

可选的，步骤8中，进行IPF预后模型的构建及内部评价，具体为：

将GSE70866数据集中的Freiburg队列和Siena队列合并为训练队列，并通过Sangerbox网站调整训练队列基因表达的批次效应，将Leuven队列作为独立验证队列。在训练队列中，对候选基因进行单变量Cox回归分析，以筛选具有IPF预后潜力的标志物，基于十折交叉验证框架和十种机器学习算法的101种组合，包括随机生存森林、弹性网络、套索回归、岭回归、逐步Cox、CoxBoost、偏最小二乘回归Cox、监督主成分分析、通用梯度回归模型和生存支持向量机。具体细节如下：(1)随机生存森林模型通过randomForestSRC v3.2.2R包实现。该模型的关键参数为ntree和mtry，其中ntree表示树的数量，mtry是用于在每个节点进行拆分的随机选择的变量数。ntree设置为1000，并对mtry使用网格搜索来寻找最佳节点大小。基于确定的节点大小，形成(ntree，mtry)的所有组合，并且具有最佳一致性指数的组合被识别为最优参数。(2)弹性网络、套索回归和岭回归通过glmnet v4.1-8R包实现。正则化参数λ通过10倍交叉验证确定，而L1-L2权衡参数α设置为0-1(间隔为0.1)。(3)逐步Cox模型通过survival v3.5-5R包实现。调用了基于Akaike信息准则的逐步选择算法，并将搜索步骤的方向模式分别设置为“同时”、“向后”和“向前”。(4)CoxBoost模型通过CoxBoostv1.5 R包实现。该模型用于通过逐个分量的基于似然的增强方法拟合Cox比例风险模型。对于CoxBoost模型，我们基于10倍交叉验证框架首先使用opticCoxBoostPenalty函数确定最佳惩罚参数(收缩量)。确定该参数后通过cv.CoxBoost函数来调整其他参数，即进行提升步骤的数量。最终通过CoxBoost函数来设置所选多变量Cox模型的维度。(5)plsRcox模型通过plsRcoxv1.7.7R包实现。使用cv.plsRcox函数来确定所需组件的数量，并使用plsRcox函数拟合偏最小二乘回归广义线性模型。(6)SuperPC模型通过SuperPC v1.12 R包实现。该模型是主成分分析的扩展，通过生成感兴趣的特征或变量的线性组合，捕捉数据集中最大变化的方向。superpc.cv函数基于十折交叉验证框架来估计监督主成分中的最优特征阈值。我们采用了“预验证”方法来避免将Cox模型拟合到小型验证集的问题。(7)GBM模型通过GBMv2.1.8.1R包实现。该模型基于十折交叉验证框架，通过cv.gbm函数选择具有最小交叉验证误差的树的索引。使用gbm函数来拟合广义提升回归模型。(8)生存支持向量机模型通过survivalsvm v0.0.5R包实现。回归方法在构建支持向量问题的不等式约束时考虑了截尾情况。所有模型都在验证队列中进行了评估，评价指标为Harrell一致性指数。选取一致性指数最大的模型(或模型组合)作为最优方案。对于内部评价，使用R包survminer v0.4.9计算低风险组和高风险组的最佳界值。使用R包survival v3.5-5构建Kaplan-Meier生存曲线，并用log-rank检验方法评价组间的预后差异。使用R包ggrisk v1.3计算风险评分、生存状态和基因表达水平之间的关系。

可选的，步骤9中，进行IPF预后模型的外部评价，具体为：

在Sangerbox网站上使用pROC v1.17.0.1计算最优模型的时间依赖性曲线下面积值。使用R包survival v3.5-5和forestplot v3.1.3将最优模型计算的风险评分与其他临床特征(包括年龄、性别和性别-年龄-生理学指标)进行独立和联合比较。使用R包pecv2023.04.12生成最优模型的校准曲线。为了进一步比较最优模型与基于其他预后生物标志物构建的预后模型的性能，我们收集了一系列与GSE70866数据集有关的已公开基因特征(涉及自噬、脂肪酸代谢、上皮-间充质转化免疫浸润和CX3CR1受体信号传导途径等)，分别在训练和验证队列中将其表达水平代入文献记录的公式，计算各模型的一致性指数，最终与我们构建的最优模型输出的一致性指数进行比较。

根据本发明提供的具体实施例，本发明公开了以下技术效果：本发明提供的特发性肺纤维化浆细胞特征基因预后模型构建方法，该方法包括进行单细胞转录组数据分析，进行组织水平转录组数据差异表达基因分析，进行免疫细胞浸润分析，进行加权基因共表达网络分析，进行浆细胞相关候选基因的筛选，进行浆细胞相关候选基因的评分，进行浆细胞相关候选基因的功能通路富集分析，进行IPF预后模型的构建及内部评价，进行IPF预后模型的外部评价，该方法建立能够区分低风险组和高风险组患者的预后模型，可应用于特发性肺纤维化患者预后情况的辅助判断，有助于新型药物靶点的研发和新型监测方案的建立。

附图说明

为了更清楚地说明本发明实施例或现有技术中的技术方案，下面将对实施例中所需要使用的附图作简单地介绍，显而易见地，下面描述中的附图仅仅是本发明的一些实施例，对于本领域普通技术人员来讲，在不付出创造性劳动性的前提下，还可以根据这些附图获得其他的附图。

图1为本发明实施例特发性肺纤维化浆细胞特征基因预后模型构建方法流程示意图；

图2为本发明实施例浆细胞相关基因获取示意图；

图3为本发明实施例组织水平转录组数据差异表达基因可视化示意图；

图4为本发明实施例加权基因共表达网络分析枢纽基因获取示意图；

图5为本发明实施例浆细胞相关候选基因获取示意图；

图6为本发明实施例浆细胞相关候选基因评分示意图；

图7为本发明实施例浆细胞相关候选基因通路富集分析示意图；

图8为本发明实施例特发性肺纤维化浆细胞特征基因预后模型构建及内部评价示意图；

图9为本发明实施例特发性肺纤维化浆细胞特征基因预后模型外部评价示意图。

具体实施方式

为使本发明的上述目的、特征和优点能够更加明显易懂，下面结合附图和具体实施方式对本发明作进一步详细的说明。

在本发明中，使用Seurat v4.3.0.1R包从GSE132771数据集中提取浆细胞标记基因。然后，利用DESeq2 v1.40.2R包筛选IPF组和对照组之间的DEGs。其次，利用CIBERSORT算法分析GSE150910数据集中所有样本的免疫浸润水平，探讨浆细胞与IPF之间的关系。再次，依据IPF诊断结果和浆细胞浸润水平，在Sangerbox网站上使用WGCNAv1.72-1R包生成GSE150910数据集的基因共表达网络并获得枢纽基因。接着，通过jvenn工具将浆细胞标记物分别与枢纽基因、GSE150910和GSE70866数据集的DEGs进行交集后取并集获取浆细胞相关候选基因，并进行细胞定位评分和功能通路富集分析。最后，对浆细胞相关候选基因在训练队列中使用单因素Cox回归分析以筛选具有IPF预后潜力的标志物，再基于十折交叉验证框架和十种机器学习算法的101种组合成功构建了基于IPF患者浆细胞相关特征基因的预后模型，并在内部和外部评价中显示出良好性能，可用于预测患者生存率的高低。

如图1所示，本发明实施例提供的特发性肺纤维化浆细胞特征基因预后模型构建方法，包括如下步骤：

步骤1：进行单细胞转录组数据分析；

步骤2：进行组织水平转录组数据差异表达基因分析；

步骤3：进行免疫细胞浸润分析；

步骤4：进行加权基因共表达网络分析；

步骤5：进行浆细胞相关候选基因的筛选；

步骤6：进行浆细胞相关候选基因的评分；

步骤7：进行浆细胞相关候选基因的功能通路富集分析；

步骤8：进行IPF预后模型的构建及内部评价；

步骤9：进行IPF预后模型的外部评价。

步骤1中，进行单细胞转录组数据分析，具体为：

在GEO数据库中检索数据集GSE132771，将每个样本的基因-细胞表达矩阵导入R包Seurat v4.3.0.1。首先将少于50个基因或多于5000个基因或线粒体表达比例超过10％的细胞剔除。接着，使用SCTransform函数对每个样本的矩阵进行归一化处理并识别出3000个高变异基因。使用FindIntegrationAnchors和IntegrateData函数进行典型相关性分析以整合所有样本的矩阵并去除批次效应。使用RunPCA函数进行主成分分析。对于细胞类群，在FindNeighbours函数中使用了前50个主成分，并使用FindClusters函数以0.6的分辨率进行聚类。使用RunUMAP函数的UMAP算法对前50个主成分降维至二维空间。使用R包SingleRv2.2.0对细胞类群注释后再进行人工校正。使用FindAllMarkers函数(关键参数：test.use＝wilcox；min.pct＝0.1；logfc.threshold＝0.25)识别每种细胞类型的DEGs。使用FindConservedMarkers函数计算IPF组和对照组之间的各细胞簇特异性标记基因，纳入标准为平均log2FC>0.25且调整后的P值<0.05(初始P值使用Bonferroni算法校正)。使用ClusterGVis R软件包对每种细胞类型的前30个DEGs进行基因本体论(GO)分析，并根据P值排名选择前五条注释。

最终，共保留了25043个细胞，其中11973个细胞来自三名对照者，13070个细胞来自三名IPF患者。对于细胞类群，共有30个细胞簇中被鉴定，经过注释后获得17种细胞类型(图2A)。我们共鉴定出肺泡Ⅰ型上皮细胞(标记物为AGER和CLIC5)、肺泡Ⅱ型上皮细胞(标记物为SFTPC、SFTPB和MUC1)、B细胞(标记物为MS4A1、CD79A和CD79B)、纤毛细胞(标记物为FOXJ1和CCDC78)、俱乐部细胞(标记物为SCGB3A2和SCGB1A1)、树突状细胞(标记物为CD1C和THBD)、内皮细胞(标记物为PECAM1和VWF)、成纤维细胞(标记物为表达COL1A1、PDGFRA和ELN)、巨噬细胞(标记物为CD68、MSR1和MRC1)、肥大细胞(标记物为MS4A2、CPA3和TPSAB1)、间皮细胞(标记物为MSLN和UPK3B)、单核细胞(标记物为CD14和S100A8)、自然杀伤细胞(标记物为KLRD1和NKG7)、浆细胞样树突状细胞(标记物为LILRB4、IRF8和LILRA4)，浆细胞(标记物为XBP1、CD27和SSR4)、平滑肌细胞(标记物为CNN1、ACTA2、TAGLN和RGS5)和T细胞(标记物为CD3E、CD8A和GZMK)(图2B)。然后，对每种细胞类型的前30个DEGs进行GO分析，结果表明经注释的细胞功能与既往报道一致(图2C)，表明分析结果的准确性。最后，共有97个基因被确定为浆细胞标记基因。

步骤2中，进行组织水平转录组数据差异表达基因分析，具体为：

最终，在GSE150910数据集中鉴定出1917个DEGs，其中1356个上调DEGs，561个下调DEGs；对于GSE70866数据集，筛选出378个DEGs，其中182个上调DEGs和196个下调DEGs(图3A-B)。

步骤3中，进行免疫细胞浸润分析，具体为：

步骤4中，进行加权基因共表达网络分析，具体为：

最终，加权基因共表达网络构建的β值设置为8(R

步骤5中，进行浆细胞相关候选基因的筛选，具体为：

通过jvenn工具将浆细胞标记物分别与加权基因共表达网络分析的枢纽基因、GSE150910和GSE70866数据集的DEGs进行交集后取并集的基因作为候选基因。

最终根据图5，黄色标记区内共30个基因被鉴定为浆细胞相关候选基因。

步骤6中，进行浆细胞相关候选基因的评分，具体为：

最终，所有细胞簇的曲线下面积值的阈值为0.2(图6A)。可以看出，浆细胞对候选基因的表达表现出更高的曲线下面积值(图6B-C)，说明候选基因定位分析的准确性。

步骤7中，进行浆细胞相关候选基因的功能通路富集分析，具体为：

最终，由图7可见，候选基因涉及的生物学功能与浆细胞涉及的蛋白质加工和免疫反应密切相关，再次说明候选基因定位分析的准确性。

步骤8中，进行IPF预后模型的构建及内部评价，具体为：

最终，基于临床表型和30个候选基因的表达谱，使用单变量Cox回归分析确定了8个预后基因。接着，在训练队列中，通过十折交叉验证框架将十种机器学习算法成对组合，并在验证队列中量化所有组合的一致性指数(图8A)。考虑到模型的简化性，基于CoxBoost和Enet(alpha＝0.7)的模型组合具有最高一致性指数(0.711)和最小基因子集(命名为PCRGS)的最优方案。简言之，使用CoxBoost算法筛选出了在88个提升步骤中具有非零系数的七个有价值特征(图8B)。Enet(alpha＝0.7)算法被进一步筛选为最佳模型(图8C)。基于PCRGS的表达水平，计算Enet(alpha＝0.7)模型中的回归系数加权(图8D)。然后，在训练和验证队列中，与低风险组相比，高风险患者的总生存期显著降低(P<0.05)。因此，所建立的预后模型成功地将IPF患者分为低风险或高风险组(图8E-F)。此外，与低风险组相比，在高风险组中观察到五个高表达的特征基因和两个低表达的特征基因(图8G-H)。

步骤9中，进行IPF预后模型的外部评价，具体为：

最终，基于PCRGS表达水平的风险评分升高的IPF患者的生存时间缩短。时间依赖性受试者操作特征曲线分析证实了PCRGS的预后能力，其中训练队列中的1、2、3、4和5年受试者操作特征曲线[95％置信区间]分别为0.72[0.62-0.82]、0.75[0.6-0.85]、0.76[0.65-0.88]、0.81[0.68-0.95]和0.87[0.73-1.00]；验证队列中的1年、2年、3年和4年受试者操作特征曲线[95％置信区间]分别为0.75[0.62-0.88]、0.72[0.56-0.88]、0.72[0.55-0.89]和0.71[0.41-1.00](图9A-B)。此外，PCRGS的预后能力明显高于年龄、性别和性别-年龄-生理学指标，而风险评分联合上述临床特征的预后能力又优于所有独立特征(图9C)，可进一步提高PCRGS的临床实用性。经过多变量Cox回归分析显示，在对上述临床特征进行调整后的PCRGS仍具有统计学意义(均P<0.05)，表明PCRGS是总体生存期的独立预后因素(图9D-E)。校准曲线也证明了PCRGS的预后准确性(图9F-G)。最后，PCRGS在验证队列中还显示出较其他预后生物标志物更稳健的预后性能(图9H)。总的来说，基于PCRGS表达水平的特发性肺纤维化预后模型构建是成功的，且具有临床应用价值。

本发明提供的特发性肺纤维化浆细胞特征基因预后模型构建方法，该方法包括进行单细胞转录组数据分析，进行组织水平转录组数据差异表达基因分析，进行免疫细胞浸润分析，进行加权基因共表达网络分析，进行浆细胞相关候选基因的筛选，进行浆细胞相关候选基因的评分，进行浆细胞相关候选基因的功能通路富集分析，进行IPF预后模型的构建及内部评价，进行IPF预后模型的外部评价，该方法建立能够区分低风险组和高风险组患者的预后模型，可应用于特发性肺纤维化患者预后情况的辅助判断，有助于新型药物靶点的研发和新型监测方案的建立。

本文中应用了具体个例对本发明的原理及实施方式进行了阐述，以上实施例的说明只是用于帮助理解本发明的方法及其核心思想；同时，对于本领域的一般技术人员，依据本发明的思想，在具体实施方式及应用范围上均会有改变之处。综上所述，本说明书内容不应理解为对本发明的限制。

完整全部详细技术资料下载

该技术已申请专利。仅供学习研究，如用于商业用途，请联系技术所有人。
专利发明人：
专利申请人：广州瑞能精准医学科技有限公司;

上一篇：耕地水土流失综合治理对象识别与优先序确定方法和系统
下一篇：一种柠檬酸的制备系统及其方法