0   引言

甲状腺癌是人类内分泌系统和头颈部最常见的恶性肿瘤之一，其发病率在国内外均呈上升趋势^[1]。许多危险因素与甲状腺癌的发病相关，包括电离辐射、雌激素、碘摄入和遗传因素等。研究表明，手术联合放化疗以及分子靶向治疗有助于控制甲状腺癌进展，延长患者的无进展生存期，提高总体生存率^[2-4]。然而由于早期缺乏典型的临床表现，甲状腺癌很容易被患者自身忽略，从而错过最佳治疗时间^[5]。分化型甲状腺癌一般预后较好，但当其失分化为碘难治性甲状腺癌时便进展迅速，病死率较高，预后较差^[6]。乳头状甲状腺癌（papillary thyroid carcinoma, PTC）是最常见的甲状腺恶性肿瘤，约占甲状腺癌的80%^[7]。据报道，在发达国家PTC的发病率呈上升趋势^[8-10]，PTC的不良预后因素包括年龄、性别、肿瘤大小等^[7]。Grogan等^[11]对269例PTC患者进行平均27年的长期随访，结果发现超过25%的患者复发，并且在11%的病例中，PTC复发发生在治疗后20多年。因此，评估患者的疾病进展情况，寻找影响PTC患者预后的相关因素值得重视^[12]。

微小RNA（miRNA）是一类短的非编码RNA分子，仅含有19~22个核苷酸。它们是高度保守的分子，通过以完全或不完全互补的方式与靶基因的3’端非翻译区（3’-UTR）结合，参与靶基因的调节^[13]。越来越多的研究证实，miRNA广泛参与人类癌症的发生和进展，包括甲状腺癌。miRNA在肿瘤发生和发展过程中参与细胞增殖、凋亡和肿瘤侵袭等多种恶性生物学行为^[14]。虽有研究报道某些miRNA可以作为PTC的分子标志物，但目前仍缺乏敏感的miRNA组合标志物以评估PTC患者的预后。因此，本研究通过使用TCGA数据库，结合生物信息学分析对PTC组织中差异表达的miRNA进行筛选，利用Cox和Lasso回归分析出与PTC患者预后相关的miRNA，构建基于miRNA表达的预测PTC患者预后的生存模型。

1   资料与方法

1.1   原始数据的下载与处理

在TCGA官方网站上下载PTC患者的miRNA测序数据和患者的临床数据。其中癌组织514例，正常甲状腺组织59例。使用Perl5.24.3软件将原始miRNA测序数据转换成miRNA表达矩阵。在分析过程中，若有多个探针检测同一miRNA表达量，则取该miRNA表达量的平均值作为该miRNA的表达值。对于患者临床资料的分析，删除生存时间未知和生存时间为0的患者临床信息。

1.2   差异miRNA的筛选

以|log foldchange|≥2，错误发现率FDR < 0.05为筛选条件，使用R3.6.0软件中的edgeR包筛选出在甲状腺癌中差异表达的miRNA。使用ggplot2和heatmap软件程序包绘制差异基因的火山图和热图。

1.3   Cox/Lasso回归分析

结合患者的生存信息，首先对上述得到的差异miRNA进行单因素Cox回归分析，计算每个miRNA与甲状腺癌患者生存的风险比（hazard ratio, HR）和P值，以P < 0.05的标准筛选出与甲状腺癌患者预后显著相关的miRNA。将这些miRNA进一步行Lasso回归分析，目的减少基因之间共线性的影响，防止后续构建的风险模型变量过度拟合。Lasso回归使用交叉验证以确定参数，得到合适的模型。再将Lasso回归得到的miRNA进行多因素Cox回归分析，计算每个miRNA的多因素回归系数，构建风险评分方程。

1.4   风险预后模型的建立与分析

根据上述多因素Cox回归分析的结果，构建基于miRNA表达的风险评分方程risk score。按照文献报道使用公式：Risk score=β1×miRNA_1EXP+β2×miRNA_2EXP+......+βn×miRNA_nEXP ^[15-16]。式中β为相应miRNA的多因素回归系数，miRNA_EXP为相应miRNA的表达量。根据risk score数值的中位值，将PTC患者分为高风险评分组和低风险评分组。利用R3.6.0软件绘制模型预测预后的列线图，并比较两组患者之间生存的差异。利用R3.6.0软件绘制模型的ROC曲线和校准曲线以评价模型的敏感度和特异性。使用Survival ROC软件程序包计算受试者工作特征曲线（ROC）下面积（AUC）的数值。

2   结果

2.1   差异miRNA的筛选

利用|logfoldchange|≥2，FDR < 0.05为筛选条件，在PTC组织中共筛选到差异表达的miRNA75个。其中，上调表达的有70个、下调表达的有5个。根据FDR数值的排序前十位差异表达的miRNA，见表 1。图 1为差异miRNA相应的火山图和热图。图中红色表示与正常甲状腺组织相比，该基因在PTC组织中表达上调；绿色表示与正常甲状腺组织相比，该基因在PTC组织中表达下调。

表  1  PTC组织中前十位差异表达的miRNA

Table  1  Top 10 differentially-expressed miRNA in papillary thyroid carcinoma (PTC) tissues

下载: 导出CSV 

| 显示表格

图  1  PTC组织中差异表达的miRNA火山图(A)和热图(B)

Figure  1  Volcano plot(A) and heat map(B) of differentially-expressed miRNA in PTC tissues

下载: 全尺寸图片 幻灯片

2.2   Cox/Lasso回归分析

对上述差异表达的miRNA进行单因素Cox回归分析，计算相应miRNA与甲状腺癌患者的HR和P值，结果共有9个miRNA与甲状腺癌患者生存相关（P < 0.05）：hsa-mir-6730、hsa-mir-4709、hsa-mir-196a-2、hsa-mir-146b、hsa-mir-6860、hsa-mir-509-3、hsa-mir-513c、hsa-mir-515-1、hsa-mir-551b。进一步使用lasso回归对这些miRNA进行筛选，见图 2A，并使用交叉验证建立模型，见图 2B。结果共有8个miRNA（hsa-mir-6730, hsa-mir-4709, hsa-mir-196a-2, hsa-mir-146b, hsa-mir-6860, hsa-mir-509-3, hsa-mir-513c, hsa-mir-515-1）纳入分析模型。

图  2  Lasso回归分析筛选miRNA(A)和交叉验证结果(B)

Figure  2  miRNA screened by Lasso regression analysis(A) and cross validation results(B)

A: the Lasso regression model and cross validation method were utilized to screen miRNA. When the number of variables was 8, the partial likelihood deviation was the minimum, corresponding to the minimum λ value; B: the regression coefficient map of miRNA in lasso model, and double dashed lines showed the 1-fold standard error of the minimum partial likelihood deviation.

下载: 全尺寸图片 幻灯片

2.3   风险评分方程的建立

对上述8个miRNA进行多因素Cox回归分析，根据相应的回归系数，建立风险评分方程。Risk score=0.41×hsa-mir-196a-2_EXP–0.14×hsa-mir-146b_EXP–0.22×hsa-mir-4709_EXP+0.83×hsa-mir-509-3_EXP –0.03×hsa-mir-513c_EXP +0.36×hsa-mir-515-1_EXP – 0.33×hsa-mir-6730_EXP –0.63×hsa-mir-6860 _EXP。根据此方程，计算每位甲状腺癌患者的risk score数值，并根据risk score数值的中位值，将甲状腺癌患者分为高风险评分组和低风险评分组。

2.4   风险预后模型分析

根据上述结果，利用R3.6.0软件绘制基于8个miRNA的组合预测甲状腺癌患者生存的列线图，见图 3。

图  3  基于8个miRNA表达预测甲状腺癌患者生存的列线图

Figure  3  Nomogram of thyroid cancer patients' survival predicted based on eight miRNA expression

下载: 全尺寸图片 幻灯片

2.5   风险预后模型的评价

利用R3.6.0软件绘制模型的ROC曲线以评价模型的敏感度和特异性，见图 4。我们所构建的miRNA模型预测患者3年生存率和5年生存率的ROC曲线下面积数值AUC分别为0.860和0.896，这表明模型具有较好的敏感度和特异性。同时，校准曲线的结果也显示构建的模型可靠，这些结果表明我们构建的miRNA预后模型可以较准确的预测甲状腺癌患者的生存，见图 5。

图  4  风险预后模型的ROC曲线

Figure  4  ROC curve of risk prognosis model

下载: 全尺寸图片 幻灯片

图  5  预测模型的校准曲线

Figure  5  Calibration curve of prediction model

下载: 全尺寸图片 幻灯片

2.6   生存分析

利用R3.6.0软件中的survival程序包分析高低风险评分组患者的生存差异，绘制生存点图，图中绿色点表示未死亡患者，红色点表示死亡患者。随着risk score数值的升高，高风险评分组患者的死亡人数显著多于低风险评分组，说明高风险评分组的患者生存率较差，见图 6。

图  6  高低风险评分组患者的生存点图

Figure  6  Survival points diagram of patients in high and low risk groups

The green dots in the figure represented the surviving PTC patients, and the red dots represented the dead PTC patients. The dotted line represented the median value of risk score. The left side of the dotted line represented the low risk score group, and the right side of the dotted line represented the high risk score group. With the increase of risk score in PTC patients, the number of red dots increased gradually, and the number of dead PTC patients increased. It showed that the high risk group had poorer survival and higher risk of death.

下载: 全尺寸图片 幻灯片

3   讨论

甲状腺癌是最常见的内分泌恶性肿瘤。大多数甲状腺癌起源于甲状腺滤泡细胞。甲状腺癌可分为乳头状癌（75%~80%）、滤泡状癌（10%~15%）、甲状腺未分化癌（0.2%~2%）和由滤泡旁C细胞来源的髓样甲状腺癌（5%~10%）^[17]。许多研究表明，miRNA与甲状腺癌的发生和进展有关。Yue等^[18]研究表明，miR-7可能通过p21活化激酶1（PAK1）调节甲状腺癌细胞的生长、迁移和侵袭；Wang等^[19]研究表明，碘可以通过下调miR-422a、上调MAPK1的表达，促进甲状腺癌的发生；Dong等^[20]研究发现miR-141可以通过调控胰岛素受体底物2（IRS2）调控甲状腺癌细胞生长和转移，提示miR-141可以作为治疗甲状腺癌患者的潜在分子靶点；此外，miR-497也被认为是一种通过抑制脑源性神经营养因子发挥作用的甲状腺癌肿瘤抑制因子^[21]。因此，miRNA可以作为甲状腺癌诊断、治疗或者预后的标志物。

本研究中，我们通过TCGA数据库构建了基于8个miRNA表达的组合预测甲状腺癌患者预后的风险模型。在这8个miRNA中，有7个miRNA被报道与甲状腺癌或者其他肿瘤相关。MiR-146b被证明在乳头状甲状腺癌中表达上调，并且与肿瘤的侵袭有关^[22-23]；Fu等^[24]利用生物信息学的方法挖掘TCGA数据库发现hsa-mir-196a-2在甲状腺癌中表达显著上调，其与肿瘤进展、淋巴结转移和局部浸润有关，并且可以作为甲状腺癌的独立预后危险因素；miR-509的过表达可以抑制三阴性乳腺癌Hs578T细胞增殖、诱导细胞凋亡和抑制细胞侵袭^[25]；Tang等^[26]利用生物信息学分析发现，miR-6860和miR-509-3与甲状腺癌患者的总体生存率相关，可作为甲状腺癌的预后因子；miR-513c被证明与神经母细胞瘤^[27]、肝细胞癌^[28]的细胞迁移、入侵和增殖等有关；has-mir-4709被证明与结肠癌^[29]、胰腺导管腺癌^[30]患者的生存相关。

综上所述，我们构建了基于8个miRNA预测甲状腺癌患者生存的风险模型，该模型具较好的敏感度和特异性，模型中高风险评分组患者的死亡风险较大。研究指出^[31]，建立多基因预后模型比患者的临床病理指标更能提供精确的预后指导，对个体化治疗方案的选择具有重要的参考价值。我们所构建的8个miRNA模型可有效预测甲状腺癌患者的预后，但仍需要进一步使用大规模的多中心的循证医学证据加以验证。