视网膜母细胞瘤发生发展的关键基因挖掘及验证

黄子珊，符馨予，周希瑗

重庆医科大学附属第二医院眼科/眼科学重庆市市级重点实验室，重庆 400010

视网膜母细胞瘤(retinoblastoma，R B)是婴幼儿最常见的眼内恶性肿瘤[1]，年发病率为1/18 000～1/16 000，可分为遗传型(约45%)和非遗传型(约55%)[2]，如果不及时治疗可能导致失明或死亡。为实现RB的早诊断、早治疗及降低病死率，寻找RB的关键基因尤为重要。加权基因共表达网络(weighted gene co-expression network analysis，WGCNA)分析目前已广泛应用于与肿瘤相关的高通量数据挖掘中[3-4]，该分析方法可寻找协同表达的基因模块，探索基因模块与关注的表型(如临床数据)之间的关系，构建共表达网络并寻找网络中的核心基因，主要通过基因之间相关系数计算、基因模块的确定、共表达网络、模块与性状关联4个步骤来寻找影响疾病的关键基因。本研究运用多种生物信息学分析方法寻找与RB发生发展相关的关键基因，并通过实时荧光定量PCR(quantitative real-time polymerase chain reaction，qRT-PCR)及免疫组化等多种方法对筛选出来的关键基因进行验证，以进一步寻找能够用于早期诊断RB的生物标志物。

1 材料与方法

1.1 数据下载及处理 从GEO(https://www.ncbi.nlm.nih.gov/gds)基因表达数据库中下载基因表达谱数据集GSE110811、GSE24673、GSE97508，3个芯片数据集共包含43例RB样本和7例正常对照组织样本。其中，GSE110811基于GPL16686平台，包含28例RB和3例正常视网膜组织样本(本研究排除了1例正常对照样本GSM3017153，原因为该数据集样本可能混有肿瘤组织[5])；GSE24673基于GPL6244平台，包含9例RB和2例健康成人视网膜组织样本；GSE97508基于GPL15207平台，包含6例RB和3例正常视网膜组织样本。

根据平台注释信息将阵列探针名转换为匹配的基因名。采用perl软件将这3个数据集合并，用R软件的sva包消除3个数据集之间的批次效应(样品在不同时间、不同分组或由不同的人处理造成的差异)[6]。将3个数据集的矩阵数据合并为1个数据集，共包含43例RB和7例正常对照组织样本。

1.2 筛选差异表达基因(differentially expressed gene，DEGs) 使用R语言的limma包分析合并数据集中RB与正常对照样本之间的DEGs矩阵，筛选条件设置为|logFC|＞1，P＜0.05[7]。

1.3 GO(gene ontology)富集分析和KEGG(Kyoto Encyclopedia of Genes and Genomes)通路分析 对DEGs进一步分析，利用R语言clusterProfiler[8]包进行GO富集分析和KEGG通路分析。校正后的P值(FDR)＜0.05为差异有统计学意义。

1.4 蛋白相互作用(protein-protein interaction，PPI)网络 将DEGs输入STRING数据库，以交互得分=0.9(最高置信度)作为阈值构建PPI网络，然后导入Cytoscape(3.7.2版本)软件中进行可视化分析和整理，运用cytoHubba插件以最大集团拓扑分析法(maximal clique centrality，MCC)[9]选取关键节点的基因[10]。

1.5 WGCNA分析 提取GSE110811肿瘤组28个样本对应的DEGs矩阵进行WGCNA分析，获取该肿瘤样本对应的临床信息。

1.5.1 共表达网络的构建与模块划分 利用R软件的WGCNA包[11]，计算各基因间的Pearson相关系数，选择适当的软阈值β构建无尺度网络；建立邻接矩阵，转换为TOM重叠矩阵；表达模式相似的基因分为一类，并构成模块；基于TOM矩阵，使用平均连锁层次聚类方法来聚类基因，根据混合动态切割树的标准设置最小模块，构建加权共表达网络以筛选hub基因[12]，以0.5作为加权共表达相关系数的阈值，导入Cytoscape(3.7.2版本)中进行可视化分析，找到位于网络中心的核心基因。

1.5.2 关键模块筛选 将样本的临床信息与模块基因进行关联分析，寻找显著相关的模块。本研究患者临床信息包括：有无远处转移、有无继发性肿瘤、RB1基因突变、单双侧眼患病、诊断年龄和随访时间等。采用Pearson相关系数评估模块与临床信息之间的相关性，相关系数越大表示该模块与临床数据相关程度越高，该模块即为关键模块，并计算各个基因与临床信息的相关系数。

1.5.3 关键模块基因GO和KEGG分析 使用DAVID数据库对关键模块的基因行GO富集分析和KEGG通路分析(P＜0.05为显著性基因富集的筛选条件)。

1.6 实时荧光定量PCR(qRT-PCR)检测关键基因的表达量

1.6.1 细胞培养及分组 人RB细胞株(Y79、WERIRB-1)分别购自上海佰晔生物科技中心和上海富衡生物科技公司；人视网膜色素上皮细胞株(ARPE-19)购自上海富衡公司。Y79及WERI-RB-1细胞培养基均利用RPMI 1640培养基(美国Hyclone公司)加20%胎牛血清(上海依科赛生物)及1.0%青链霉素溶液(上海碧云天公司)配制，ARPE-19细胞培养基为含10%胎牛血清(上海依科赛生物)及1.0%青链霉素溶液的DMEM-F12(美国Hyclone公司)培养基；所有细胞均在37 ℃、5%CO2培养箱中培养。Y79、WERI-RB-1细胞作为实验组，ARPE-19细胞作为对照组。

1.6.2 qRT-PCR实验步骤 收集细胞，PBS洗2或3次，加入Trizol裂解液(cat：15596026，美国Invitrogen公司)，依次使用氯仿、异丙醇和75%乙醇提取总RNA；利用超微量分光光度计对RNA样品进行浓度和纯度测定，利用反转录试剂盒(日本TaKaRa公司)并按照操作说明书配制反应液，分两步进行实验将RNA反转录为cDNA；使用SYBR Green qPCR Master Mix试剂盒(MCE公司)配制反应体系，依次加入无酶水、SYBR Green染料、上下游引物、cDNA模板后混匀，用ABI 7500荧光定量PCR仪进行扩增反应。设置热循环条件为：95 ℃ 5 min；95 ℃ 10 s、60 ℃ 30 s、95 ℃ 15 s、60 ℃1 min、95 ℃ 30 s、60 ℃ 15 s，共40个循环[13]。以β-actin作为内参照，5个关键基因引物序列见表1(由上海生工有限公司合成)。实验重复3次。采用2-ΔΔCt法计算目的基因的相对表达量。

表1 qRT-PCR引物序列Tab.1 Primer sequences for qRT-PCR

1.7 免疫组化检测关键蛋白表达量

1.7.1 临床标本收集及分组 纳入2016年1月－2019年1月重庆医科大学附属第二医院收治的RB患者石蜡组织标本20例，所有患者术前均未进行放化疗，获取10例因车祸伤等原因进行眼球摘除患者的正常视网膜作为对照组。本研究经重庆医科大学附属第二医院伦理委员会审批[(2019)270]。

1.7.2 免疫组化实验步骤 石蜡切片常规脱蜡后，将组织切片置于乙二胺四乙酸(EDTA)抗原修复缓冲液(pH9.0)中进行抗原修复，放入3% H2O2溶液避光孵育25 min；滴加3%牛血清白蛋白(BSA)室温封闭30 min；分别加入着丝粒蛋白K(CENPK)抗体(北京博奥森公司)、兔抗人胞质分裂蛋白调控因子-1(PRC1)抗体(Abcam公司)4 ℃孵育过夜[14]，加入辣根过氧化酶(HRP)标记的山羊抗兔二抗；用3,3-二氨基联苯胺(DAB)显色、苏木精复染，脱水后中性树脂封片，组织切片在显微镜下进行检查。免疫组化的结果分级标准(0～12分)：0～4分为低表达组(阴性)，5～12分为高表达组(阳性)[15]。上述试剂除一抗外均购买自武汉赛维尔生物公司。

1.8 统计学处理 采用GraphPad Prism 8.3.0软件进行qRT-PCR结果作图及统计分析，多组间关键基因表达量的比较采用单因素方差分析，进一步两两比较采用Dunnett-t检验；采用SPSS 25.0软件对免疫组化结果进行分析，RB组织与对照组中关键蛋白的表达差异采用单侧Fisher确切概率法进行比较。P＜0.05为差异有统计学意义。

2 结 果

2.1 DEGs的筛选 共得到1254个DEGs，其中表达上调基因422个，表达下调基因832个，具体分布情况见图1。

图1 视网膜母细胞瘤组织与正常组织差异表达基因的火山图及热图Fig.1 Volcano and heat maps of differentially expressed genes in RB tissues and normal tissues

2.2 DEGs的功能注释和通路富集 GO富集分析结果显示，DEGs主要与蛋白质异二聚体活性、阳离子跨膜转运蛋白活性、染色质结合、无机阳离子跨膜转运活性和细胞黏附分子结合等生物学功能密切相关(图2A)。KEGG通路分析结果显示，DEGs主要富集在细胞周期、光传导、DNA复制和卵母细胞减数分裂等通路上(图2B)。

图2 视网膜母细胞瘤差异表达基因的GO富集分析(A)和KEGG通路分析(B)Fig.2 GO enrichment (A) and KEGG pathway analysis (B) of differentially expressed genes in RB

2.3 PPI网络构建及核心基因筛选 PPI网络由3567条边、611个节点构成，以MCC＞2.09×1013为标准得到79个基因[9]，即为网络核心基因(图3)。

图3 视网膜母细胞瘤差异表达基因PPI网络构建及核心基因筛选Fig.3 PPI networks of DEGs and screening core genes in RB

2.4 WGCNA分析结果

2.4.1 共表达网络构建与枢纽基因筛选 利用pickSoftThreshold函数筛选合适的软阈值β，使无尺度拓扑拟合指数R2＞0.8以构建共表达网络(图4A)，β=6时共表达网络中基因之间的连接性满足无标度网络分布(图4B)，利用MCC算法以得分＞4作为标准[9]，得到位于网络中心的前13个基因，即为枢纽基因(图4C)。

图4 视网膜母细胞瘤WGCNA网络构建及其枢纽基因的筛选Fig.4 Screening of WGCNA network and key genes in RB

2.4.2 基因模块的划分与关键模块筛选 利用动态混合剪切法获得11个共表达基因模块，其中灰色模块是无法聚集到其他模块的基因，基于模块基因构建模块聚类树(图5A)，模块的拓扑重叠热图反映模块与模块间基因的相关程度(图5B)，基因的聚类树状图展示模块划分过程(图5C)。分别计算各个模块与临床性状之间的Pearson相关系数及P值(图6A)，以Pearson相关系数的绝对值＞0.55且P＜0.05作为筛选条件，筛选出5个重要模块，分别是blue、pink、turquoise、red和brown模块。Blue、pink和turquoise模块与诊断年龄呈明显正相关(r=0.65，P=2×10-4；r=0.58，P=0.001；r=0.59，P=0.001)，其中pink模块还与单双侧眼患病和随访时间呈明显负相关(r=-0.57，P=0.002；r=-0.63，P=4×10-4)；而red和brown模块与继发性肿瘤呈明显正相关(r=0.57，P=0.001；r=0.55，P=0.002)。其中blue模块与诊断年龄的相关系数最高，将blue模块包含的基因导入STRING数据库构建PPI网络(图6B)。各个重要模块包含的基因见表2。本研究还将基因与临床性状关联，分别计算出每个基因与临床性状的相关性，选取与每项临床性状相关系数绝对值最高的基因，得到hmgb3、pcdhb3、ids、polh、birc5、guca1b共6个基因(表3)，其中pcdhb3与诊断年龄的相关系数最高(r=0.8411)。

图5 视网膜母细胞瘤WGCNA分析模块划分过程Fig.5 WGCNA analysis module partition process in RB

表2 视网膜母细胞瘤研究中各个重要模块包含的基因Tab.2 Genes contained in each important module of research on RB

表3 视网膜母细胞瘤与临床性状相关程度较高的基因列表Tab.3 The genes with the higher degree of correlation for each clinical trait in RB

2.4.3 GO富集分析和KEGG通路分析 Pink和red模块包含基因数目较少，无法进一步进行富集分析，因此将blue、turquoise和brown模块中的基因导入DAVID数据库中，以校正后的P值(FDR)＜0.05作为筛选条件。Blue模块GO分析结果显示，对于生物过程主要富集在细胞分裂、有丝核分裂等过程，对于细胞组分主要富集在细胞核、细胞质、微管等结构，对于分子功能主要富集在ATP结合、微管结合等功能方面；KEGG分析结果显示主要富集在细胞周期、卵母细胞减数分裂等通路。Turquoise模块基因主要与细胞核、核质、ATP结合等相关，brown模块基因主要富集在细胞核、核质等细胞成分上。

2.5 关键基因筛选 将DEGs PPI的核心基因、WGCNA网络的枢纽基因及重要模块包含基因取交集(图6C)，得到染色体结构维持蛋白4(structural maintenance of chromosome 4，smc4)、微小染色体维持蛋白6(minichromosome maintenance complex component 6，mcm6)、着丝粒蛋白K(centromere protein K，cenpk)、驱动蛋白家族成员15(kinesin family member 15，kif15)、胞质分裂蛋白调控因子-1(protein regulator of cytokinesis 1，prc1)等5个关键基因。上述5个关键基因在RB中均为高表达基因，其logFC均＞1，关键基因所属模块及差异表达信息见表4。

图6 模块与临床性状相关性热图及blue模块PPI网络Fig.6 Heat map of correlation between module and clinical traits and PPI network of blue module

表4 视网膜母细胞瘤关键基因所属模块及差异表达信息Tab.4 The modules and differential expression information of key genes in RB

2.6 关键基因在细胞株中的表达情况 采用qRTPCR检测5个关键基因的相对表达量，结果显示，prc1和cenpk的相对表达量在Y79、WERI-RB-1细胞中明显高于对照组(P＜0.05，图7)。

图7 视网膜母细胞瘤关键基因的相对表达量Fig.7 Relative expression levels of key genes in RB

2.7 RB及正常视网膜组织中PRC1、CENPK蛋白的表达情况 免疫组化结果表明，20例RB中有10例(50%)PRC1蛋白表达阳性，而10例正常视网膜中仅有1例(10%)呈阳性，故RB组织中PRC1蛋白阳性率高于正常视网膜，且差异有统计学意义(P＜0.05)。CENPK蛋白在RB和正常视网膜中的表达差异无统计学意义(P＞0.05，图8)。

图8 PRC1及CENPK蛋白在视网膜母细胞瘤组织及正常视网膜中的表达(免疫组化染色, ×200)Fig.8 Expression of PRC1 and CENPK protein in RB tissue and normal retina (Immunohistochemical staining, ×200)

3 讨 论

RB发生率居小儿恶性肿瘤的第2位[1]，严重威胁患儿生命。RB发病多在5岁以下，约2/3的患儿在3岁前患病，近年来其发病率呈上升趋势，且如果治疗不及时，可能导致失明或死亡。因此，为实现RB的早诊断、早治疗，寻找其关键基因及生物标志物显得越来越重要。

目前，高通量数据资源越来越丰富，许多学者运用多种生物信息学方法进行大样本数据挖掘[16]，寻找导致某种疾病发生的关键基因，以进一步发现与疾病诊断及预后相关的生物标志物，此方式较传统研究方法更省时省力、节约经费且科学有效。然而，以往的生物信息学分析多使用单个数据集[17]和单一方法分析DEGs[18]，而RB发病率较低且组织样本获取难度大，故组织样本较为稀缺，单个数据集样本量较少，可能导致分析结果不准确。因此，本研究将3个RB数据集中的原始数据进行合并以增加样本量，运用差异表达分析联合加权基因共表达网络分析，将DEGs划分成不同模块并与临床性状相关联，筛选出了与RB发生发展密切相关的5个可能基因。其中，SMC4对哺乳动物的DNA修复至关重要，可能参与了基因表达沉默状态的维持和DNA修复等过程，与肺腺癌、神经胶质瘤及结直肠肿瘤发生相关，且与乳腺癌预后相关[19-22]。MCM6是高度保守的微型染色体维持蛋白之一，是启动真核基因组复制所必需的，Liu等[23]发现其可促进肝细胞癌转移，且可作为肝细胞癌早期复发的一种新的血清生物标志物[24]。CENPK与着丝粒功能和有丝分裂进程密切相关，有研究发现，它可能参与了肝细胞癌的恶性化进展[25]。Lee等[26]发现，CENPK在卵巢癌细胞中表达上调且可作为卵巢癌的一种新的肿瘤标志物。KIF15可促进黑色素瘤生成[27]，也是促进胰腺癌增殖的重要调控因子[28]。PRC1基因位于15q26.1，与细胞增殖、胞质分裂、纺锤体及微管的形成密切相关。

本研究的生物信息学分析结果显示，这5个基因均为RB组织中表达上调的基因。WGCNA分析结果显示它们均为重要模块中的基因，其中smc4属于blue模块，mcm6、prc1属于brown模块，cenpk、kif15则是turquoise模块的组成部分。由于blue模块与诊断年龄呈明显正相关，brown及turquoise模块与继发性肿瘤呈明显正相关，因此，smc4可能与RB的发病进程相关，而mcm6、prc1、cenpk、kif15可能与RB继发性肿瘤的发生存在一定相关性。但由于目前这些基因在RB中的研究较少，其在RB中的功能仍需进一步探索。另外，通过对上述关键模块基因的GO及KEGG分析发现，此类基因主要参与细胞核、微管等结构形成，与细胞周期及有丝分裂过程等通路密切相关，因此，推测这些关键模块中的差异表达基因可能是通过调控细胞生物学过程及信号通路来影响RB的发病进程的。此外，本研究还将基因与临床性状关联，分别计算出每个基因与临床性状的相关性，选取与每项临床性状相关系数绝对值最高的基因，得到了pcdhb3等共6个基因，其中pcdhb3与临床性状(诊断年龄)的相关系数最高。有研究发现，PCDHB3是一种新型的结直肠癌肿瘤抑制因子，且可能是晚期结直肠癌的预后标志物[29]，因此，此类基因在RB中的功能也值得进一步研究。

为验证上述筛选出的关键基因是否在RB细胞株中高表达，本研究体外培养了Y79、WERI-RB-1及ARPE-19细胞进行qRT-PCR实验，结果表明prc1及cenpk这两个基因在RB细胞中高表达(P＜0.05)。此外，通过前面的生物信息学分析可以发现，在这5个基因中cenpk、prc1在RB中的表达明显上调，而且它们在肿瘤的发生中发挥重要作用。结合qRT-PCR实验结果，我们挑选cenpk、prc1这两个基因进行进一步实验，即在RB组织样本中利用免疫组化实验在蛋白质水平进行验证，结果表明PRC1蛋白在RB组织中的表达量明显增高(P＜0.05)。PRC1在恶性肿瘤中的作用研究尚少，Wang等[30]发现，PRC1高表达在肝癌细胞中增加了癌细胞的化疗耐药性。Liao等[31]对2株RB细胞(HXO-RB44和WERI-Rb-1)进行siRNA转染使prc1表达量降低后，发现RB细胞的增殖和血管生成受到了抑制，该过程是通过抑制Wnt/β-catenin通路实现的，表明PRC1蛋白对RB的发生具有重要影响。

综上所述，本研究运用差异表达分析联合加权基因共表达网络分析提取出了与RB相关的关键基因，其中RB细胞及组织中的PRC1表达量明显高于对照组，且Liao等[31]的体外实验表明沉默prc1基因可抑制RB细胞的增殖和血管生成，这些结果均提示PRC1可能对RB的发生发展具有一定的促进作用。此外，本研究对探究RB的发病机制及潜在的治疗靶点提供了新视角，但还需进行进一步体内实验来验证这些关键基因的生物学功能。