引用本文： 王凡, 李培, 许芳, 李风森. 特发性肺间质纤维化中长链非编码 RNA 及转录因子的生物信息学挖掘. 中国呼吸与危重监护杂志, 2020, 19(6): 554-562. doi: 10.7507/1671-6205.201909064 复制
肺间质纤维化是一种以成纤维细胞/肌成纤维细胞的异常激活，以及胶原蛋白从肺泡到肺间质基质中的过度沉积为特征的持续性肺组织损伤，其中特发性肺间质纤维化（idiopathic pulmonary fibrosis，IPF）是间质性肺炎中最常见的类型，其 3 年到 5 年的中位生存率仅为 20% 到 30%[1-2]。由于 IPF 的复杂性和异质性，该病在全球范围内的发病率稳步上升，且与高龄呈正相关[3]。尽管近 5 年来 IPF 的药物治疗（如吡非尼酮、尼达尼布）取得了一定的进展，但由于 IPF 病程多变、不可预测、个体差异较大，到目前为止治疗效果并不理想[4]。IPF 的发病机制尚不明确，但公认的病理过程是由炎症、吸烟、病毒、环境等因素所致的肺泡上皮损伤，并随着肺组织结构的重构而使肺功能衰竭逐渐加重[5]。有研究表明 IPF 的发病与肺组织损伤区域的纤维细胞募集以及上皮间充质转化有关。某些遗传易感性，如Ⅱ型肺泡上皮细胞功能障碍也会增加 IPF 的易感性。此外，目前新兴的研究方法，如蛋白质组和转录组等分析可揭示部分与 IPF 发病相关的基因和通路，包括 转化生长因子-β、Wnt、基质金属蛋白酶、白细胞介素-13 等[6-8]。然而，IPF 在纤维化病理中的潜在机制和纤维化发展的详细过程仍不为人所知，因此肺间质纤维化的进程是不可预测的。
长链非编码 RNA（long non-coding RNAs，lncRNAs）是一类在哺乳动物基因组中被普遍转录，在表观遗传调控、转录调控以及转录后调控这 3 个层面发挥调控基因表达的作用，其相关研究已经受到越来越多的关注[9]，为 IPF 的潜在分子机制提供了一个新的视角[10]。然而，在 IPF 领域中，涉及转录组的相关研究对于低丰度转录本的分析结果不太准确，而且涉及不同阵列平台的转录本的研究较少等诸多因素导致假阳性结论[11-12]。因此，基于数据库挖掘的进一步生物信息学分析，不仅可以为转录调控提供一个全面的视角，还可以为转录后调控提供一个全面的视角，并可能揭示与 IPF 发病机制和进展相关的更可靠、更一致的调控因子。本研究基于 NCBI 基因表达芯片数据库（Gene Expression Omnibus，GEO）[13]的数据挖掘，进行深度生物信息学综合分析，从而获得在 IPF 中差异表达的 lncRNAs 和转录因子（transcription factors，TFs），通过 IPF 中差异表达的基因及其相关基因表达网络，为 IPF 的发病及进展的调控机制提供了一个新的视角，旨在为减轻或逆转肺纤维化提供有效治疗靶点。
1 资料及方法
1.1 数据来源及注释
GEO 是免费的芯片数据库，其中包括微阵列/基因图谱、二代测序、杂交阵列和芯片数据。本研究的所有数据均来自 GEO 中的三个数据集 GSE2052、GSE44723 和 GSE24206。GSE2052 的微阵列数据是基于 GPL1739 平台，其样本来源于 15 个 IPF 和 11 个正常肺组织；GSE44723 数据是基于 GPL570 平台，其样本来源于 10 个 IPF 和 4 个正常肺组织；GSE24206 数据是基于 GPL570 平台，其样本来源于 17 个 IPF 和 6 个正常肺组织[14-17]。
1.2 筛选 IPF 中差异表达的基因
从上述 3 个数据集下载 CEL 格式的原始数据用于后续生物信息分析。采用 R 软件（3.6.0 版本）中的 limma 包进行差异基因的分析[18]。采用 Cluster 3.0 进行无监督分层聚类分析将数据进行归一化处理以及 log2 转换，对来自多个探针组的基因进行预处理后，采用 t 检验计算 P 值，当│logFC│>1 以及 P<0.05 时作为截取值，选择差异表达的基因[19]。绘制表达基因的火山图，采用韦恩图（Venn 图）方法，获得三个数据集中共同异常表达的基因。
1.3 差异基因的生物功能富集分析
利用在线生物信息学数据库（The Database for Annotation，Visualization and Integrated Discovery，DAVID v6.7）（
1.4 构建蛋白质–蛋白质互作网络
对 GO 功能和 KEGG 通路富集分析的 DEGs 进行蛋白质–蛋白质互作网络（protein-protein interaction network，PPI 网络）分析并筛选出核心基因（Hub gene），为 IPF 的发病机制提供一个新的视角。通过基因交互检索在线数据库（STRING，11.0 版本）（
1.5 筛选关键基因
采用 Cytoscape 中的 MCODE（1.4.2 版本）插件，基于所确定基因网络的拓扑结构，对密集连接域进行聚类，从而筛选基因网络中的关键基因。Cytoscape 中的 ClueGO（2.5.1 版本）插件对所筛选的关键基因进行 GO 和通路富集功能分析。并进一步采用 Cytoscape 中的 BiNGO（3.0.3 版本）插件构建所筛选关键基因的生物通路关系网络[25-27]。采用超几何检验方法，获得 IPF 中的差异表达的 lncRNAs 及 TFs，其中 P<0.05 为差异有统计学意义。
1.6 实时荧光定量聚合酶链反应对筛选关键基因的验证
参考 2016 年《特发性肺纤维化诊断和治疗中国专家共识》[28]，排除其他已知的可能导致继发性肺间质纤维化的因素，如职业暴露、药物性肺损伤、结缔组织疾病；未行肺组织活检者，但胸部 HRCT 应呈典型普通型间质性肺炎改变。符合以上两条即可诊断为 IPF。根据以上诊断标准，选取的 16 例 IPF 病例均为 2017 年 10 月至 2019 年 3 月新疆维吾尔自治区中医院呼吸科门诊就诊及呼吸科病房收治的患者，根据公式：
，式中 
，
，采用双侧检验，
，
，
估算出两组样本量分别为 16 例，随机将健康体检者设置为对照组。IPF 组中男 10 例，女 6 例，年龄最小 56 岁，最大 80 岁，平均年龄（70.2±21.1）岁，病程最短 2 年，最长 8 年。对照组中男 8 例，女 8 例，年龄最小 45 岁，最大 68 岁，平均年龄（71.1±15.8）岁。经新疆维吾尔自治区中医院伦理委员会审查并批准，批准号 2018XE0109-1，所有参与者均签署知情同意书。
提取 IPF 患者及对照组的外周血单个核细胞，随后采用 TRlzol 法提取 RNA。使用 PrimeScriptTM RT reagent Kit 逆转录 RNA，每个样本的 cDNA 以 GAPDH 作为内参，目的和内参基因引物序列见表 1。使用 StepOnePlus™Real-time PCR 系统进行实时荧光定量聚合酶链反应（real-time fluorescence quantitative polymerase chain reaction，RT-qPCR）。结果表示为 3 次重复的均值，采用 2-∆∆ct 方法进行量化。采用 SPSS 22.0 和 graphpad prism 6.0 进行统计分析和图表绘制。计量资料以均数±标准差（
±s）表示，计数资料采用 χ2 检验比较差异。三组以上数据比较，对于方差齐性时采用方差分析，组间两两比较采用 SNK 和 LSD 法，方差不齐时采用非参数检验；两组数据比较，方差齐时采用 t 检验，方差不齐的数据资料则采用 t’ 检验。P<0.05 为差异有统计学意义。
表1
基因引物列表
2 结果
2.1 鉴定 IPF 中差异表达的基因
GEO 中三个数据集 GSE2052、GSE44723 以及 GSE24206 的原始数据经归一化处理（图 1a～c）。与正常肺组织相比，IPF 中共有 8483 个异常表达的基因，其中包括 988 个表达上调和 7495 个表达下调的基因。采用韦恩图方法，获得三个数据集中共同异常表达的 29 个基因（图 1d）。
图1
从三个 GEO 数据集（GSE2052、GSE44723 和 GSE24206）获得 DEGs
a～c. 三个数据集的 DEGs 分布的火山图，绘制上调基因（红点）和下调基因（蓝点）。没有显著变化的基因被标记为灰点。d. 从这三个数据集中使用 Venn 图识别 29 个常见的变化数据集。（http://www.ehbio.com/ImageGP/index.php/Home/）。不同的颜色区域表示不同的数据集。重叠区域表示共同表达的差异基因
2.2 DEGs 的 GO 和 KEGG 富集分析
结果如表 2 所示。前 6 位 GO 功能分析结果显示，在 IPF 中 DEGs 中显著富集的 BPs 主要集中在细胞黏附、细胞增殖调控、微管发育、骨骼系统发育等。前 6 位显著富集的 MFs 主要集中在钙离子结合、细胞因子结合、趋化因子活性、同源蛋白结合、趋化因子受体结合、微管结合等。前 6 位显著富集的 CCs 主要集中在胞外区域部分、胞外空间、胞外基质、蛋白质胞外基质、质膜部分以及质膜类固有成分等。前 6 位显著富集的 KEGG 通路主要集中在细胞外基质–受体交互、细胞因子–细胞因子受体交互、黏着斑、补体系统、谷胱甘肽代谢、苯丙氨酸代谢类。
表2
IPF 中 DEGs 的生物功能富集分析
2.3 分子模块的构建和富集分析
将所筛选出的 DEGs，采用 STRING 在线数据库构建 PPI，使用 Cytoscape 中的 ClusterONE 插件获得 18 个分子模块（图 2a）。进一步采用 DAVID 分析这 18 个分子模块所涉及基因的功能和通路富集信息（图 2b～d）。前 10 个分子模块中显著富集的 BPs 主要集中在蛋白多聚化、纤毛基底体–质膜对接、高尔基体囊泡转运、受体介导的内吞作用、Ⅰ型干扰素信号通路、染色体分离、核苷酸生物合成、链间交联修复、核苷酸切除修复、核分裂等。前 10 个分子模块中显著富集的 CCs 主要集中在泛素连接酶复合物、微管组织中心部分、微管相关复合物、MHC 蛋白复合物、染色体中心体区、鸟苷环化酶复合物、DNA 修复复合物等。前 10 个分子模块中显著富集的 MFs 主要集中在泛素蛋白转移酶活性、细胞骨架结构成分、微管运动活性、网格蛋白结合、肽抗原结合、着丝粒结合、磷酸转移酶活性、DNA 的催化活性等（表 3）。这 18 个分子模块中的 KEGG 通路富集结果集中在泛素介导的蛋白酶解、剪接体、嘌呤代谢、谷胱甘肽代谢、核糖体、细胞周期等。
图2
IPF GO 与 KEGG 分析
a. 采用 STRING，构建 PPI 网络；b～d. 采用 GO 分析确定模块基因的生物功能；e. 采用 KEGG 分析，识别模块基因的生物学通路信息
表3
IPF 中模块基因的显著 GO 富集分析
2.4 关键基因的筛选
采用 Cytoscape 中的 BiNGO 插件构建 IPF 中 DEGs 的生物网络，结果显示大部分的生物系统集中在细胞大分子代谢、核酸磷酸二酯键水解、RNA 加工等（图 3a）。采用 Cytoscape 的另一个插件 ClueGO 注释和可视化生物功能，结果显示大部分生物功能集中在谷胱甘肽代谢、核苷酸切除修复、泛素介导的蛋白酶解、细胞色素 P450 的生物代谢、核糖体类等方面（图 3b）。
图3
IPF 关键基因的鉴定与生物学功能分析
a. 利用 Cytoscape 的 BiNGO 插件，对显著 DEGs 进行生物网络构建；b. 利用 Cytoscape 的 ClueGO 插件，GO 和 KEGG 分析确定显著的差异基因；c. 利用超几何检验鉴定模块基因中的 ncRNAs；d. 利用超几何检验鉴定模块基因中的 TFs
通过超几何检验方法获得 IPF 中差异表达的关键 lncRNAs 及 TFs。在其中有 13 个模块共富集 30 个 microRNAs 和 4 个 lncRNAs，其中的 10 个分子模块中共富集 44 个 TFs（图 3c-3d）。根据富集分数，相关的非编码 RNA 主要富集在与泛素介导蛋白酶解 m1 模块，剪接体 m2 模块，细胞周期 m14 及 m18 模块，以及内吞作用 m12 模块，其中的 lncRNAs 有 MALAT1（模块连接数=14，P=7.6×10–3），FENDRR（模块连接数=18，P=2.5×10–3），RNU1-1（模块连接数=23，P=0）以及 TUG1（模块连接数=17，P=4.03×10–7）。根据富集分数而确定的 TFs，其生物功能主要聚集在与 GPR 信号通路相关 m3 模块，细胞外基质–受体交互 m4 模块，谷胱甘肽代谢 m5 模块，刺激神经配体–受体交互 m9 模块，胞内吞 m12 模块，细胞粘附 m13 模块，核苷酸切除修复 m17 模块，同源重组 m16 模块，以及细胞周期 m14、m18 模块，其中的 TFs 有 E2F1（模块连接数=5，P=3×10–4），TP53（模块连接数=6，P=2×10–4），YBX1（模块连接数=4，P=1.24×10–5），E2F4（模块连接数=3，P=2×10–4），SP1（模块连接数=7，P=4.3×10–3），BRCA1（模块连接数=3，P=2.8×10–3），CREB1（模块连接数=5，P=4.74×10–5），CIITA（模块连接数=5，P=4.19×10–7）。在所构建的分子网络中，根据节点颜色和 P 值等作为选择候选基因的依据，选择 MALAT1、FENDRR，以及 E2F1、YBX1 作为本研究的初步候选 lncRNAs 及 TFs。这些候选基因可能在细胞周期调控、甲基化、乙酰转移酶活性和细胞周期等层面发挥重要作用。
2.5 RT-qPCR 验证 IPF 中 lncRNAs 和 TFs
RT-qPCR 结果显示，lncRNA MALAT1、E2F1 在 IPF 患者外周血中高表达，而在对照组中低表达（P<0.001），YBX1 在 IPF 患者外周血中低表达，而在对照组中高表达（P<0.001），FENDRR 在两组患者表达差异不显著（P>0.05）。具体见表 4。
表4
候选基因的 RT-qPCR 结果
3 讨论
肺间质纤维化是以病因不明的限制性气道重塑为主要特征，尽管越来越多的证据表明活化的成纤维细胞可能对 IPF 的发病机制和进展有重要影响，但参与 IPF 调控的潜在分子机制仍不清楚。近 10 年来，表观遗传学涉及在没有改变 DNA 序列，而基因功能出现可逆的、可遗传的变化，并已逐渐成为研究热点。大量证据表明，包括 lncRNAs 在内的表观遗传学过程可导致限制性肺结构重构。lncRNAs 作为一种多功能的适配体分子，可以与 DNA、蛋白质以及特定的 RNA 相结合，通过调控 mRNA 的衰减、剪接、基因环化等，在调控基因表达方面发挥多种功能[29-32]。虽然高通量测序技术已经揭示出越来越多的 lncRNAs 与多种疾病相关，但 lncRNAs 参与 IPF 的分子机制仍缺乏足够的认识。lncRNAs 微阵列分析和 RNA 测序等高通量测序技术有助于探索 IPF 基因组水平的改变，揭示准确的诊断标志物和治疗方法。然而，在单个数据集的分析中可能只会显示单方面的结果，并且可能会出现一定的假阳性率[14]。因此，我们整合了来自 IPF 和正常肺组织的三个微阵列数据集的数据分析。
大量涉及 IPF 发病机制的研究发现，成纤维细胞和肌成纤维细胞过度积累，并伴有肺泡细胞外基质的沉积，是 IPF 病理过程的主要特征[33]。在这些过程中，主要的病理特征表现在肌成纤维细胞持续激活及成纤维细胞不断抵抗凋亡、α2β1 整合素及病理通路的异常激活等[34-36]。本研究整合了 GSE2052、GSE44723、GSE24206 等 3 个 mRNA 芯片数据集的原始数据，经过一系列生物信息学分析获得 18 个与 IPF 相关的 DEGs 模块，随后进行 GO 和 KEGG 富集分析等，与既往机制研究结果一致，本研究发现差异基因的生物功能主要集中在细胞粘附、细胞骨架发育、细胞外基质改变、细胞因子相互作用、多聚泛素化等方面。经过一系列生物信息分析，并通过收集 IPF 的血液标本及对照组的血液标本进行验证，最终获得与 IPF 发病相关的关键基因中 MALAT1、E2F1 在 IPF 组中高表达，而在对照组中低表达，YBX1 在 IPF 组中低表达，而在对照组中高表达。转移相关肺腺癌转录本 1（MALAT1），也称为核富集丰富转录本 2（NEAT2），涉及以下生物学功能类别：核糖核酸蛋白复合物的分子支架、基因转录调控和细胞周期调控。大量研究表明 lncRNA MALAT1 在多种实体肿瘤、感染性肺损伤、心肌或肾缺血再灌注损伤、心肌纤维化、肝纤维化、硅诱导的肺纤维化等多种疾病的发病机制和进展中均有作用[37-44]。E2F1 是 E2F 转录因子家族中的成员，参与调节细胞周期，介导细胞增殖和凋亡。相关文献检索结果仅揭示了尼达尼布治疗后 E2F1 和 miRNA 表达谱的变化，未揭示可能的 E2F1 相关分子机制[45]。YBX1 为 Y-Box 结合蛋白 1，是 Y-BOX 结合蛋白家族中的成员，参与调控转录翻译、premRNA 剪接等[46]。
高通量测序技术和分子机制的大规模研究为 IPF 的起源和发展提供了线索，但进一步阐明 IPF 的发病机制仍需进一步研究。考虑到单个微阵列数据集分析可能出现假阳性结果、样本量有限以及潜在的异质性，我们整合了来自三个数据集（GSE2052、GSE24206 和 GSE44723）的数据，所获得的 DEGs 所涉及的 IPF 发病分子机制仍需进一步探索。本研究的局限性在于，虽然通过对多个数据集的数据挖掘和详细的生物信息学分析，提取了与 IPF 相关的 DEGs，并根据分子相关网络获得 lncRNAs 和 TFs 等关键基因，虽经 RT-qPCR 技术获得在 mRNA 水平的表达验证，但后期仍需进行蛋白免疫印迹、免疫组织化学分析等相关实验来验证这些基因在 IPF 中的表达。除了通过分析三个数据集来识别与 IPF 相关的 DEGs 和关键基因外，还需要通过新的技术和生物信息学分析方法来探索更加可靠和准确的 DGEs。第二，必须在体内和体外进行进一步的功能丧失和功能获得实验证实。本研究为 IPF 的发病机制和进展提供了可靠而全面的视角，通过生物信息学分析可获得大量与 IPF 发病相关的潜在分子，并能通过后期的实验验证，为 IPF 的诊断和治疗提供更详细的依据。
利益冲突：本研究不涉及任何利益冲突。
肺间质纤维化是一种以成纤维细胞/肌成纤维细胞的异常激活，以及胶原蛋白从肺泡到肺间质基质中的过度沉积为特征的持续性肺组织损伤，其中特发性肺间质纤维化（idiopathic pulmonary fibrosis，IPF）是间质性肺炎中最常见的类型，其 3 年到 5 年的中位生存率仅为 20% 到 30%[1-2]。由于 IPF 的复杂性和异质性，该病在全球范围内的发病率稳步上升，且与高龄呈正相关[3]。尽管近 5 年来 IPF 的药物治疗（如吡非尼酮、尼达尼布）取得了一定的进展，但由于 IPF 病程多变、不可预测、个体差异较大，到目前为止治疗效果并不理想[4]。IPF 的发病机制尚不明确，但公认的病理过程是由炎症、吸烟、病毒、环境等因素所致的肺泡上皮损伤，并随着肺组织结构的重构而使肺功能衰竭逐渐加重[5]。有研究表明 IPF 的发病与肺组织损伤区域的纤维细胞募集以及上皮间充质转化有关。某些遗传易感性，如Ⅱ型肺泡上皮细胞功能障碍也会增加 IPF 的易感性。此外，目前新兴的研究方法，如蛋白质组和转录组等分析可揭示部分与 IPF 发病相关的基因和通路，包括 转化生长因子-β、Wnt、基质金属蛋白酶、白细胞介素-13 等[6-8]。然而，IPF 在纤维化病理中的潜在机制和纤维化发展的详细过程仍不为人所知，因此肺间质纤维化的进程是不可预测的。
长链非编码 RNA（long non-coding RNAs，lncRNAs）是一类在哺乳动物基因组中被普遍转录，在表观遗传调控、转录调控以及转录后调控这 3 个层面发挥调控基因表达的作用，其相关研究已经受到越来越多的关注[9]，为 IPF 的潜在分子机制提供了一个新的视角[10]。然而，在 IPF 领域中，涉及转录组的相关研究对于低丰度转录本的分析结果不太准确，而且涉及不同阵列平台的转录本的研究较少等诸多因素导致假阳性结论[11-12]。因此，基于数据库挖掘的进一步生物信息学分析，不仅可以为转录调控提供一个全面的视角，还可以为转录后调控提供一个全面的视角，并可能揭示与 IPF 发病机制和进展相关的更可靠、更一致的调控因子。本研究基于 NCBI 基因表达芯片数据库（Gene Expression Omnibus，GEO）[13]的数据挖掘，进行深度生物信息学综合分析，从而获得在 IPF 中差异表达的 lncRNAs 和转录因子（transcription factors，TFs），通过 IPF 中差异表达的基因及其相关基因表达网络，为 IPF 的发病及进展的调控机制提供了一个新的视角，旨在为减轻或逆转肺纤维化提供有效治疗靶点。
1 资料及方法
1.1 数据来源及注释
GEO 是免费的芯片数据库，其中包括微阵列/基因图谱、二代测序、杂交阵列和芯片数据。本研究的所有数据均来自 GEO 中的三个数据集 GSE2052、GSE44723 和 GSE24206。GSE2052 的微阵列数据是基于 GPL1739 平台，其样本来源于 15 个 IPF 和 11 个正常肺组织；GSE44723 数据是基于 GPL570 平台，其样本来源于 10 个 IPF 和 4 个正常肺组织；GSE24206 数据是基于 GPL570 平台，其样本来源于 17 个 IPF 和 6 个正常肺组织[14-17]。
1.2 筛选 IPF 中差异表达的基因
从上述 3 个数据集下载 CEL 格式的原始数据用于后续生物信息分析。采用 R 软件（3.6.0 版本）中的 limma 包进行差异基因的分析[18]。采用 Cluster 3.0 进行无监督分层聚类分析将数据进行归一化处理以及 log2 转换，对来自多个探针组的基因进行预处理后，采用 t 检验计算 P 值，当│logFC│>1 以及 P<0.05 时作为截取值，选择差异表达的基因[19]。绘制表达基因的火山图，采用韦恩图（Venn 图）方法，获得三个数据集中共同异常表达的基因。
1.3 差异基因的生物功能富集分析
利用在线生物信息学数据库（The Database for Annotation，Visualization and Integrated Discovery，DAVID v6.7）（
1.4 构建蛋白质–蛋白质互作网络
对 GO 功能和 KEGG 通路富集分析的 DEGs 进行蛋白质–蛋白质互作网络（protein-protein interaction network，PPI 网络）分析并筛选出核心基因（Hub gene），为 IPF 的发病机制提供一个新的视角。通过基因交互检索在线数据库（STRING，11.0 版本）（
1.5 筛选关键基因
采用 Cytoscape 中的 MCODE（1.4.2 版本）插件，基于所确定基因网络的拓扑结构，对密集连接域进行聚类，从而筛选基因网络中的关键基因。Cytoscape 中的 ClueGO（2.5.1 版本）插件对所筛选的关键基因进行 GO 和通路富集功能分析。并进一步采用 Cytoscape 中的 BiNGO（3.0.3 版本）插件构建所筛选关键基因的生物通路关系网络[25-27]。采用超几何检验方法，获得 IPF 中的差异表达的 lncRNAs 及 TFs，其中 P<0.05 为差异有统计学意义。
1.6 实时荧光定量聚合酶链反应对筛选关键基因的验证
参考 2016 年《特发性肺纤维化诊断和治疗中国专家共识》[28]，排除其他已知的可能导致继发性肺间质纤维化的因素，如职业暴露、药物性肺损伤、结缔组织疾病；未行肺组织活检者，但胸部 HRCT 应呈典型普通型间质性肺炎改变。符合以上两条即可诊断为 IPF。根据以上诊断标准，选取的 16 例 IPF 病例均为 2017 年 10 月至 2019 年 3 月新疆维吾尔自治区中医院呼吸科门诊就诊及呼吸科病房收治的患者，根据公式：，式中 ，，采用双侧检验，，， 估算出两组样本量分别为 16 例，随机将健康体检者设置为对照组。IPF 组中男 10 例，女 6 例，年龄最小 56 岁，最大 80 岁，平均年龄（70.2±21.1）岁，病程最短 2 年，最长 8 年。对照组中男 8 例，女 8 例，年龄最小 45 岁，最大 68 岁，平均年龄（71.1±15.8）岁。经新疆维吾尔自治区中医院伦理委员会审查并批准，批准号 2018XE0109-1，所有参与者均签署知情同意书。
提取 IPF 患者及对照组的外周血单个核细胞，随后采用 TRlzol 法提取 RNA。使用 PrimeScriptTM RT reagent Kit 逆转录 RNA，每个样本的 cDNA 以 GAPDH 作为内参，目的和内参基因引物序列见表 1。使用 StepOnePlus™Real-time PCR 系统进行实时荧光定量聚合酶链反应（real-time fluorescence quantitative polymerase chain reaction，RT-qPCR）。结果表示为 3 次重复的均值，采用 2-∆∆ct 方法进行量化。采用 SPSS 22.0 和 graphpad prism 6.0 进行统计分析和图表绘制。计量资料以均数±标准差（±s）表示，计数资料采用 χ2 检验比较差异。三组以上数据比较，对于方差齐性时采用方差分析，组间两两比较采用 SNK 和 LSD 法，方差不齐时采用非参数检验；两组数据比较，方差齐时采用 t 检验，方差不齐的数据资料则采用 t’ 检验。P<0.05 为差异有统计学意义。
表1
基因引物列表
2 结果
2.1 鉴定 IPF 中差异表达的基因
GEO 中三个数据集 GSE2052、GSE44723 以及 GSE24206 的原始数据经归一化处理（图 1a～c）。与正常肺组织相比，IPF 中共有 8483 个异常表达的基因，其中包括 988 个表达上调和 7495 个表达下调的基因。采用韦恩图方法，获得三个数据集中共同异常表达的 29 个基因（图 1d）。
图1
从三个 GEO 数据集（GSE2052、GSE44723 和 GSE24206）获得 DEGs
a～c. 三个数据集的 DEGs 分布的火山图，绘制上调基因（红点）和下调基因（蓝点）。没有显著变化的基因被标记为灰点。d. 从这三个数据集中使用 Venn 图识别 29 个常见的变化数据集。（http://www.ehbio.com/ImageGP/index.php/Home/）。不同的颜色区域表示不同的数据集。重叠区域表示共同表达的差异基因
2.2 DEGs 的 GO 和 KEGG 富集分析
结果如表 2 所示。前 6 位 GO 功能分析结果显示，在 IPF 中 DEGs 中显著富集的 BPs 主要集中在细胞黏附、细胞增殖调控、微管发育、骨骼系统发育等。前 6 位显著富集的 MFs 主要集中在钙离子结合、细胞因子结合、趋化因子活性、同源蛋白结合、趋化因子受体结合、微管结合等。前 6 位显著富集的 CCs 主要集中在胞外区域部分、胞外空间、胞外基质、蛋白质胞外基质、质膜部分以及质膜类固有成分等。前 6 位显著富集的 KEGG 通路主要集中在细胞外基质–受体交互、细胞因子–细胞因子受体交互、黏着斑、补体系统、谷胱甘肽代谢、苯丙氨酸代谢类。
表2
IPF 中 DEGs 的生物功能富集分析
2.3 分子模块的构建和富集分析
将所筛选出的 DEGs，采用 STRING 在线数据库构建 PPI，使用 Cytoscape 中的 ClusterONE 插件获得 18 个分子模块（图 2a）。进一步采用 DAVID 分析这 18 个分子模块所涉及基因的功能和通路富集信息（图 2b～d）。前 10 个分子模块中显著富集的 BPs 主要集中在蛋白多聚化、纤毛基底体–质膜对接、高尔基体囊泡转运、受体介导的内吞作用、Ⅰ型干扰素信号通路、染色体分离、核苷酸生物合成、链间交联修复、核苷酸切除修复、核分裂等。前 10 个分子模块中显著富集的 CCs 主要集中在泛素连接酶复合物、微管组织中心部分、微管相关复合物、MHC 蛋白复合物、染色体中心体区、鸟苷环化酶复合物、DNA 修复复合物等。前 10 个分子模块中显著富集的 MFs 主要集中在泛素蛋白转移酶活性、细胞骨架结构成分、微管运动活性、网格蛋白结合、肽抗原结合、着丝粒结合、磷酸转移酶活性、DNA 的催化活性等（表 3）。这 18 个分子模块中的 KEGG 通路富集结果集中在泛素介导的蛋白酶解、剪接体、嘌呤代谢、谷胱甘肽代谢、核糖体、细胞周期等。
图2
IPF GO 与 KEGG 分析
a. 采用 STRING，构建 PPI 网络；b～d. 采用 GO 分析确定模块基因的生物功能；e. 采用 KEGG 分析，识别模块基因的生物学通路信息
表3
IPF 中模块基因的显著 GO 富集分析
2.4 关键基因的筛选
采用 Cytoscape 中的 BiNGO 插件构建 IPF 中 DEGs 的生物网络，结果显示大部分的生物系统集中在细胞大分子代谢、核酸磷酸二酯键水解、RNA 加工等（图 3a）。采用 Cytoscape 的另一个插件 ClueGO 注释和可视化生物功能，结果显示大部分生物功能集中在谷胱甘肽代谢、核苷酸切除修复、泛素介导的蛋白酶解、细胞色素 P450 的生物代谢、核糖体类等方面（图 3b）。
图3
IPF 关键基因的鉴定与生物学功能分析
a. 利用 Cytoscape 的 BiNGO 插件，对显著 DEGs 进行生物网络构建；b. 利用 Cytoscape 的 ClueGO 插件，GO 和 KEGG 分析确定显著的差异基因；c. 利用超几何检验鉴定模块基因中的 ncRNAs；d. 利用超几何检验鉴定模块基因中的 TFs
通过超几何检验方法获得 IPF 中差异表达的关键 lncRNAs 及 TFs。在其中有 13 个模块共富集 30 个 microRNAs 和 4 个 lncRNAs，其中的 10 个分子模块中共富集 44 个 TFs（图 3c-3d）。根据富集分数，相关的非编码 RNA 主要富集在与泛素介导蛋白酶解 m1 模块，剪接体 m2 模块，细胞周期 m14 及 m18 模块，以及内吞作用 m12 模块，其中的 lncRNAs 有 MALAT1（模块连接数=14，P=7.6×10–3），FENDRR（模块连接数=18，P=2.5×10–3），RNU1-1（模块连接数=23，P=0）以及 TUG1（模块连接数=17，P=4.03×10–7）。根据富集分数而确定的 TFs，其生物功能主要聚集在与 GPR 信号通路相关 m3 模块，细胞外基质–受体交互 m4 模块，谷胱甘肽代谢 m5 模块，刺激神经配体–受体交互 m9 模块，胞内吞 m12 模块，细胞粘附 m13 模块，核苷酸切除修复 m17 模块，同源重组 m16 模块，以及细胞周期 m14、m18 模块，其中的 TFs 有 E2F1（模块连接数=5，P=3×10–4），TP53（模块连接数=6，P=2×10–4），YBX1（模块连接数=4，P=1.24×10–5），E2F4（模块连接数=3，P=2×10–4），SP1（模块连接数=7，P=4.3×10–3），BRCA1（模块连接数=3，P=2.8×10–3），CREB1（模块连接数=5，P=4.74×10–5），CIITA（模块连接数=5，P=4.19×10–7）。在所构建的分子网络中，根据节点颜色和 P 值等作为选择候选基因的依据，选择 MALAT1、FENDRR，以及 E2F1、YBX1 作为本研究的初步候选 lncRNAs 及 TFs。这些候选基因可能在细胞周期调控、甲基化、乙酰转移酶活性和细胞周期等层面发挥重要作用。
2.5 RT-qPCR 验证 IPF 中 lncRNAs 和 TFs
RT-qPCR 结果显示，lncRNA MALAT1、E2F1 在 IPF 患者外周血中高表达，而在对照组中低表达（P<0.001），YBX1 在 IPF 患者外周血中低表达，而在对照组中高表达（P<0.001），FENDRR 在两组患者表达差异不显著（P>0.05）。具体见表 4。
表4
候选基因的 RT-qPCR 结果
3 讨论
肺间质纤维化是以病因不明的限制性气道重塑为主要特征，尽管越来越多的证据表明活化的成纤维细胞可能对 IPF 的发病机制和进展有重要影响，但参与 IPF 调控的潜在分子机制仍不清楚。近 10 年来，表观遗传学涉及在没有改变 DNA 序列，而基因功能出现可逆的、可遗传的变化，并已逐渐成为研究热点。大量证据表明，包括 lncRNAs 在内的表观遗传学过程可导致限制性肺结构重构。lncRNAs 作为一种多功能的适配体分子，可以与 DNA、蛋白质以及特定的 RNA 相结合，通过调控 mRNA 的衰减、剪接、基因环化等，在调控基因表达方面发挥多种功能[29-32]。虽然高通量测序技术已经揭示出越来越多的 lncRNAs 与多种疾病相关，但 lncRNAs 参与 IPF 的分子机制仍缺乏足够的认识。lncRNAs 微阵列分析和 RNA 测序等高通量测序技术有助于探索 IPF 基因组水平的改变，揭示准确的诊断标志物和治疗方法。然而，在单个数据集的分析中可能只会显示单方面的结果，并且可能会出现一定的假阳性率[14]。因此，我们整合了来自 IPF 和正常肺组织的三个微阵列数据集的数据分析。
大量涉及 IPF 发病机制的研究发现，成纤维细胞和肌成纤维细胞过度积累，并伴有肺泡细胞外基质的沉积，是 IPF 病理过程的主要特征[33]。在这些过程中，主要的病理特征表现在肌成纤维细胞持续激活及成纤维细胞不断抵抗凋亡、α2β1 整合素及病理通路的异常激活等[34-36]。本研究整合了 GSE2052、GSE44723、GSE24206 等 3 个 mRNA 芯片数据集的原始数据，经过一系列生物信息学分析获得 18 个与 IPF 相关的 DEGs 模块，随后进行 GO 和 KEGG 富集分析等，与既往机制研究结果一致，本研究发现差异基因的生物功能主要集中在细胞粘附、细胞骨架发育、细胞外基质改变、细胞因子相互作用、多聚泛素化等方面。经过一系列生物信息分析，并通过收集 IPF 的血液标本及对照组的血液标本进行验证，最终获得与 IPF 发病相关的关键基因中 MALAT1、E2F1 在 IPF 组中高表达，而在对照组中低表达，YBX1 在 IPF 组中低表达，而在对照组中高表达。转移相关肺腺癌转录本 1（MALAT1），也称为核富集丰富转录本 2（NEAT2），涉及以下生物学功能类别：核糖核酸蛋白复合物的分子支架、基因转录调控和细胞周期调控。大量研究表明 lncRNA MALAT1 在多种实体肿瘤、感染性肺损伤、心肌或肾缺血再灌注损伤、心肌纤维化、肝纤维化、硅诱导的肺纤维化等多种疾病的发病机制和进展中均有作用[37-44]。E2F1 是 E2F 转录因子家族中的成员，参与调节细胞周期，介导细胞增殖和凋亡。相关文献检索结果仅揭示了尼达尼布治疗后 E2F1 和 miRNA 表达谱的变化，未揭示可能的 E2F1 相关分子机制[45]。YBX1 为 Y-Box 结合蛋白 1，是 Y-BOX 结合蛋白家族中的成员，参与调控转录翻译、premRNA 剪接等[46]。
高通量测序技术和分子机制的大规模研究为 IPF 的起源和发展提供了线索，但进一步阐明 IPF 的发病机制仍需进一步研究。考虑到单个微阵列数据集分析可能出现假阳性结果、样本量有限以及潜在的异质性，我们整合了来自三个数据集（GSE2052、GSE24206 和 GSE44723）的数据，所获得的 DEGs 所涉及的 IPF 发病分子机制仍需进一步探索。本研究的局限性在于，虽然通过对多个数据集的数据挖掘和详细的生物信息学分析，提取了与 IPF 相关的 DEGs，并根据分子相关网络获得 lncRNAs 和 TFs 等关键基因，虽经 RT-qPCR 技术获得在 mRNA 水平的表达验证，但后期仍需进行蛋白免疫印迹、免疫组织化学分析等相关实验来验证这些基因在 IPF 中的表达。除了通过分析三个数据集来识别与 IPF 相关的 DEGs 和关键基因外，还需要通过新的技术和生物信息学分析方法来探索更加可靠和准确的 DGEs。第二，必须在体内和体外进行进一步的功能丧失和功能获得实验证实。本研究为 IPF 的发病机制和进展提供了可靠而全面的视角，通过生物信息学分析可获得大量与 IPF 发病相关的潜在分子，并能通过后期的实验验证，为 IPF 的诊断和治疗提供更详细的依据。
利益冲突：本研究不涉及任何利益冲突。
