1. 前言
恶性胸腔积液(MPE)是肺癌等恶性肿瘤发生转移的重要临床表现，其早期的精准鉴别对患者的治疗方案和预后影响较大[1]。现有的诊断技术存在着敏感度或特异度较低的局限性[2]，故继续探寻高诊断效能的靶标仍很重要。铜稳态的破坏已被证实参与肿瘤生长以及迁移，但是目前未见铜代谢相关基因(CMRGs)在MPE中的相关报道，对良恶性胸腔积液性质判断的研究尚不明确。转录组测序(RNA-Seq)作为分析肿瘤分子机制的技术日渐成熟，可结合生物信息学分析探究疾病潜在靶标。本研究着手分析良性和恶性胸腔积液的转录组测序信息，探索CMRGs在良恶性胸腔积液中的表达差异，期待可为患者的诊断寻找潜在的生物标志物，为临床精准诊断提供理论依据。
2. 材料与方法
2.1. 材料
实验对象
本实验将2023~2025年内蒙古科技大学包头医学院第一附属医院检验科收集的85例胸腔积液标本纳入研究。通过结合患者的病例资料，通过严格的纳入和排除标准进行筛选后剩余66例进行实验，包括MPE 33例和BPE 33例。本课题已通过医院医学伦理会审核。
病例纳入标准：(1) 出现胸腔积液的患者并经过影像学检查确认；(2) 患者的病理诊断明确以及临床资料完善；(3) 既往未接受化疗或放射治疗的患者。病例排除标准：(1) 临床诊断不明确的患者标本；(2) 患者的病因复杂，合并多种疾病者；(3) 已经经过手术、化疗、放射治疗等抗肿瘤治疗的患者。
病例分组依据：(1) 对照组入组标准：① 由于肺炎、肺脓肿或者结核等感染导致的胸腔积液；② 由于上腔静脉回流障碍、充血性心衰、缩窄性心包炎或低白蛋白血症等导致的漏出性积液；③ 发育性淋巴管引流异常导致的胸腔积液。(2) 恶性测定组入组标准：严格按照新版恶性胸腔积液诊断和治疗专家共识作为判断标准，再结合患者临床资料和病理结果加强诊断依据，确保入组的精确性。
2.2. 实验方法
2.2.1. RNA-Seq测序
将收集到的3例恶性胸腔积液和3例良性胸腔积液送检进行转录组测序。
2.2.2. 差异基因表达分析
使用RSEM软件对基于的表达定量计算，DESeq2软件对两组之间基因进行差异表达分析，计算每个基因的差异表达倍数和调整后的显著性P值。DEGs的筛选条件：差异表达倍数|log2FoldChange| > 1，P < 0.05，此时认为基因存在差异性。使用R包pheatmap绘制表达量聚类热图，采用R语言ggplots2软件包绘制DEGs的火山图。
2.2.3. 富集分析
使用R软件中的phyper函数，进行差异表达的GO注释及KEGG注释，当调整P值小于0.05则认为其差异显著。观察其生物学功能和富集通路现象。
2.2.4. 铜代谢相关基因获取
从分子特征数据库(MSigDB) (https://www.gsea-msigdb.org/gsea/msigdb)下载铜代谢相关基因，同时补充了铜死亡基因，用于后续分析。
2.2.5. 受试者工作特征曲线(ROC)分析
UCSC XENA (https://xenabrowser.net/datapages/)下载经Toil流程统一处理的TCGA和GTEx的TPM格式的RNAseq数据。提取肺癌对应TCGA的数据和GTEx中对应的正常组织数据，总共含有1013个肺癌和397个正常组织。使用R包pROC包进行对数据进行ROC分析，结果用ggplot2包进行可视化。
2.2.6. 相关性分析
使用Spearman相关性系数算法计算基因表达相关性，根据预设筛选标准，当基因对的相关系数绝对值超过0.9且其P值小于0.01时，即可判定两者存在统计学显著的共表达关系。使用igraph、ggplot2包绘制相关性桑基图和网络图。
2.2.7. RT-qPCR方法检测基因表达
将收集到的30例恶性胸腔积液和30例良性胸腔积液用于后续验证实验，从−80℃冰箱中取出收集的标本，重新溶解并混合，然后用新的无酶EP管分装标本，每管标本约100 μl左右，按顺序编号，放在冰上备用。提取总RNA。RNA逆转录成cDNA：冰上操作。将逆转录条件设置为：45℃，30 min；95℃，5 min进行逆转录。进行RT-qPCR实验，用96孔板加样，测定β-actin、MT1E、MT1X、MT2A和LOC105369559的mRNA表达量，每个样本设置3个复孔测。反应参数设置为：① 95℃，30 sec，1个循环；② 95℃，5 sec，60℃，30 sec，40个循环；③ 95℃，15 sec，60℃，60 sec，95℃，15 sec，1个循环。参照选择β-actin，PCR各引物序列如下表1：
Table 1. PCR primer sequences
表1. PCR引物序列表
基因名称 |
|
引物序列 |
MT1E |
F |
5'-ACCTCCTGCAAGAAGAGTGAGT-3' |
R |
5'-CCAGCCTTGAGTTCCCTCCCAAT-3 |
MT1X |
F |
5'-CTCGAAATGGACCCCAACTG-3' |
R |
5'-TGTACTGTGGCTGGGAAACG-3' |
MT2A |
F |
5'-ATGGATCCCAACTGCTCCTG-3' |
R |
5'-AGCAGCAGCTTTTCTTGCAG-3' |
LOC105369559 |
F |
5'-AGAGGGAGATTGGCCAGCTA-3' |
R |
5'-TGAAAACTCCACGTGCTGGT-3' |
β-actin |
F |
5'-TCCCTGGAGAAGAGCTACGA-3' |
R |
5'-AGCACTGTGTTGGCGTACAG-3' |
2.3. 统计方法
相关性分析采用sperman相关性系数算法计算，认为在相关系数绝对值大于0.9，P < 0.01时有统计学意义。将RT-qPCR的上机数据采用2−ΔΔCT处理及分析数据并记录，并应用GraphPad Prism9.0进行作图及两独立样本t检验估计正态分布变量的统计显著性，P < 0.05有统计学意义。
3. 结果
3.1. RNA-Seq质量统计
为了探究铜代谢相关基因在良恶性胸腔积液中的表达情况，通过RNA-Seq得到原始数据，首先需分析测序结果的可靠性。将得到的下机数据(表2)进行统计分析，其中样本A1、A2和A3为BPE组，样本B1、B2和B3为MPE组。将数据整理、过滤后，得到每个样本过滤前总的reads数分布在141.82~144.31百万条，过滤掉低质量片段后总的reads数分布在131.91~134.19百万条之间，clean reads比例分布在92.49%~93.24%之间，各样本的Q30比例分布在93.27%~94.55%之间。以上均代表测序质量较好，可以用于后续分析。
3.2. 差异表达基因的筛选
使用DESeg2软件进行良恶性胸腔积液之间的差异表达分析，筛选条件为：表达差异倍数|log2FoldChange| > 1，显著性P < 0.05，对筛选到的显著DEGs集进行统计，得到以下结果：总共测得47,137个基因，包含25,195个mRNA和19,483个lncRNA，其中共有470个DEGs，其中上调基因411个，下调基因59个。图1为差异基因的火山图，其中红色代表上调的差异基因，绿色代表下调的差异基因，灰色代表没有统计学差异的基因。
Table 2. Reads filtering statistics
表2. Reads过滤统计表
sample |
Total Raw Reads (M) |
Total Clean Reads (M) |
Clean Reads Q30 (%) |
Clean Reads Ratio (%) |
A1 |
141.82 |
131.93 |
94.46 |
93.02 |
A2 |
144.31 |
134.13 |
94.02 |
92.95 |
A3 |
141.82 |
132.23 |
93.27 |
93.24 |
B1 |
144.31 |
133.47 |
94.3 |
92.49 |
B2 |
144.31 |
134.19 |
94.55 |
92.99 |
B3 |
144.31 |
134.15 |
94.38 |
92.96 |
注：Total Raw Reads (M)：过滤前reads数(百万条)；Total Clean Reads (M)：过滤后reads数(百万条)；Clean Reads Q30 (%)：过滤后reads质量⼤于Q30⽐例；Clean Reads Ratio (%)：过滤后reads⽐例。
得出的470个DEGs中包括422个mRNA和48个lncRNA，进而，为了系统地识别差异表达基因的生物学联系并识别与肿瘤相关通路，我们对筛选出来的DEGs进行了GO和KEGG富集分析。GO富集分析从分子功能(MF)、细胞的组件作用(CC)、参与的生物学过程(BP)，这三个方面对不同组别的差异表达基因进行进一步的分析。
Figure 1. Differential gene volcano plot
图1. 差异基因火山图
3.3. 差异表达基因的富集分析
根据对良恶性胸腔积液的DEGs进行GO分析，DEGs在MF、BP和CC三个方面，共富集了2804个条目，当调整P值小于0.05时表示显著，其中DEGs在MF方面显著富集的是：钙离子结合、α-淀粉酶活性、α-淀粉酶活性(释放麦芽六糖)等功能(图2(A))；在BP方面富集的是：O-聚糖加工、双细胞紧密连接组件、表皮发育、细胞蛋白代谢等过程(图2(B))；而在CC方面显著富集的是：细胞外空间、细胞外外泌体、顶端质膜、细胞连接等结构(图2(C))。
注：(A) MF富集气泡图；(B) BP富集气泡图；(C) CC富集气泡图
Figure 2. GO enrichment bubble plot
图2. GO富集气泡图
将差异性表达的基因进行KEGG数据库通路分析，从而确定DEGs发挥的主要生物学功能。两组的DEGs进行KEGG信号通路富集分析后，共富集了263个通路，其中显著富集的通路有7个(图3)，了解到两组的DEGs显著富集在细胞粘附分子、花生四烯酸代谢、紧密连接、矿物质吸收等途径。
3.4. 差异表达基因与铜代谢相关基因的分析
从差异基因的GO和KEGG富集分析中我们发现其中涉及了细胞紧密连接或细胞粘附功能条目，有研究报道肿瘤的进展、迁移与其有着密切联系[3] [4]。同样地，我们观察到KEGG富集到了矿物质吸收通路，表明MPE与BPE中金属离子调控通路有显著差异。基于铜稳态受通路中的基因调控，同时铜代谢相关的基因(CMRGs)在肿瘤增殖、血管生成和转移中起着重要作用[5] [6]。使用分子特征数据库(MSigDB)筛选并加入铜死亡相关基因共154个CMRGs，然后把DEGs与CMRGs取交集，共得到了3个基因(图4)。分别是MT1E、MT1X和MT2A，说明这3个铜代谢相关基因在良恶性胸腔积液中的表达存在差异。
Figure 3. KEGG enrichment bubble plot
图3. KEGG富集气泡图
Figure 4. Venn diagram of DEGs and CMRGs
图4. DEGs和CMRGs维恩图
进一步探究MT1E、MT1X和MT2A在肺癌和正常人群中的受试者工作特征曲线(receiver operating characteristic curve, ROC)来观察其诊断性能，在TCGA和GTEx的数据集中进行验证，总共含有1013个肺癌组织和397个正常肺组织测序表达量信息，通过与临床常用筛查肺癌的肿瘤标记物癌胚抗原(CEACAM5) [7] [8]进行比较判断。使用R包分析得到ROC曲线(图5)以及曲线下面积(AUC)分别为MT1E (AUC = 0.694)，MT1X (AUC = 0.817) MT2A (AUC = 0.840)，CEACAM5 (AUC = 0.768)，AUC取值范围一般在0.5和1之间，AUC越接近于1，说明选择的变量在想要预测的结局的诊断效果越好。结果表明MT1E、MT1X和MT2A在肺腺癌患者的诊断中具有较好的效能，其中MT2A与MT1X的诊断效能高于CEACAM5。
注：(A) MT2A诊断肺癌的ROC曲线；(B) MT1X诊断肺癌的ROC曲线；(C) MT1E诊断肺癌的ROC曲线；(D) CEACAM5诊断肺癌的ROC曲线。横坐标X轴为1-特异性，也称为假阳性率，X轴越接近零准确率越高；纵坐标Y轴称为敏感度，也称为真阳性率(敏感度)，Y轴越大代表准确率越好。
Figure 5. ROC curve of target genes for diagnosing lung cancer
图5. 目的基因诊断肺癌的ROC曲线
3.5. 差异表达的CMRGs与免疫浸润的相关性
通过Timer数据库计算差异表达的CMRGs MT1X、MT1E和MT2A在LUAD中与B细胞、CD4+ T细胞、CD8+ T细胞、中性粒细胞、巨噬细胞和树突状细胞的浸润程度数据，并估算其与不同类型免疫细胞浸润丰度之间的关系(见图6)。可以观察到MT2A表达与中性粒细胞和树突状细胞浸润呈正相关(P < 0.05)，MT1X和MT1E与免疫细胞浸润的相关性系数较小。
3.6. 差异表达的lncRNA和铜代谢相关基因的相关性分析
根据现有研究显示，长链非编码RNA(lncRNA)在肺癌、胃癌、乳腺癌等多种恶性肿瘤的病理进程中发挥重要调控作用[9] [10]。构建了以归一化信号强度为基础的lncRNA-mRNA共表达网络模型。根据预设筛选标准，当基因对的相关系数绝对值超过0.9且其P值小于0.01时，即可判定该lncRNA与mRNA间存在统计学显著的共表达关系，绘制出共表达桑基图(图7)和网络图(图8)，共得到了82个网络节点，与CMRGs的mRNA连接关系最密切的lncRNA为LOC105369559，表明在恶性胸腔积液中LOC105369559参与CMRGs的调控，并且其扮演了一个重要的角色。
Figure 6. The correlation between CMRGs and immune cell infiltration
图6. CMRGs与免疫细胞浸润的相关性
Figure 7. Sankey diagram of co-expression correlations between differential lncRNAs and CMRGs
图7. 差异lncRNA和CMRGs共表达相关性桑基图
3.7. RT-qPCR验证关键基因
RNA-Seq的准确性是面向整体的数据，其可能会受到测序深度或者测序误差等的影响，导致对某个基因的表达量测量出现错误结果，所以后续对结果进行RT-qPCR验证可以提高筛选结果的准确性。本研究使用RT-qPCR对MT1E、MT1X、MT2A的mRNA表达水平和LOC105369559的lncRNA表达水平进行验证，并得到图9的结果。RT-qPCR得出目的基因在MPE和PE中表达水平趋势和测序一致，其中MT1E的P = 0.0252，MT1X的P = 0.0016，MT2A的P = 0.0442，LOC105369559的P = 0.0024，P < 0.05提示有统计学意义。
Figure 8. Co-expression correlation network of differential lncRNAs and CMRGs
图8. 差异lncRNA和CMRGs共表达相关性网络图
注：*代表P < 0.05，**代表P < 0.01。
Figure 9. RT-qPCR validation of key genes
图9. RT-qPCR验证关键基因
4. 讨论
通过对DEGs富集分析得到显著富集于细胞粘附分子(cell adhesion molecules, CAMs)、花生四烯酸(AA)代谢等通路。同样以上通路对调控肿瘤的进展起到关键作用。Wei等[11]发现线粒体核糖体蛋白L19 (MRPL19)的高表达通过调控细胞粘附分子通路，促进肿瘤微环境改变、上皮–间质转化(EMT)进程以及免疫逃逸，从而驱动肺癌的恶性进展使得预后不良，同时其表达与miR-26a-5p相关，阻断miR-26a-5p可能会限制肺癌细胞的增殖、侵袭和迁移[12]；GATA6通过激活囊性纤维化跨膜传导调节因子(CFTR)来抑制AA代谢途径的激活从而发挥抑制肺癌转移的作用[13]；上述研究结果发现调节肺癌中的CAMs途径和AA代谢途径可以影响肺癌的进展，本研究在MPE中观察到这些通路的显著富集，那么针对CAMs和AA通路中的关键因子开发的方案可能会有利于MPE患者的治疗和预后情况。
本研究的富集结果显示矿物质吸收通路在MPE中的富集，同样引起了我们的注意，既往该通路在MPE中的研究报道较少，随着其中CMRGs对肿瘤的调节作用渐渐引起关注，其重要性被逐步证实，成为近年来的研究焦点，探究CMRGs在MPE中的表达情况可能对MPE的诊断和患者的治疗提供新的生物标志物。我们发现其在MPE中的表达下调与前人在多种肿瘤中探究到的结果有一致和相反的情况。MT1E的下调与前列腺癌的进展相关，可能的机制是异常的MT1E启动子甲基化导致的[14]。MT1E也被证明可以影响EMT表型的变化，其过表达后可以通过增加上皮钙黏蛋白表达，降低神经钙黏蛋白和波形蛋白的表达从而抑制肝癌细胞的迁移和侵袭能力[15]。后续可以探究MPE中的MT1E甲基化水平和MT1E对EMT的影响进一步阐明其和肿瘤的作用机制。在透明细胞肾癌中MT1X表达上调，其可能是通过促进TGF-β1的过表达抑制免疫细胞的增殖，此时肿瘤细胞更具有侵袭性[16]。与本研究在MPE中观察的MT1X表达情况相反，可能是由于肿瘤的异质性导致了这一结果。也有报道称体外实验中肝癌细胞MT1X的表达下调，其过表达延缓了细胞的G1/S周期的进展并促进了细胞凋亡，又在裸鼠体内观察到其上调抑制了肿瘤生长和肺转移MT1X可能是通过抑制肝癌中的NF-κB途径的激活来诱导细胞周期停滞和细胞凋亡[17]。这与本研究探究出的MT1X在MPE中的表达情况相同。过表达的MT2A可以减少磷酸化的核心激酶MST1/2和LATS1/2进而减少下游效应因子YAP/TAZ的磷酸化，抑制Hippo信号通路从而促进细胞的增殖、迁移和抵抗凋亡等过程[18]，MT2A的上调可以使磷酸化IκB-α与CyclinD1显著降低从而抑制NF-κB通路的激活，有效减少急性髓系白血病的异常增殖[19]。可见MT1E、MT1X和MT2A在多种肿瘤类型中进行了研究，但其在MPE中的情况的尚不明确，本研究得出的MT1E、MT1X和MT2A在MPE中mRNA的表达水平相对BPE中降低的结论，为这一方向的研究空缺提供研究数据。
接着，我们通过与癌胚抗原这个常用指标的诊断效能对比，得到筛选出的3个CMRGs的判断性能，由于MPE和BPE测序数据集较少，我们观察在肺癌和正常肺组织数据集中的诊断性能来判断其性能趋势。得到的结果是相对于目前临床上的常用标志物癌胚抗原，MT2A和MT1X的AUC皆大于癌胚抗原，诊断性能更具优势，有望通过更进一步的研究与验证后应用于MPE的鉴别诊断中。
lncRNA在基因表达调控、染色质重塑和细胞周期调控上起着重要作用，近些年越来越多的相关研究表明lncRNA与癌症等疾病的关联紧密[20] [21]。为了探究本研究得到的MPE中差异表达的lncRNA是否与CMRGs存在相互影响的关系，试图找出关键lncRNA。我们通过计算其与CMRGs共表达相关性，绘制出网络图，得到连接最多的lncRNA为LOC105369559，意味LOC105369559可能在MPE中对于铜稳态的调控起关键作用。本研究中得到的CMRGs共表达枢纽lncRNA目前有报道的相关文献很少，LOC105369559对MPE以及其他肿瘤的调控机制仍待研究。
RNA-seq所得到的结果可能会由于误差出现错误，对筛选出的mRNA和lncRNA进一步通过RT-qPCR验证分析，得到的差异趋势与测序结果的方向相同，提高了结果的准确性。
综上所述，本研究通过对MPE以及BPE进行转录组学测序，DEGs分析，发现了MT1E、MT1X、MT2A和LOC105369559在MPE中表达显著下调，对PE的性质判断有提示性意义，后续可以通过研究分子机制进一步探索其作用条件。
5. 结论
细胞粘附分子途径以及花生四烯酸代谢途径与恶性胸腔积液进展相关。恶性胸腔积液差异表达的lncRNA LOC105369559与铜代谢相关基因关联性最强。差异表达的铜代谢相关基因MT1E、MT1X、MT2A和lncRNA LOC105369559在恶性胸腔积液中的表达显著下调。MT1X和MT2A在引起恶性胸腔积液的肺癌诊断中优于传统指标。
NOTES
*通讯作者。