N6-甲基腺苷相关lncRNAs是预测肾癌患者预后和免疫浸润的潜在生物标志物①

肖胜英 闫志广 曾福仁 邱 俊 卢义晨 朱小东

（湖南省人民医院（湖南师范大学附属第一医院肿瘤科），长沙 410016）

据美国癌症协会推测，2021 年将出现超过76 080 例肾癌（renal cell carcinoma，RCC）新病例，约13 780 例将死于该病［1］。随着诊断技术、手术和分子治疗的进步，肾癌患者的预后得到显著改善。然而，晚期肾癌患者的预后仍然很差，五年生存率约为11%［2］。据报道，分子标志物的发现和应用可为癌症患者提供预后价值［3］。

N6-甲基腺苷（N6-methyladenosine，m6A）修饰是在真核生物的mRNA、miRNA、长链非编码RNA（long non-coding RNA，lncRNA）和其他RNA 分子上发现的最常见和最丰富的修饰之一。这些修饰是可逆的，可影响RNA 分子的转录、加工、翻译和代谢［4］。m6A 修饰受m6A 调节器调控，包括“写入器”（甲基转移酶）、“读取器”（识别蛋白）和“擦除器”（去甲基化酶）［5］。m6A 调节多种细胞的关键生物过程，包括干细胞更新、细胞分化和对DNA 损伤的反应［6］。m6A 调节因子的异常表达与肿瘤发生、免疫调节异常和恶性肿瘤进展有关［6-7］。有报道称，甲基转移酶样因子14（METTL14）的表达水平与肾癌的病理和临床分期呈负相关，但与总生存期（overall survival，OS）呈正相关［8］。因此，m6A 调节的lncRNAs可能对癌症的增殖和转移至关重要［4］。研究证明，m6A 调节因子的失调与肾癌的发生有关［9］。然而，m6A 调节因子对lncRNA 的作用机制尚不明确。因此，了解肾癌中m6A 相关lncRNAs 的具体机制有助于识别肾癌的治疗靶点和与其预后相关的生物标志物。

本研究从癌症基因组图谱（TCGA）数据库中提取了14 087 个lncRNAs 表达谱，通过生物信息学和统计学分析鉴定了m6A 相关的lncRNAs。结果表明，27个m6A相关lncRNAs在肾癌中具有预后价值。此外，本研究构建了一个新的预后模型，即9 个m6A相关的lncRNAs的预后签名（m6A-LPS），该模型用于预测肾癌患者的OS。此外，本研究分析了不同亚组间肿瘤微环境（tumor microenvironment，TME）和免疫相关性的差异，以阐明免疫细胞浸润及其与临床预后的相关性。

1 资料与方法

1.1 资料 于2021 年4 月从TCGA 数据门户（https：//portal.gdc.cancer.gov/）下载肾癌患者的RNA-seq 转录组数据和相应的临床数据。共得到539 个肾癌样本和72 个邻近正常组织数据，纳入标准：①经病理证实的肾癌；②同时获得的关于mRNA表达谱和OS 的信息。最后纳入530 例肾癌患者和相应临床病理信息（年龄、生存时间、生存状态、性别、分级和TNM 分期）进一步分析。使用caret 包将530 例肾癌患者随机分配到训练集（266 例患者）或验证集（264 例患者）队列。表1 为两个队列的临床病理特征基线资料。根据已发表的文献确定了23个m6A 相关基因，包括“写入器”［METTL3、METTL14、METTL16、WTAP、VIRMA（KIA1499）、RBM15、RBM15B 和ZC3H13］、“擦除器”（FTO 和ALKBH5）和“读取器”（YTHDC1、YTHDC2、FMR1、IGF2BP1、IGF2BP2、IGF2BP3、LRPPRC、YTHDF1、YTHDF2、YTHDF3、HNRNPC、HNRNPA2B1 和RBMX），并从TCGA数据集中提取相应表达矩阵。

表1 TCGA训练和验证集中患者的临床病理特征［例（%）］Tab.1 Clinicopathological features of patients in TCGA training and validation cohorts［n（%）］

1.2 方法

1.2.1 lncRNA 注释 Genome Reference Consortium Human Build 38（GRCh38）的lncRNA 注释文件来自GENCODE 网站（https：//www.gencodegenes.org/hu‐man/），用于注释TCGA数据集中的lncRNAs。

1.2.2 生物信息学分析 首先，以Pearson相关分析挖掘数据集中的m6A 相关lncRNAs（|Pearsonr|>0.3和P<0.001）。与已发表的23 个m6A 相关基因中的一个或多个lncRNA 的表达值相关（|Pearsonr|>0.3和P<0.001）即为m6A 相关的lncRNA。进行单变量Cox 回归分析以过滤数据集中的与m6A 预后相关lncRNA。“glmnet”R 包被用于最小绝对收缩和选择算子（LASSO）Cox 回归（惩罚参数通过10 倍交叉验证估计）［10］，共发现9 个m6A 相关lncRNAs 可预测肾癌患者预后。以LASSO 回归算法获得回归系数和风险评分方程：风险评分=系数×m6A相关lncRNA 表达水平总和。在该方程的基础上，分别计算训练集和验证集中每个患者的风险评分。并以中位风险评分为分界点将患者分为高危组和低危组。使用“ConsensusClusterPlus”包阐明肾癌中m6A 相关lnc-RNAs 的生物学特性。然后，将肾癌患者分为不同亚型。使用JAVA 程序对MSigDB 的Hallmark 基因集“h.all.v6.2.symbols.gmt”进行基因集富集分析（GSEA），探索肾癌不同亚型间存活率差异的原因。两种亚型间的富集途径通过误检率<0.05 和归一化富集评分（NES）确定。采用R 包“estimate package”计算每个患者的免疫评分［11］。通过CIBERSORT 网站（https：//cibersort.stanford.edu/）得到每个样本中22 种免疫细胞类型比例，采用1 000 次排列算法，将CIBERSORTP<0.05 的样本进行亚组间的免疫浸润水平差异比较，亚组按照聚类亚型和风险评分分组。

1.3 统计学处理 R4.0.2 和SPSS24.0 软件用于统计学检验。采用Mann-WhitneyU检验比较癌组织与正常组织中m6A 相关lncRNAs 的表达水平。两个及两个以上的亚组采用学生t检验和单因素方差分析比较。应用卡方检验比较训练集和验证集中的分类变量。两组间差异比较采用log-rank 检验。Pearson 相关性检验用于分析亚型、风险评分、临床病理特征、程序性死亡配体1（PD-L1）和免疫浸润水平间的相关性。Cox 回归模型用于单变量和多变量分析，以确定风险评分与其他临床特征相结合的独立预后价值。受试者工作特征曲线（ROC）用于评估m6A 相关lncRNAs 的预后签名在三年和五年OS 的预测效率。P<0.05表示差异具有统计学意义。

2 结果

2.1 鉴定肾癌患者中m6A 相关的lncRNAs 使用“GENCODE”网站上下载的文件，基于ensemble ID，在TCGA 数据集中共识别14 087 个lncRNAs 用于后续分析。从TCGA数据集中提取了23个m6A相关基因的表达矩阵。通过Pearson相关分析寻找m6A相关的lncRNAs，得到239个与m6A显著相关的lncRNAs，随后将其与预后信息结合，进行单变量Cox 回归分析（P<0.05）。最后，在TCGA 数据集中筛出27 个与肾癌患者OS显著相关的m6A相关lncRNAs（图1A），其表达水平在肾癌和正常组织中差异很大（图1B）。AL008718.3、AC009948.2、AF117829.1、LINC01409、AL133243.3、AC084018.1、AC005253.1、AC090589.3、RUSC-1AS1、AC012615.6、RAD51-AS1、AL162586.1、AL928654.2、AC006435.2、AL136295.7、LINC00115、AC245052.4、AC114730.3、ARHGAP27P1-BPTFB1-KPNA2P3、AL135999.1、PTOV1-AS2、LINC00342 和AC00414在肾癌组织中的表达水平均高于正常邻近组织（P<0.05）。AC018752.1、COL18A1-AS1、RPL34-AS1 和SNHG10 在肾癌组织中的表达水平明显低于正常组织（P<0.05）。

图1 m6A 相关lncRNAs 的表达及其在肾癌患者预后中的作用Fig.1 Expressions of m6A-related lncRNAs and their role in prognosis of RCC patients

2.2 TCGA 数据集中m6A-LPS 的构建和验证 进一步分析530 例肾癌患者和相应的临床病理信息。为构建m6A-LPS 以预测肾癌患者的OS，基于TCGA训练集中的27 个预后相关的m6A 相关lncRNAs进行LASSO Cox 分析并生成了包含9 个m6A 相关lncRNAs 的m6A-LPS，即AC018752.1、COL18A1-AS1、RPL34-AS1、AL008718.3、AC009948.2、AL133243.3、AL162586.1、AL136295.7和LINC00342。使用LASSO算法得到的系数计算TCGA 训练集和验证集的风险评分，公式如下：风险评分=1.659 2×AL008718.3 表达水平−0.460 4×AC018752.1 表达水平−1.044 3×COL18A1-AS1表达水平+0.637 1×AC009948.2表达水平−1.8806×RPL34-AS1 表达水平+0.461 1×AL133243.3 表达水平+0.193 2×AL162586.1 表达水平+0.118×LINC00342 表达水平−0.085 8×AL136295.7 表达水平。根据中位风险评分将患者分为高风险组和低风险组。图2A、B 显示了TCGA训练集和验证集中9 个m6A-LPS 的风险评分、热图、OS、OS 状态和表达谱分布。在高危人群中，与m6A相关的高风险基因AC009948.2、LINC00342、AL008718.3、AL162586.1、AL133243.3和AL136295.7高表达，而保护性m6A相关基因如AC018752.1、COL18A1-AS1 和RPL34-AS1 在低风险人群中表达上调。Kaplan-Meier生存曲线显示，在TCGA 训练集和验证集中，具有高风险评分的肾癌患者临床结果较差（OS率低，OS时间短，P<0.000 1，图2C、D）。结果表明，m6A-LPS 对TCGA 训练集中肾癌患者的OS有较好的预测能力（三年ROC AUC=0.737，五年AUC=0.779，图3A、B）。在TCGA验证集中，9 个风险签名的三年和五年AUC 值分别为0.703 和0.724（图3C、D）。以上结果表明m6A-LPS 对肾癌患者具有稳健的OS预测能力。

图2 m6A 相关lncRNAs 风险评分及其对肾癌患者预后的影响Fig.2 m6A-related lncRNAs risk score and their influence on prognosis of RCC patients

图3 ROC曲线Fig.3 ROC curve

2.3 m6A-LPS的分层分析 为更好地评估m6A-LPS的预后能力，进行分层分析以确认m6A-LPS 是否保留其预测每个亚组OS的能力。与低风险患者相比，高风险肾癌患者在G1~2 和G3~4、Ⅰ~Ⅱ期和Ⅲ~Ⅳ期亚组的OS较差（图4A~D）。还证明m6A-LPS具有预测年龄≤65 岁或>65 岁患者（图4E、F）及不同性别患者（图4G、H）OS 的能力。这些数据表明m6A-LPS可作为肾癌患者的潜在预测因子。

2.4 m6A-LPS 是肾癌患者的独立预后因素 应用单变量和多变量Cox 分析评估m6A-LPS 是否是肾癌患者的独立预后因素。根据TCGA 数据集中的肾癌患者数据，单变量Cox 分析显示m6A-LPS 与OS 显著相关［风险比（HR）：2.261，95%置信区间（CI）：1.884~2.713，P<0.001，图4I］，而多变量Cox 分析进一步表明m6A-LPS 是OS 的独立预测因子（HR：2.058，95%CI：1.655~2.560，P<0.001，图4J）。因此，m6A-LPS 是验证集中肾癌患者OS 的独立预测因子（单变量分析HR：1.255，95%CI：1.149~1.371，P<0.001；多变量分析HR：1.206，95%CI：1.082~1.344，P<0.001，图4K、L）。以上结果表明m6A-LPS是一个独立的预后指标，可能有助于肾癌的临床预后评估。

图4 模型验证的生存曲线和多变量和单变量独立预后分析Fig.4 Survival curve for model validation and multivariate and univariate analysies of independent prog⁃nostic analysis

2.5 m6A 相关预后lncRNAs 的共识聚类与肾癌患者的特征和生存期显著相关 根据m6A 相关预后lncRNAs 的表达水平与模糊聚类检验比例的相似性，k=2 在k=2~9 间具有最好的聚类稳定性。根据m6A相关预后lncRNAs的表达水平，将530例肾癌患者分为两个亚型，即聚类1（n=313）和聚类2（n=217），见图5A。27个m6A相关预后lncRNAs在聚类1中的表达水平大多低于聚类2（图5B）。聚类1的OS比聚类2 长（P<0.001，图5C）。分析27 个m6A 相关预后lncRNAs 与PD-L1 间的相关性，结果表明，AL008718.3、AL133243.3、SNHG10、LINC00115、AF117829.1、AC084018.1、AC005253.1 和PD-L1 的表达水平呈正相关（P<0.05，图5D）。

图5 TCGA 队列中Cluster1/2 亚型肾癌的不同临床病理特征和生存率以及m6A 预后lncRNAs 与靶基因的表达和关系Fig.5 Differential clinicopathological features and survival of RCC in Cluster1/2 subtypes in TCGA cohort and expression and relationship of m6A prognostic lncRNAs and target genes

2.6 m6A 相关预后lncRNAs 的共识聚类及其22 种免疫细胞浸润水平 评估聚类1 与下调的m6A 相关预后lncRNAs 表达和聚类2 与上调的m6A 相关预后lncRNA 表达间的22 种免疫细胞浸润水平差异（图6A）。聚类1 显示幼稚B 细胞和巨噬细胞M0 浸润水平更高（图6A、B），而聚类2 中CD8 T 细胞和滤泡辅助T 细胞浸润水平更高（图6A、C、D）。随后进行GSEA 以阐明导致两个亚组间免疫细胞浸润差异的潜在调节机制。结果表明，肿瘤的恶性特征，包括哺乳动物雷帕霉素靶标（mTOR）信号（NES=1.827，标准化P=0.007）、血管内皮生长因子（VEGF）信号（NES=1.672，标准化P=0.019）和Notch 信号（NES=1.739，标准化P=0.037）在聚类2中活跃（图6E）。因此，两个亚组间免疫细胞浸润差异可能与mTOR、VEGF 和Notch 信号通路的活跃有关。

图6 免疫细胞浸润和TME分析Fig.6 Analysis of immune cell infiltration and TME

2.7 与肾癌分级、分期和免疫评分相关的预后风险评分 评估临床特征和风险评分间的关系。图7A为TCGA数据集中高风险和低风险组中9个m6A-LPS的表达水平。其中AC018752.1、COL18A1-AS1 和RPL34-AS1在高危组中的表达水平通常低于低危组。AC009948.2、LINC00342、AL008718.3、AL162586.1、AL133243.3和AL136295.7在低危组中表达水平较低。M 分期（P<0.001）、T 分期（P<0.001）、TNM 分期（P<0.001）、分级（P<0.001）、聚类亚型（P<0.001）和免疫评分（P<0.001）在高风险组和低风险组间差异有统计学意义。进一步探索风险评分、亚型、分类和免疫评分间的关系。风险评分随着组织学分级和分期的增加而增加（P<0.05，图7B、C）。聚类1的风险评分明显低于聚类2（P<0.001，图7D）。与低免疫评分组相比，高免疫评分组的风险评分较高（P<0.001，图7E）。提示肾癌患者的风险评分与分级、分期、聚类亚型和免疫评分显著相关。

图7 TCGA训练集中预后风险评分与临床病理特征和免疫评分相关Fig.7 Prognostic risk scores correlated with clinicopatho⁃logical features and immunoscore in TCGA training cohort

2.8 22 种免疫细胞的浸润丰度与风险评分的关系 分析22 种免疫细胞的浸润水平与风险评分间的关系，以评估9 个m6A-LPS 对肾癌免疫微环境的影响。风险评分与浆细胞、CD4 记忆T 细胞、CD8 T细胞、调节性T 细胞和滤泡辅助性T 细胞的浸润水平呈正相关（P<0.001，图8A~E）。风险评分和单核细胞浸润水平（P<0.001）、静息肥大细胞（P<0.001）、CD4 记忆静息T 细胞和巨噬细胞M2 呈负相关（P<0.001，图8F~I）。以上结果表明m6A-LPS 与肾癌TME有关。

图8 22种免疫细胞类型的浸润丰度与风险评分间的关系Fig.8 Relationship between risk score and infiltration abundances of 22 immune cells

3 讨论

众所周知，病理分期是临床上预测肾癌预后的关键依据［12］。然而，相同分期的患者预后不同，表明目前的分期系统在反映肾癌的异质性方面不甚准确［13］。因此，应探索潜在的预测和治疗生物标志物。本研究探讨了m6A 相关lncRNAs 对来自TCGA数据集的530 例患者的预后意义。在TCGA 数据集中共有27 个m6A 相关lncRNAs 表现出预后价值，其中9个用于建立m6A-LPS预测肾癌患者的OS。根据中位风险评分，肾癌患者被分为低风险组和高风险组。高风险组表现出较差的临床结果。高风险组和低风险组还与不同的聚类亚型、免疫评分、等级和TNM 分期相关。多变量Cox 回归分析证实m6ALPS 为OS 的独立危险因素。此外，本研究中的模型可独立于影响患者预后的其他临床预后因素。因此，该模型可应用于不同的临床亚组。

m6A 修饰mRNA 和lncRNAs 使其具有多种功能，调控mRNA 的剪接、稳定、核输出和翻译，是真核生物中最丰富的修饰。此外，与m6A 相关的lncRNAs 研究也引起众多癌症领域专家的关注［14］。ALKBH5 通过降低lncRNA 核旁斑点组装转录本1的甲基化加速胃癌的转移和侵袭［15］。DKK1 的lncRNA 激活调节因子通过与热休克蛋白家族A（Hsp70）成员1A 和Y-box 结合蛋白1 形成三元复合物稳定m6A 甲基化，从而促进头颈部鳞状细胞癌进展［16］。lncRNAs 的m6A 修饰可介导基因转录抑制、改变结构、调节lncRNAs 的稳定性影响肿瘤的发生发展［17-19］。lncRNAs 和m6A 是肾癌发生的重要调节因子［9］。目前已初步了解lncRNA 和m6A 在肾癌中的功能，然而，lncRNA 和m6A 在肾癌进展中的病理作用机制尚不明确，其预后生物标志物的研究甚少。本研究基于m6A 相关的lncRNA 构建了肾癌的独立预后模型。

本研究展示了肾癌中m6A 相关lncRNAs 的表达情况、预后价值和差异的免疫浸润丰度。筛选了9 个m6A 相关lncRNAs 构建模型以预测肾癌患者的OS。其中RPL34-AS1 是结直肠癌和胃癌的抑癌基因［20］。RPL34-AS1 在食管癌中通过下调RPL34 表达、在乳头状甲状腺癌中通过竞争性结合miR-3663-3p/RGS4、在宫颈癌中通过MDM2-p53 途径抑制细胞的增殖、迁移和侵袭［21-23］。本研究结果也提示RPL34-AS1 是肾癌的抑癌基因。COL18A1-AS1在肾透明细胞癌中是抑癌基因，并可能作为癌症诊断和治疗的预后生物标志物［24-25］。胆管癌患者中COL18A1-AS1 高表达可提高存活率［26］。这些报道为本研究结果提供支持，即COL18A1-AS1可能为肾癌的抑癌基因，有望成为其诊断和治疗的预后生物标志物。AL162586.1是一种关键的免疫相关lncRNA，其表达与N 期呈正相关，可能通过剪接途径影响前列腺癌患者预后［27］。本研究结果也认为AL162586.1可能为肾癌患者的致癌基因。此外，其他lncRNAs是首次被发现。这些lncRNAs可能是预测肾癌患者预后的潜在生物标志物。因此，本研究结果可能有助于识别与m6A 相关的lncRNAs，为深入了解其在肾癌中发生、发展的潜在作用提供见解。

进一步分析样本中不同免疫细胞的浸润水平以确定免疫细胞浸润和TME 在肾癌中的作用。与聚类2 相比，聚类1 中的患者具有更长的OS。B 细胞、幼稚细胞和巨噬细胞M0的浸润水平在聚类1中更高，而聚类2中CD8 T细胞和滤泡辅助T细胞的浸润水平更高。GSEA 结果表明，肿瘤的恶性功能特征包括mTOR、VEGF 和Notch 信号通路，在聚类2 中显著富集。ZHOU 等［28］观察到METTL14 表达下调可上调AKT 和mTOR 表达水平，与肾癌患者的不良预后相关。越来越多的证据表明，肿瘤中VEGF 水平升高可能导致适应性和先天免疫反应的抑制，血清或肿瘤的升高与肾癌患者预后不良有关［29］。mTOR 和VEGF 在致癌过程中发挥关键作用，包括细胞存活、增殖和预后［30］。Notch 信号通路影响TME 并调节癌症生物学过程，包括癌症的转移和进展［31］。Notch3 调节细胞周期进程和缺氧诱导因子2α，后者参与肾癌的增殖［32］。这些研究支持本文结果，即具有更短OS 的聚类2 显著富集mTOR、VEGF和Notch 信号通路，表明mTOR、VEGF 和Notch 信号通路的mRNA可作为m6A甲基化修饰的靶标。因此推测m6A 和mTOR、VEGF 和Notch 信号通路的甲基化共同参与了肾癌两个亚组间TME 的差异调节。随后进行风险评分与免疫细胞的相关分析，以评估其间关系。CD4 记忆T 细胞、CD8 T 细胞、Treg 和滤泡辅助T细胞与风险评分呈正相关。TME中浸润的活化CD4 记忆T 细胞、CD8 T 细胞、Treg 和滤泡辅助性T细胞的减少可能不利于肾癌患者的预后。本推论与何天基等［33］的结论大致相同，即大量CD4+T 细胞和瘤内CD8+T 细胞与短生存期和高肿瘤分级相关。TME 中Tregs 的诱导可促进肿瘤生长，且TM 中Tregs 的分化与癌症亚型患者总体生存率低有关。Tregs可抑制抗原呈递细胞的成熟、细胞因子的分泌以及细胞毒性颗粒酶和穿孔素的产生，还可促进血管生成、肿瘤生长、增殖和肿瘤向转移性疾病的转变［34-35］。以上研究结论同样支持本研究结果，因此，通过敲除或阻断Tregs、活化的CD4 记忆T 细胞和CD8 T细胞的相关信号通路或基因可能改善TME的抗肿瘤活性和免疫抑制。

总之，本研究建立了一个m6A 相关的lncRNA模型，为预测肾癌患者预后提供了生物标志物和治疗靶点，并为未来对lncRNA 的m6A修饰机制的研究提供见解。然而，本研究也存在一些缺陷和局限性，m6A 相关lncRNA 的生物学机制尚未完全阐明。因此，课题组将在后续工作中收集临床数据并扩大样本量，通过外部试验验证模型的准确性。