基于生物信息学的胆囊癌差异表达谱中关键蛋白调控基因分析

韩伟光，莘玮，苏水霞，王青

（空军军医大学第二附属医院 普通外科，陕西 西安 710038）

胆囊癌（gallbladder cancer，GBC）作为胆道系统中最常见的恶性肿瘤，具有高度侵袭性，其病死率高且预后不良[1-2]。GBC占胆道恶性肿瘤的80%～95%，其总发生率为0.8%～1.2%[3]。以往研究[4]表明，超过85%的GBC患者被诊断时患有胆结石，证实胆结石的存在与GBC发生密切相关，是GBC的一个重要危险因素。然而，胆结石患者中少于3%会发展为GBC，这表明需要额外的因素来驱动GBC的进展。因此，解析GBC发生的潜在分子生物学机制对GBC的诊断及治疗具有重要的临床意义。

在组织学上，GBC主要分为腺癌和鳞癌两种，分别占大约90%和1%～12%[5-6]。手术切除是GBC患者可能治愈的一线选择，但由于其在早期发病阶段通常无明显症状，因此往往难以进行早期诊断[7-8]。当出现黄疸和疼痛等明显症状时，患者大多已经处于晚期，并伴有浸润和转移，普遍错过了接受根治性手术的机会[9-10]。一些非手术方法，如化学疗法、放射疗法和分子靶向疗法等可为GBC患者提供姑息性缓解，减少肿瘤生长和抑制肿瘤转移[11]。然而，在过去的几十年中，GBC的总生存期并没有得到明显的改善[12-14]。因此，对GBC潜在发病机制的分子理解以及研究显得尤为重要，同时发现特异的标记物用于GBC的早期诊断和靶向治疗也面临着巨大的挑战。

随着高通量测序技术的发展，基因组与转录组层面的数据为研究GBC疾病的分子生物学机制提供了契机。目前已有研究[15-16]表明基因组的不稳定性和变异性可能影响GBC疾病的发生发展。通过二代测序技术分析已鉴定出与GBC疾病发生有关的许多基因突变，包括TP53和ERBB3[17-18]。本研究通过转录组层面的GBC数据来筛选差异表达基因，旨在发现GBC中异常表达的基因与相关生物学过程。此外对出现的关键蛋白调控基因进行评估，为GBC患者的治疗提供新思路，有助于对GBC疾病进行更深入认识。

1 资料与方法

1.1 材料

为研究GBC患者差异基因的表达特征，本文在美国国家生物技术信息中心（National Center for Biotechnology Information，NCBI）的基因表达汇编（Gene Expression Omnibus，GEO）数据库中下载人类GBC相关的转录组公共数据，数据集编号为GSE100363和GSE139682，数据信息和样本分类见表1。

表1 GBC 相关转录组数据信息和样本分类情况Table 1 Data information of GBC-related transcriptome and sample classification

1.2 方法

1.2.1 转录组数据预处理GSE100363 数据集的预处理过程：⑴去除每段测序读段的接头序列并过滤低质量序列；⑵使用bowtie2 软件[19]将序列与参考基因组hg19 进行比对；⑶在RSEM 软件[20]中计算基因表达的FPKM（fragments per kilobase million mapped reads）值。GSE139682 数据集的预处理过程：⑴对序列进行过滤，除去碱基质量低于13 的20%的读段；⑵使用HISAT2.1.0 软件[21]将序列比对到参考基因组hg38；⑶根据Chepelev等[22]报道的方案计算基因表达的RPKM（reads per kilobase per million mapped reads）值。

1.2.2 筛选转录组数据差异表达基因在GBC 和正常样本中，同时使用双总体Student t 检验和差异倍数FC（fold change）法来识别差异表达基因，本文定义符合阈值为P＜0.05 且|log2FC| ≥1的基因为差异表达基因（differentially expressed genes，DEGs）。

1.2.3 差异表达基因的GO 注释为深入了解GBC 中差异表达基因的相关生物学过程，本文使用R软件中clusterProfiler、org.Hs.eg.db、DOSE 等程序包对差异表达基因进行GO（gene ontology）功能注释[23]，包括基因参与的生物学过程（biological process，BP）、基因所处的细胞组分（cellular component，CC）和基因执行的分子功能（molecular function，MF）。在此采用多重检验校正的方法Benjamini-Hochberg 对P 值进行校正，设置满足校正P＜0.05 的GO 术语具有显著性。

1.2.4 差异表达基因的蛋白质互作网络构建与模块挖掘随后采用STRING[24]蛋白质互作数据库中的蛋白质互作信息对GBC 差异表达基因构建蛋白质互作（protein-protein interaction，PPI）网络，并用Cytoscape 软件对蛋白质互作网络进行可视化，设置最低互作分数为0.4。同时使用Cytoscape 中的图论聚类算法MCODE（molecular complex detection）在构建的蛋白互作网络中进行模块挖掘，设置度的阈值为2，节点阈值为0.2。

1.2.5 关键蛋白调控基因的表达分析和预测能力评估采用Graphpad Prism 8.3.0 软件对最终保留的关键蛋白调控基因表达量进行两独立样本t 检验，P＜0.05 为差异具有统计学意义。并且对关键蛋白调控基因绘制受试者工作特征曲线（receiver operating characteristic curve,ROC）。ROC 曲线下面积（area under curve，AUC）越接近1.0，其诊断效果越好。

2 结果

2.1 GBC 差异表达基因分析结果

在两个转录组数据中，同时采用t检验和FC的差异基因分析方法，设置P＜0.05且|log2FC|≥1的阈值条件。GSE100363数据集筛选出644个（2.54%）差异表达基因，其中上调基因282 个（43.79%），下调基因362 个（56.21%）。GSE139682数据集筛选出2206个（5.28%）差异表达基因，上调基因829个（37.58%），下调基因1377个（62.42%）。两个数据集中可重复差异基因144 个（图1 A）。共同上调基因20 个（图1B），共同下调基因120个（图1C），表2列举了在各自数据集中差异表达最大的5个共同上调基因和5个共同下调基因。由于这140个共同差异表达基因可能在GBC发生发展中发挥一定作用，故随后的分析都基于此。

图1 GBC 表达基因的Venn 图 A：两个数据集的共同差异表达基因Venn 图；B：两个数据集共同差异表达的上调基因Venn 图；C：两个数据集共同差异表达的下调基因Venn 图Figure1 The Venn diagram of differentially expressed genes in GBC A:The Venn diagram of differentially expressed gene in GBC in both data sets;B:The Venn diagram of up-regulated genes in both data sets;C:The Venn diagram of down-regulated genes in both data sets

表2 GBC 组织样本在各自数据集中差异表达最大的5 个共同上调和5 个共同下调基因Table 2 The 5 common up-regulated and 5 common down-regulated genes with the largest differential expression in GBC tissue samples in each dataset

2.2 差异表达基因的功能富集

为研究差异表达的140个基因的相关生物学功能，对这些基因进行GO注释，分别富集到2249条生物学过程术语、247 条细胞组分术语和295 条分子功能术语。其中符合阈值校正P＜0.05 的生物学过程术语454 条，从最显著的5 条，可以看出这140个表达异常的基因主要影响前脑发育、端脑发育和神经发生的正调控等生物学过程，体现了GBC相关的通路出现了不同程度干扰。在这5条通路中，CDON（cell adhesion associated,oncogene regulatged）、FEZ1（fasciculation and elongation protein zata 1）、GLI3（GLI family zinc finger 3）、SLIT2（slit homolog 2）、SOCS2（suppressor of cytokine signaling 2）基因出现频率最高，表明这5个基因主要参与了GBC的生物学发生过程（图2）。符合阈值的2条细胞组分术语，其中9个基因参与突触后膜组成，4个基因参与横小管组成。另外在阈值的设定下没有符合的分子功能术语（图3）。

2.3 差异表达基因PPI 网络分析及模块挖掘

对这140个差异表达基因在STRING数据库中筛选互作值≥0.4的互作对，之后用Cytoscape软件对蛋白互作网络进行可视化（图4）；用Cytoscape软件中的MCODE插件在PPI网络中进行模块挖掘（表3），PPI网络图中不同的颜色对应表3中不同的模块。结果表明SFRP1基因与GLI3、CCND2互作形成模块1；SEMA3D基因与SLIT2、NTN1互作形成模块2；ANK2基因与GNAO1、PTPRD互作形成模块3。在不同基因表达的蛋白互作下，3个种子基因对于GBC的生物学过程可能显得更为重要。

图2 差异表达基因的GO 生物学过程注释（纵坐标表示富集到的不同GO 生物学过程，横坐标表示注释到某条术语的基因数目；不同色阶表示为不同校正P 值，颜色越红表示校正P 值越小，其差异越显著，越蓝表示校正P 值越大越不显著）Figure 2 Annotation of GO biological process of differentially expressed genes (the vertical axis showing enrichment of different biological processes,the horizontal axis showing the number of genes annotated to a certain term;different color gradations standing for different adjusted P values,namely the darker the red color,the smaller the adjusted P value and the more significant the difference,and the darker the blue color,the bigger the adjusted P value and the less significant the difference )

图3 差异表达基因的GO 细胞组分注释（纵坐标表示富集到的不同GO 细胞组分，横坐标表示注释到某条术语的基因数目；图中展示的两条术语校正P 值在0.036附近）Figure 3 Annotation of GO cellular component of differentially expressed genes (the vertical axis showing enrichment of different cellular components,the horizontal axis showing the number of genes annotated to a certain term;the adjusted P values of the two terms displayed in above picture are around 0.036)

图4 差异表达基因的PPI 网络图Figure 4 PPI network diagram of differentially expressed genes

表3 PPI 网络中对应的模块信息Table 3 The corresponding module information in PPI network

2.4 关键蛋白调控基因表达分析情况

将GBC组织和正常胆囊组织相比，采用两独立样本t检验，发现GBC组织的关键蛋白调控基因SFRP1、SEMA3D、ANK2的表达均明显下降，其差异均具有统计学意义（均P＜0.05）（图5）。这一结果表明调控这3种基因的高表达可能降低GBC的发生。SFRP1和SEMA3D基因的差异情况略优于ANK2基因。

图5 关键蛋白调控基因表达情况分析Figure 5 Analysis of expressions of the key protein regulating genes

2.5 评估关键蛋白调控基因对GBC 的预测能力

通过绘制ROC曲线来评估关键蛋白调控基因诊断和预测GBC的效果。结果发现SFRP1和ANK2基因对GBC具有一定的预测能力，其AUC值均＞0.9。SFRP1基因的预测效果略优。SFRP1、ANK23、SEMA3D联合应用的AUC值最大，为0.9643（图6）。

图6 关键蛋白调控基因诊断能力分析Figure 6 Analysis of diagnostic ability of the key protein regulatory genes

3 讨论

GBC是最常见的胆道恶性肿瘤，由于缺乏特异性症状和可靠敏感的疾病标志，大多数GBC患者被诊断即为晚期，因此迫切面临对分子机制的确切理解与研究。本研究旨在利用转录组层面的高通量数据，使用T检验和倍数变换法筛选差异表达基因，结果发现在两个数据集中符合阈值为P＜0.05且|log2FC|≥1的共同出现的差异基因140个（上调20个，下调120个）。通过GO中3种类型的功能注释，在阈值为校正P＜0.05的条件下，发现差异基因主要与前脑发育、端脑发育和神经发生的正调控等生物学过程有关，与突触后膜和横小管这两种细胞组分密切相关，而没有发现与分子功能显著相关的注释信息。随即对差异表达基因构建蛋白互作网络以发现存在明显相互作用的基因。值得注意的是3个种子节点的基因SFRP1、SEMA3D、ANK2，它们与正常胆囊组织相比，在GBC组织中表达量均为明显下降趋势，可能在GBC发生中起着尤为重要的作用。其中，SFRP1基因相比于其他两个基因的差异表达量最显著且诊断效果也最明显。

既往研究证实由SFRP1基因所编码的SFRP1蛋白作为Wnt信号蛋白的拮抗剂，在Wnt/catenin信号通路的调控中发挥重要作用。此外，SFRP1基因在胆道系统中也被报道是高度相关的甲基化基因，基于其甲基化的分子标记物已成功地用作常规、预后及预测指标，以更好地管理各种癌症患者，一些研究还描述了SFRP1作为癌症的预后和预测生物标志物的价值[25-26]。在人类中，高度相关甲基化基因的生物标记，允许从相应的正常组织中区分出肿瘤组织，这些标志物还与独特的癌症表型有关，通过进一步优化分析SFRP1基因可能会增加对GBC早期的发现与诊断，并在人体液体中表现相对可靠的监测潜力[27-28]。

同时，针对目前发现的在差异最显著的5 条GO生物学过程中出现的频率最多的基因也有研究证实它们的发生与GBC有关。有研究表明，GBC的许多抑癌基因中已经描述了杂合子丢失，包括染色体1p34-36（p73），5q21（APC），8p21-23（PRLTS和FEZ1）和17p13（p53基因）等[29]。最新报道[30]发现，参与GBC发生发展的长链非编码RNA（long-chain non-coding RNA，lncRNA）SNHG6的沉默和miR-26b-5p的上调可促进细胞凋亡、抑制细胞生长和上皮间质转化，并促进GLI3和E-钙黏着蛋白表达的上调，进而影响GBC的生物学过程。另外转录因子GLI3是Hedgehog（Hh/HH）信号通路的成员，可以全长（Gli3-FL/GLI3-FL）或阻遏物（Gli3-R/GLI3-R）的形式存在。GLI3调节各种对于癌细胞生长和发展至关重要的生物过程，锚定非依赖性生长、增殖和迁移的肿瘤细胞系[31]。另外SLIT2基因的甲基化在多种肿瘤中也有不同程度的研究[32-33]。

此外，本文也存在如下不足之处。首先，由于是基于生物信息学方法进行的研究，所以部分结果还需要进一步的功能实验进行验证，因而可以为之后的临床应用提供更可靠的依据。 其次，由于GBC疾病本身发生机制的异质性和研究样本数量的有限性，快速发展的高通量单细胞层面测序技术可能会解决这一问题。最后，因为GBC疾病的发生是受多方面因素影响的，只从转录组层面看存在一定的片面性，本文尚未结合基因组、蛋白质组、代谢组等多组学数据信息，多组学数据之间的相关性有待研究。同时环境的因素也可能是GBC发生与进展的重要原因。因此，本研究结果可能适用于部分GBC患者，并非适用于所有GBC患者。

综上，本研究通过分析两个转录组层面数据集取交集的相关差异基因的GO生物学过程和细胞组分，关键蛋白调控基因SFRP1可能在GBC中发挥不可忽视的作用。研究结果具有一定的可靠性，可以为更深刻的认识GBC分子机制提供新的见解与方向。