不同病理类型甲状腺癌共同调控网络的转录组学整合研究△

谢丽君

昆明医科大学第一附属医院核医学科，昆明 653002

甲状腺癌是全球常见的恶性肿瘤，甲状腺乳头状癌（papillary thyroid carcinoma，PTC）是主要的甲状腺癌类型，治愈率超过90%；相反，滤泡状甲状腺癌（follicular thyroid carcinoma，FTC）易发生血行转移且复发率高，甲状腺未分化癌（anaplastic thyroid carcinoma，ATC）、甲状腺髓样癌（medullary thyroid carcinoma，MTC）是最具侵袭性的内分泌恶性肿瘤，ATC 和MTC 对放射性碘消融和传统化疗无效。因此，识别PTC、ATC、MTC 和FTC 间主调节因子的差异和共同调控的信号网络可能对甲状腺癌的治疗有重要意义。

近年来，多项研究阐明了不同病理类型甲状腺癌的遗传复杂性。Yoo 等[1]比较了12 例晚期分化型甲状腺癌（differentiated thyroid carcinoma，DTC）和13 例ATC 患者的mRNA 表达情况发现，异常的甘油脂类代谢、脂肪酸代谢、CDKN2A基因表达与甲状腺癌的侵袭性明显相关。Nicolson 等[2]对17例FTC 和4 例正常甲状腺组织的转录组进行测序发现，组蛋白去乙酰化酶1（histone deacetylase 1，HDAC1）是一种能在肿瘤侵袭前发挥增强表达作用的组蛋白修饰物。尽管既往已取得了重要进展，但PTC、ATC、MTC 和FTC 间的共同调控网络和特异性调控因子仍然未知。

本研究从GSE27155、GSE29265 和GSE53157共3 个数据集的组织样本中筛选基因表达谱，共包括9 例ATC、43 例PTC、2 例MTC 和4 例FTC 患者的数据，通过PTC 分别与FTC、ATC 和MTC 比较的交集来确定差异表达基因（differential expressed gene，DEG）；此外，采用基因富集分析注释4 种不同病理类型甲状腺癌共同特异基因的生物学功能和信号通路；构建转录因子调控网络以确定甲状腺癌的共同主调控因子，并分析其与甲状腺癌患者预后的关系，现报道如下。

1 资料与方法

1.1 DEG 筛选

通过基因表达综合数据库（Gene Expression Omnibus，GEO）进行DEG 筛选，在GEO 数据库中以“thyroid cancer”“Homo sapiens”“tissue”为关键词搜索基因序列，排除微小RNA、线粒体DNA、非体内甲状腺癌组织取材标本等相关数据，最终选取GSE27155、GSE29265 和GSE53157 共3 个数据集，其中GSE27155 的数据来自26 例PTC 和2 例MTC 患者，GSE29265 的数据来自10 例PTC 和9 例ATC 患者，GSE53157 的数据来自7 例PTC 和4 例FTC 患者。利用基因芯片差异分析软件GEO2R（https://www.ncbi.nlm.nih.gov/geo/geo2r/）分析各个数据集，选取P＜0.05、|logFC|＞1 的基因为候选差异基因，最终将3 个数据集中选取的差异基因取交集，筛选出DEG。利用微生信在线分析软件（http://www.bioinformatics.com.cn/）对筛选出的DEG 进一步分析以筛选出共同表达的差异基因。

1.2 基因本体（gene ontology，GO）功能注释和京都基因与基因组百科全书（Kyoto Encyclopedia of Genes and Genomes，KEGG）通路富集分析

利用2021 更新版生物学信息注释及可视化数据库（the database for annotation,visualization and integrated discovery，DAVID）（https://david.ncifcrf.gov/tools.jsp）中在线分析工具以人源基因为背景进行GO 功能注释和KEGG 通路富集分析，设定P＜0.05。

1.3 蛋白相互作用网络分析

利用String 数据库（https://cn.string-db.org/）构建蛋白质相互作用网络，并采用Cytoscape 3.7.2 进行可视化分析，通过CytoHubba 插件对构建的蛋白质相互作用网络进行关联度分析。通过分析网络结构，根据关联积分值，获得整个网络中可能形成的蛋白簇和关键节点蛋白，并在Cytoscape 中进行可视化。关键节点筛选标准为最大团中心性（maximal clique centrality，MCC）算法，同时采用Mcode 插件进行蛋白簇和关键节点蛋白分析，随后采用ClueGO 插件（http://apps.cytoscape.org/apps/cluego）对CytoHubba 生成的模块进行通路相互作用网络构建。

1.4 转录因子的预测

通过Cytoscape 3.7.2 软件中的iRegulon 插件对1.3 中获得的关键节点蛋白进行转录因子的预测。预测转录因子时，参数被设置为默认值，此外，选择标准化富集分数（normalized enrichment score，NES）＞5.0 的转录因子作为重要因子的调控网络。

1.5 不同病理类型甲状腺癌共同转录因子的Kaplan-Meier 分析

为评估1.4 中获得的转录因子对甲状腺癌患者预后的预测价值，使用Kaplan-Meier plotte 计算表达值的最佳临界值[3]，将来自癌症基因组图谱（The Cancer Genome Atlas，TCGA）（https://portal.gdc.cancer.gov/）的353 例甲状腺癌组织分为高表达组与低表达组，进行Kaplan-Meier 生存分析（www.kmplot.com），组间比较采用Log-rank 检验。对于Kaplan-Meier 曲线，P值和95%CI 的风险比（hazard ratio，HR）通过Log-rank 检验和多因素Cox 回归分析得出，以P＜0.05 为差异有统计学意义，其中不同组间HR 代表高表达组相对于低表达组的风险系数，若HR＞1 代表该基因为危险因素，若HR＜1 则代表该基因为保护因素。

2 结果

2.1 不同病理类型甲状腺癌DEG 的筛选

通过对GSE27155、GSE29265 和GSE53157 共3个数据集中PTC 与MTC、PTC 与ATC、PTC 与FTC比较分析，分别确定了2394、1104 和2805 个独立的DEG，共同交集筛选获得191 个共同表达的DEG，其中上调基因143 个，下调基因48 个。（图1）

图1 不同病理类型甲状腺癌的DEG筛选图

2.2 GO 功能注释分析

选取2.1 中所筛选出的共同表达的DEG，利用DAVID 在线分析工具共得出96 项重要富集的基因功能，其中生物进程主要富集于细胞形态调控、细胞黏附的负调节、细胞连接组件的正向调节等；细胞成分主要富集于细胞前缘、丝状膜、外侧质膜等；分子功能主要富集于其他取代磷酸基团的磷酸转移酶活性、Poly（A）结合、细胞间黏附介质活性等。（图2）

图2 191个共同表达的DEG的GO功能注释分析

2.3 KEGG 通路富集分析

选取2.1 中所筛选出的共同表达的DEG，利用DAVID 在线分析工具共得出重要富集的5 条KEGG 通路，主要包括碱基切除修复、膦酸盐和次膦酸盐代谢、致病性大肠杆菌感染、肿瘤中的蛋白聚糖、细胞黏附分子通路。涉及的主要基因包括多聚二磷酸腺苷核糖聚合酶4[poly（ADP-ribose）polymerase 4，PARP4]、8-氧鸟嘌呤DNA 糖苷酶（8-oxoguanine DNA glycosylase，OGG1）、DNA 聚合酶δ2 辅助亚基（DNA polymerase delta 2,accessory subunit，POLD2）、胆碱磷酸转移酶1（choline phosphotransferase 1，CHPT1）、胆碱/乙醇胺磷酸转移酶1 （choline/ethanolamine phosphotransferase 1，CEPT1）、连接蛋白（claudin，CLDN）3、肌球蛋白Ⅹ（myosin Ⅹ，MYO10）、FAS、肌球蛋白重链10（myosin heavy chain 10，MYH10）、CLDN1、促分裂原活化的蛋白激酶12（mitogen-activated protein kinase 12，MAPK12）、erb-b2 受体酪氨酸激酶3（erb-b2 receptor tyrosine kinase 3，ERBB3）、多配体蛋白聚糖4（syndecan 4，SDC4）、整合素β3（integrin subunit beta 3，ITGB3）、MET、钙黏蛋白3（cadherin 3，CDH3）和蛋白酪氨酸磷酸酶受体M（protein tyrosine phosphatase receptor type M，PTPRM）。（图3）

图3 191个共同表达的DEG的KEGG通路富集分析

2.4 蛋白相互作用网络可视化分析

基于前述的DEG 富集通路的蛋白构建相互作用网络，对蛋白-蛋白相互作用进行可视化分析，结果显示，甲状腺癌共同表达的DEG 形成的蛋白相互作用网络中含有24个下调基因及97个上调基因，并共同富集于重要的蛋白-蛋白相互作用通路上。

CytoHubba 插件的MCC 分析显示有24 个节点、26 条边，主要富集在蛋白质二硫键异构酶活性、鸟苷酸交换因子活性的负调节、磷脂酸磷酸酶活性的调节、富含脯氨酸区域结合、核纤层蛋白结合、磷脂酰肌醇3-激酶催化亚单位结合、SUMO 接合酶活性、细胞骨架-核膜锚定活性，关键基因包括凝溶胶蛋白（gelsolin，GSN）、RAS p21 蛋白激活因子（RAS p21 protein activator 1，RASA1）、真核生物翻译起始因子2α亚单位（eukaryotic translation initiation factor 2 subunit alpha，EIF2S1）、UBX 结构域蛋白7（UBX domain protein 7，UBXN7）、ITGB3、Yes 相关转录调控因子1（Yes1 associated transcriptional regulator 1，YAP1）、脂蛋白1（lipin 1，LPIN1）、核膜含血影蛋白重复蛋白2（spectrin repeat containing nuclear envelope protein 2，SYNE2）、泛素结合酶E2C 结合蛋白E2I（ubiquitin conjugating enzyme E2I，UBE2I）、肌营养不良蛋白（dystrophin，DMD）。（图4A）

Mcode 模块1 分析包括6 个节点、6 条边，主要富集在单酚单加氧酶活性的调节和非别肽酶活性的正调节信号通路，关键基因包括肿瘤关联钙信号转导因子2（tumor associated calcium signal transducer 2，TACSTD2）、CLDN3、CLDN1、CDH3。（图4B）

Mcode 模块2 分析包括27 个节点、53 条边，主要富集于二酰甘油胆碱转移酶活性、磷脂酰肌醇转移酶活性、富含G 链的端粒DNA 结合、错配修复、细胞骨架-核膜锚定活性、磷脂酰磷酸酶活性调节通路，关键基因包括磷脂酸胞苷酸转移酶1（CDP-diacylglycerol synthase 1，CDS1）、CEPT1、LPIN1、CHPT1、SYNE2、端粒重复序列结合因子2（telomeric repeat binding factor 2，TERF2）、端粒重复序列结合因子2 结合蛋白（TERF2 interacting protein，TERF2IP）、POLD2。（图4C）

Mcode 模块3 分析包括9 个节点、10 条边，主要富集于磷脂酰肌醇-5-磷酸结合、磷脂酰肌醇磷酸激酶活性通路，关键基因包括磷脂酰肌醇-4-磷酸3-激酶催化亚单位2α（phosphatidylinositol-4-phosphate 3- kinase catalytic subunit type 2 alpha，PIK3C2A）、分拣微管连接蛋白5（sorting nexin 5，SNX5）、肿瘤蛋白D52 样蛋白1（tumor protein D52 like 1，TPD52L1）。（图4D）

图4 蛋白相互作用网络的模块分析和基因功能注释图

2.5 主要转录因子和网络预测

使用iRegulon 插件对结果2.4 中的25 个关键基因进行转录因子预测，结果显示，有14 个重要的转录因子高度富集于蛋白网络中，包括EBF 转录因子1（EBF transcription factor 1，EBF1）、髓样锌指蛋白1（myeloid zinc finger 1，MZF1）、E2F 转录因子1（E2F transcription factor 1，E2F1）、信号转导与转录激活因子（signal transduction and activator of transcription，STAT）6、核因子κB 亚基（nuclear factor kappa B subunit，NFKB）1、NFKB2、RELA、churchill结构域蛋白1（churchill domain containing 1，CHURC1）、衔接因子相关蛋白复合体3β1（adaptor related protein complex 3 subunit beta 1，AP3B1）、含锌指和BTB 结构域7A（zinc finger and BTB domain containing 7A，ZBTB7A）、Ikaros 家族锌指蛋白2（IKAROS family zinc finger 2，IKZF2）、上游转录因子2 与c-fos 相互作用（upstream transcription factor 2,c-fos interacting，USF2）、STAT1、ovo 样锌指蛋白2（ovo like zinc finger 2，OVOL2），NES=5.119。对应重要的靶基因包括细胞分裂周期14A（cell division cycle 14A，CDC14A）、Rho GTP 酶活化蛋白44（Rho GTPase activating protein 44，ARHGAP44）、倒置形成素蛋白2（inverted formin 2，INF2）、磷酸酶和肌动蛋白调控因子2（phosphatase and actin regulator 2，PHACTR2）、碳水化合物磺基转 移 酶 15（carbohydrate sulfotransferase 15，CHST15）、CD47、UBE2I。

2.6 关键转录因子与甲状腺癌患者预后的关系

评估2.5 中共同转录因子与甲状腺癌患者预后的关系，以掌握甲状腺癌的重要调节因子，使用Kaplan-Meier plotter 法计算上述基因高表达组与低表达组甲状腺癌患者的无复发生存率（relapse-free survival，RFS），结果显示，E2F1（HR=6.51，P=0.003）、NKFB1（HR=4.13，P=0.000）和NFKB2（HR=2.35，P=0.029）均是甲状腺癌患者预后的危险因素，IKZF2（HR=0.27，P=0.023）和MZF1（HR=0.31，P=0.002）均是甲状腺癌患者预后的保护因素。（图5）

图5 不同转录因子表达情况甲状腺癌患者的RFS曲线

3 讨论

不管是预后良好的PTC，还是预后极差的ATC和MTC，以及易发生血行转移或远处转移的FTC，诊疗的相关问题已引起临床医师和研究者的高度关注。既往学者一直利用不同的技术关注不同病理类型甲状腺癌的基因组或转录组改变[4]。但由于ATC 和MTC 的罕见性和异质性常导致不同研究结果出现误差。本文从3 个独立数据集中生成了一个大数据集以探索PTC、FTC、ATC 及MTC 间调控网络的差异和共同存在的重要转录因子。

对FTC、ATC 及MTC 特定DEG 的识别可全面了解肿瘤的侵袭性。本研究通过比较发现，与PTC组织相比，FTC、ATC、MTC 组织中均存在独特表达的基因。基因集富集分析发现，细胞周期、细胞有丝分裂、细胞增殖、免疫信号失调与不同病理类型甲状腺癌发生发展密切相关。重要模块分析和基因功能注释分析也表明关键模块也与细胞周期、细胞增殖有关。

有研究指出，细胞周期基因在甲状腺癌中被高度激活[5-6]，本研究提供了重要模块转录因子调控网络的预测，确定EBF1、MZF1、E2F1、STAT6、NFKB1、NFKB2、RELA、CHURC1、AP3B1、ZBTB7A、IKZF2、USF2、STAT1、OVOL2基因均是细胞周期和免疫调节相关模块的潜在主调控因子。研究显示，E2F1对控制G1/S 期基因表达至关重要[7-8]。本研究分析结果显示，E2F1是甲状腺癌中重要的特异因子，与Yang 等[9]的研究结果相似。此外，本研究证实了E2F1的潜在转录靶点可能是肿瘤细胞周期失调的机制。研究表明，转录因子NFKB可协同E2F靶基因的转录并调节细胞周期[10]，NFKB的异常表达与多种肿瘤有关，其中就包括甲状腺癌[11-12]。IKZF2是造血特异性转录因子，参与调节淋巴细胞发育，可控制细胞的自我更新和分化等功能[13]，也具有肿瘤抑制因子的特性[14]。本研究同样发现，IKZF2是甲状腺癌的重要转录因子之一，但IKZF2如何在不同的细胞类型中发挥不同的作用尚不清楚。多数肿瘤患者死亡是由于肿瘤细胞转移到其他器官，而侵袭是转移形成的先决条件，有研究指出，MZF1是致癌因子，是驱动病变恶性转化的重要节点[15]。本研究结果显示，MZF1是甲状腺癌患者预后的保护因素，表明其能调控肿瘤中特定靶基因的激活，但对于MZF1与甲状腺癌全基因组相关性的研究仍待进一步深入研究。本研究筛选出的GSN、RASA1和EIF2S1等25 个关键基因的表达具有潜在的促进甲状腺癌发生发展的作用。由此可见，E2F1、NKFB1、NFKB2、IKZF2、MZF1主调控因子及其潜在靶点所形成复杂的细胞周期基因调控网络可能是治疗效果差或易复发甲状腺癌快速进展的原因。

综上所述，本研究整合转录组学研究揭示了不同病理类型甲状腺癌共同的恶性特征及相关性，E2F1、NKFB1、NFKB2、IKZF2、MZF1被认为是调控网络中共同表达的DEG 的主要调控因子，这为甲状腺癌的恶性机制提供新的见解。但本研究仍存在不足之处，未来仍需考虑一些因素如本文中的3个数据集中的样本不对称性及恶性程度高的甲状腺癌样品量少所带来的潜在影响。尽管存在局限性，但这些发现仍可支持与甲状腺癌相关的主要调控因子和网络在未来研究中的特殊价值。