用基因互作网络识别胰腺导管腺癌基因生物标记

王新格，田卫东

(复旦大学 生命科学学院 生物统计学与计算生物学系，上海 200438)



用基因互作网络识别胰腺导管腺癌基因生物标记

王新格，田卫东

(复旦大学 生命科学学院 生物统计学与计算生物学系，上海 200438)

胰腺导管腺癌(PDAC)是全球高致死率癌症中的一种.PDAC基因生物标记的识别可以通过构建基因互作网络完成.利用蛋白质互作网络来分析与研究基因表达芯片数据，构建出PDAC基因互作网络并对其划分基因模块，进而筛选出在模块中的PDAC差异性表达基因.通过筛选在癌症样本和正常组织样本中共表达的基因对，并利用STRING蛋白质互作网络评估基因功能相关性，构建出具有PDAC特异性的基因互作网络.利用iNP算法进行网络模块化分，在每一个模块中，模块内基因都具有强的共表达特性和模块功能相关性.通过筛选，获得了34个基因模块，其中20个在癌症样本中表达明显上调，14个在癌症样本中表达明显下调.从这些模块中又筛选出在PDAC样本中表达上调的10个基因生物标记，如DMBT1、DSC3等和表达下调的10个基因生物标记，如DLG5、NRCAM等.

胰腺导管腺癌； 基因生物标记； 基因互作网络； 功能富集分析

胰腺导管腺癌(Pancreatic Ductal Adenocarcinoma， PDAC)是一种致死率很高的癌症.作为消化道常见的恶性肿瘤，其多发生于胰腺头部.其常见的致病原因为饮酒和糖尿病.由于PDAC在早期没有明显的症状，因此在发现时一般患者已处于病程的晚期.其病情发展快且恶化迅速，最常见为上腹部饱胀不适、疼痛.患者平均存活时间为6个月，其中超过5年存活的患者只占6%[1].2011年美国的统计情况为新增病例约4万，同时致死人数约3万人.由于PDAC的恶性程度高，诊断和治疗都十分困难.只有20%的病人可以通过手术切除来获得治愈的可能.手术切除后，80%的患者存活时间少于2年[2-3].然而，约12%的术后患者可以存活5年，这同所患癌症所在阶段和手术效果相关.一线药物吉西他滨是治疗胰腺癌的典型药物，但其疗效只能达到12%[4].在最新的药物实验当中，像奥沙利铂、依立替康、氟二氧嘧啶和甲酰四氢叶酸等药物能够显著延长受试者的存活时间[5].目前PDAC的活体检测可以通过穿刺抽吸胰腺肿块做细胞学检查进行.但此诊断常用于晚期病人，对于早期诊断并不便捷有效.鉴于PDAC死亡率和发病率近年来明显上升，对它的致病机理的研究以及找到新的有效的药物靶点都是临床医学迫切需要解决的问题.

一种识别和发现PDAC潜在药物靶点的方法就是通过分析PDAC基因芯片表达数据来筛选致病胰腺癌基因[6-8].一些研究已经筛选出可能同PDAC相关的致病基因，并对其致病机理做了进一步分析和研究[9-11].本研究通过对PDAC的基因芯片表达数据中共表达基因对的筛选，结合STRING蛋白质互作网络，构建出胰腺癌特异性基因互作网络.通过iNP[12]算法对网络进行模块化分，找出在PDAC中表达显著差别的基因网络.最终筛选出20个模块中的基因生物标记，这些标记可以作为预测和区分PDAC样本和正常样本的参考指标.

1 数据和方法

1.1 数据收集和处理

本研究通过GEO query[13]工具，从GEO数据库[14]中下载PDAC相关的基因芯片表达数据.选取Affymetrix Human Gene 1.0 ST Array芯片为分析平台，芯片代码为GDS4336[15].该芯片含有90个样本，其中45个为PDAC样本，另外45个为对应的邻近的正常组织样本.通过Bioconductor中limma包[16]对芯片表达值进行正态化后，我们可以整理得到2个数据集，分别为PDAC基因表达数据集和其对应邻近正常组织基因表达数据集.数据概况详见表1.对于每一个数据集，计算其中每一对基因表达值之间的皮尔森相关系数(Pearson Correlation Coefficient， PCC)并排序，从而选取排名前1%和后1%的基因对进行后续分析.依据PCC排序结果筛选出在2个数据集中具有明显相关性的基因对.例如，按照PCC从大到小的顺序，前1%的基因对为在该数据集中具有明显的正相关表达关系.相反地，后1%的基因对则表明在该数据集中具有明显的负相关表达关系.筛选后，将从2个数据集中筛选出的基因对合并，并对应到STRING在线数据库[17]中进行功能相关性评估.将功能相关性打分高于500的基因对筛选出后，构建了一个关于PDAC的特异性基因互作网络.接着，运用iNP基因互作网络模块划分算法将此网络划分为更为清晰的若干基因模块，以利于后续进一步分析.

表1 PDAC样本和正常样本基因芯片数据分布

1.2 识别PDAC特异性基因模块

经过iNP算法将整个关于PDAC的特异性基因网络划分为若干小的基因模块后，将分析各个模块中的基因是否在癌症和正常样本中表达值有明显区别.由于每个模块是由一组功能显著相关的基因所组成，我们可以将一个模块看成是一条表达通路.证明一条通路在癌症和正常样本中是否有显著区别的方法有很多，本研究采用一个简单的方法，中值对比检验法.每个模块中的每一个基因都可以找到其在45个PDAC样本和45个正常组织中的表达值，因此每个基因都有一对分别来自癌症样本和正常样本的基因表达中值.运用t检验进行统计显著性分析，由于每个模块包含基因众多，本研究通过对伪发现率(False Discovery Rate， FDR)的控制，定义P<0.01才构成统计显著性.除此条件外，log10(表达值)>2才构成有意义的选取.经过以上步骤，我们可以从一个关于PDAC的特异性基因大网络中精选出具有统计显著性的模块集合.获得模块后，通过绘制ROC曲线，对每个模块在PDAC样本中表达值是上调还是下调进行深度分析.利用ROC曲线下的面积AUC计算此调控方向.AUC数值分布从0到1，其中0和1代表该模块在癌症样本中绝对的上调或者下调，而0.5则代表模块中的基因表达为随机.本研究中计算出的AUC数值如果大于0.5，则代表该模块在癌症中表达上调，相反地，则代表模块在癌症样本中表达下调.通过Cytoscape 3.3.0软件绘制模块的结构.

1.3 功能富集分析

给定一基因模块，对其中所包含基因进行功能富集分析.先下载GO基因功能注释文件[18]和GSEA[19]数据库中的MSigDB分子标签数据库.本研究利用Fisher算法检测GO注释和MSigDB分子标签的富集情况.统计显著性初始阈值选为0.01.由于功能富集出来的基因集重合过多，本研究通过聚类筛选的方法筛选去除冗余数据.聚类时通过计算所有基因集合的关联度(重合基因数/所有基因数)来形成新的基因网络.进而，利用iNP算法划分网络为若干基因模块，其中基因集之间仍有高度重合.我们通过整合模块特异性和基因表达值差异给模块中基因打分.一个基因共表达差异值等于该基因在该模块中所有互作共表达差异值的均值.一个互作共表达差异值等于两个互作基因在癌症样本和正常样本中的差异表达值.在筛选步骤中，将所有功能富集分析出的基因集对应到划分出的基因模块中，同时筛选出模块内部基因富集统计显著性最强的基因集作为该模块代表基因集.可以取排名前10位的基因作为生物标记的候选对象.最终，本研究收集出所有具有功能富集统计显著性的基因模块进行研究.

2 结果和分析

2.1 构建PDAC功能相关基因网络

本研究整理获得两个基因表达数据集，分别对应正常胰腺组织(样本个数为45)和PDAC组织(样本个数为45).基因表达数据下载于GEO DataSets数据库(http：/www.ncbi.nlm.nih.gov/go’s/，记录号GDS4336)，数据集概况见表1.经过基因芯片表达值正态化后，计算每个数据中的每一对基因其皮尔森相关系数.皮尔森相关系数数值介于-1至1.1代表基因对表达呈正相关，-1则代表负相关.对各个数据集中基因对皮尔森相关系数从大到小排序后，挑选出排名前1%的基因对定义为表达正相关，排名后1%则代表基因对表达值负相关.合并两个数据集中筛选出的具有强烈表达相关性的基因对后，将其对应到STRING蛋白质互作网络中，从而获得一个由筛选出基因对所构建的关于PDAC的功能性相关网络.其中利用STRING蛋白质互作网络进行功能性相关评估，选取分数大于500的基因对去绘制此网络.经过上述步骤，本研究构建出一个PDAC功能相关的基因网络，此网络包含4460个基因作为网络节点和 24633 条边.

2.2 划分PDAC功能相关基因网络中的癌症相关基因模块

本研究运用iNP基因模块划分算法划分上述分析方法所获得的PDAC功能相关基因网络.通过基因模块划分的方法，可以更加清晰地将一个大的基因网络划分为小的基因富集的模块，从而有利于基因表达数据的分析.通过iNP算法，本研究所获得的PDAC功能相关基因网络可以被划分为748个基因模块，其模块化选取阈值为0.556(一般阈值为0.3即可认为是基因模块).每个模块中的基因数目大于3个对于后续分析才有意义，经过筛选，所得基因模块个数为452个，各个模块中基因数目分布情况见图1.可以看出，90%以上的模块所含有的基因数目在50个以下.

对于每一个划分出的基因模块，其各模块内部所构成基因相互间都有极强的表达相关性.为了检验一个基因模块在癌症和正常样本中表达值是否显著区分，我们采用中值差异显著性检验方法.每个模块中的基因都可以找到其在45个PDAC样本和45个正常组织中的表达值，计算其在癌症和正常样本中的中值对，运用t检验进行统计显著性分析.

由于每个模块包含基因众多，要使模块统计显著性满足条件，本研究通过对伪发现率的控制在0.05以内，计算并定义出P<0.01才构成统计显著性.除此外，其还要满足表达值的log10值转换大于2这一筛选条件.经过筛选，本研究筛选出在癌症样本和正常组织样本中具有表达显著性区别的基因模块个数为34.在这34个基因模块中，对比于正常组织样本，其中20个模块在PDAC样本中显著上调，14个相对显著下调.

利用R语言中pheatmap包可以将模块差异表达情况显示如图2.图2左边14个列代表代表在癌症样本中基因表达下调的模块，右边20列则代表在癌症样本中基因表达下调的基因模块.所有模块都是经过P值筛选检验后具有显著性的模块.图2中右侧颜色渐进条表示各个模块的表达中值.

2.3 PDAC差异表达基因模块的功能富集分析

上述结果所得到的34个在PDAC差异表达的基因模块中，我们想知道能否通过差异表达基因模块去区分癌症和正常样本.本研究利用模块表达值(模块中所有基因表达值的中值)去划分癌症样本和正常样本，并画出预测的ROC曲线图.ROC曲线图中曲线下方面积，即AUC的数值.其阈值为0到1，0.5代表随机预测.愈接近0或1则代表可以该模块在癌症样本和正常样本中差异表达愈显著.本研究中接近1则代表相对于正常组织样本，该模块在PDAC中上调，反之接近0，则代表下调.上调的20个模块中，其AUC中值为0.8315，其中有3个模块AUC值高于0.9.这表示这3个模块在癌症样本中表达值显著的上调.下调的14个模块中，AUC模块中值为0.21，其中最低AUC为0.117.这表示该模块在癌症样本中表达值显著的下调.上调模块和下调模块AUC数值分布见图3.

对模块整体分析后，本研究进一步挑选出这些模块中PDAC的显著表达差异基因，作为gene biomarker来帮助临床PDAC的鉴定.比如，两个互作基因A和B，如果在癌症样本中表达正相关，然而在正常样本中表达负相关，则共表达差异值为1-(-1)=2.如果两个互作基因在样本中相关性随机，则赋值1.如果在两种样本中都正相关或者负相关，则为0.通过计算模块中所有互作共表达差异值，我们可以通过计算均值来代表该基因在模块中的共表达差异值.最终用该数值乘以所属模块在癌症中的特异性便可计算该基因在癌症样本中的表达特异性.通过上述方法，我们筛选出上调模块中333个差异表达基因，和下调模块中270个差异表达基因.我们选取上调模块和下调模块中前10个差异表达的基因做重点分析.

如图4所示.条形图中灰色条长度代表每个基因能够区分癌症和正常样本的分数高低.基因分数大于0则代表该基因在癌症样本中表达上调，分数小于0则代表该基因在癌症样本中表达下调.各基因所属模块标注于灰色条形部分.

从挑选出的模块及其对应基因来看(图4)，其中2个上调基因来自M478模块.图5中可以看出该模块在癌症样本和正常样本中的基因互作网络，以及对应模块AUC的数值.图5中网络节点为基因，网络的边为互作PCC数值，网络边的粗细则代表PCC数值的大小.可看出M478模块AUC高达0.926，证明该对比与正常组织样本，该模块在胰腺癌样本中显著上调.其中DMBT1基因在张光斌等[20]的研究中证实，其过表达对癌细胞系EC9706具有生物学行为的影响.

其中1个上调的基因来自M57模块.该模块富集了表皮生长因子相关功能.该模块AUC的数值高达0.881，证明该对比与正常组织样本，该模块在胰腺癌样本中显著上调.其中DSC3基因打分最高，在正常样本中DSC3基因同DSG3基因表达呈现正相关.然而在如图5(见第646页)所示，在PDAC样本中DSC3基因和DSG3基因并无互作.这说明，DSC3基因很有可能通过打乱表皮生长因子相关功能而引起PDAC的产生.Hamidov等[21]研究结果显示，DSC3作为15个通过临床病理学筛选出的gene biomarker中的一个基因，其过量的表达可以通过影响细胞间粘附而对肿瘤产生产生影响.而在下调模块的基因当中M466中，其模块AUC为0.156，说明该模块在癌症样本中显著下调.M466模块中的DLG5基因，在2005年Taniuchi[22]等人的研究中发现，通过小RNA干扰RAB6KIFL基因在PDAC细胞系中的表达，细胞系中的细胞大幅下降.而DLG5作为蛋白质复合物的连接体，其与RAB6KIFL基因互作可以对胰腺癌的发生产生影响.

3 结 论

基因表达谱芯片为癌症基因检测带来极大便捷.通过对基因表达芯片的数据分析，比对癌症样本和正常样本中基因表达差异，所筛选的基因生物标识对癌症的致病机理以及辅助临床诊断具有重要的参考作用.在本研究中，我们通过计算筛选出能够用来区分PDAC和正常组织的基因生物标记，可以为临床检测作参考.通过整理在癌症样本中和正常组织样本中具有明显正相关和负相关共表达基因对，以及评估它们的功能相关性，我们构建出一个PDAC的特异性基因互作网络.通过网络模块识别划分算法，将该网络划分为若干基因模块.在每一个模块中，模块内基因都具有共表达特性和模块功能相关性.通过分析各个模块在癌症样本和正常样本中的表达中值，我们可以筛选出具有癌症表达特异性的基因模块.我们筛选出了34个基因模块，其中20个为在癌症样本中明显上调的模块，另外14个为在癌症样本中表达值明显下调的模块.对筛选出的每个模块中差异表达的基因进行整合，得到20个可以作为生物标记的基因，它们可以用来识别PDAC样本和正常样本.这些标记可以为PDAC发生分子机制的研究、早期诊断和分子靶向治疗提供线索.

[1] JEMAL A， SIEGEL R， XU J，etal. Cancer Statistics [J].CA：ACancerJournalforClinicians， 2010，60(5)：277-300.

[2] YEH J J. Prognostic signature for pancreatic cancer：Are we close? [J].FutureOncol， 2009，5(3)：313-321.

[3] COLLISSON E A， SADANANDAM A， OLSON P，etal. Subtypes of pancreatic ductal adenocarcinoma and their differing responses to therapy [J].NatMed， 2011，17(4)：500-503.

[4] OETTLE H， POST S， NEUHAUS P，etal. Adjuvant chemotherapy with gemcitabine vs observation in patients undergoing curative-intent resection of pancreatic cancer：A randomized controlled trial [J].JAMA， 2007，297(3)：267-277.

[5] CONROY T， DESSEIGNE F， YCHOU M，etal. FOLFIRINOX versus gemcitabine for metastatic pancreatic cancer [J].NEnglJMed， 2011，364(19)：1817-1825.

[6] GOGGINS M. Identifying molecular markers for the early detection of pancreatic neoplasia [J].SeminOncol， 2007，34(4)：303-10.

[7] KOLBERT C P， KOLBERT C P， CHARI S，etal. Microarray technologies for gene transcript analysis in pancreatic cancer [J].TechnolCancerResTreat， 2008，7(1)：55-59.

[8] GRUTZMANN R， BORISS H， AMMERPOHL O，etal. Meta-analysis of microarray data on pancreatic cancer defines a set of commonly dysregulated genes [J].Oncogene， 2005，24(32)：5079-5088.

[9] CAMPAGNA D， COPE L， LAKKUR S S，etal. Gene expression profiles associated with advanced pancreatic cancer [J].IntJClinExpPathol， 2008，1(1)：32-43.

[10] KIM H N， CHOI D W， LEE K T，etal. Gene expression profiling in lymph node-positive and lymph node-negative pancreatic cancer [J].Pancreas， 2007，34(3)：325-334.

[11] STRATFORD J K， BENTREM D J， ANDERSON J M，etal. A six-gene signature predicts survival of patients with localized pancreatic ductal adenocarcinoma [J].PLoSMed， 2010，7(7)：e1000307.

[12] SUN S， DONG X， FU Y，etal. An iterative network partition algorithm for accurate identification of dense network modules [J].NucleicAcidsRes， 2012，40(3)：e18.

[13] DAVIS S， MELTZER P S. GEOquery：A bridge between the Gene Expression Omnibus(GEO) and BioConductor [J].Bioinformatics， 2007，23(14)：1846-1847.

[14] BARRETT T， TROUP D B， WILHITE S E，etal. NCBI GEO：Mining tens of millions of expression profiles—database and tools update [J].NucleicAcidsRes， 2007，35(Database issue)：760-765.

[15] ZHANG G， SCHETTER A， HE P，etal. DPEP1 inhibits tumor cell invasiveness， enhances chemosensitivity and predicts clinical outcome in pancreatic ductal adenocarcinoma [J].PLoSOne， 2012，7(2)：e31507.

[16] SMYTH G K. Linear models and empirical bayes methods for assessing differential expression in microarray experiments [J].StatApplGenetMolBiol， 2004，3：Article3.

[17] JENSEN L J， KUHN M， STARK M，etal. STRING 8-a global view on proteins and their functional interactions in 630 organisms [J].NucleicAcidsRes， 2009，37(Database issue)：412-416.

[18] ASHBURNER M， BALL C A， BLAKE J A，etal. Gene ontology：tool for the unification of biology. The Gene Ontology Consortium [J].NatGenet， 2000，25(1)：25-29.

[19] LIBERZON A， SUBRAMANIAN A， PINCHBACK R，etal. Molecular signatures database(MSigDB) 3.0 [J].Bioinformatics， 2011，27(12)：1739-1740.

[20] 张光斌，贺 涛，张 宁.DMBT1过表达对食管癌细胞系EC9706生物学行为的影响 [J].世界华人消化杂志，2009，17(17)：5.

[21] HAMIDOV Z， ALTENDORF-HOFMANN A， CHEN Y，etal. Reduced expression of desmocollin 2 is an independent prognostic biomarker for shorter patients survival in pancreatic ductal adenocarcinoma [J].JClinPathol， 2011，64(11)：990-994.

[22] TANIUCHI K， NAKAGAWA H， NAKAMURA T，etal. Down-regulation of RAB6KIFL/KIF20A， a kinesin involved with membrane trafficking of discs large homologue 5， can attenuate growth of pancreatic cancer cell [J].CancerRes， 2005，65(1)：105-112.

Identification of Gene Biomarkers for Distinguishing Pancreatic Ductal Adenocarcinoma from Normal Tissue Using a Network-based Approach

WANG Xinge， TIAN Weidong

(Department of Biostatistics and Computational Biology， School of Life Sciences，FudanUniversity，Shanghai200438，China)

Pancreatic ductal adenocarcinoma(PDAC) is one of the most lethal cancers in the world. In this study， we have developed a network-based approach to identify gene biomarkers that can distinguish PDAC and normal tissues. By identifying strongly coexpression gene pairs in normal tissues and PDAC samples and applying functional association information from STRING network online， a PDAC-specific gene association network was constructed. Based on this network， gene modules was obtained in which genes are highly functionally associated using iNP algorithm. Then gene modules which are differentially expressed between PDAC and normal samples were identified. Finally， genes inside those modules were selected with discriminating coexpression patterns between PDAC and normal samples and those genes can be used as candidate biomarkers to facilitate clinicopathologic diagnose.

pancreatic ductal adenocarcinoma； gene biomarker； gene coexpression network； function enrichment analysis

0427-7104(2016)05-0642-07

2016-01-14

王新格(1989—)，女，硕士研究生；田卫东，男，教授，通讯联系人，E-mail：weidong.tian@fudan.edu.cn.

Q 332

A