结直肠腺瘤的加权基因共表达网络构建与分析

高亚东，屈亚威，刘海峰

作为一种常见的消化道恶性肿瘤，结直肠癌（colorectal cancer，CRC）的发病率及病亡率近年来在全球范围内呈上升趋势。其中绝大多数散发的CRC是由其重要的癌前病变即结直肠腺瘤（colorectal adenoma，CRA）发展而来[1]，因此探究CRA的发生发展机制显得尤为重要。而CRA的发生发展机制较为复杂，一般涉及多个基因间的相互作用，传统的仅关注单个基因或者少数几个基因的方法由于不能提供基因表达全景信息，极大地限制了疾病发生发展机制的研究。下一代测序技术的出现实现了可同时研究多个基因表达情况的功能。例如转录组测序、基因芯片等技术已经被广泛用于识别与CRA发生发展相关的生物标志物[2,3]，然而，基于这些技术的相关研究多数关注寻找差异基因，忽略了基因间可能的高度相关性，而这些具有相似表达模式的基因往往在功能上也是相关的。基于此，本研究采取一种称之为加权基因共表达网络分析（weighted gene co-expression network analysis，WGCNA）[4]的方法对来自于基因表达综合数据库（gene expression omnibus，GEO）的CRA基因表达谱进行分析，构建CRA的基因共表达网络，识别共表达模块，并结合临床信息探讨其与模块间的相关性，最后筛选出与临床信息密切相关的枢纽基因。

1 材料与方法

1.1 数据来源及数据预处理 本研究基因表达谱数据（GSE41657）及其临床信息来源于美国国家生物信息中心下的GEO数据库，平台编号为GPL6480，包含12个正常结直肠黏膜组织、51个结直肠腺瘤组织及25个腺癌组织样本基因表达谱数据，本文选取其中51个结直肠腺瘤组织样本基因表达谱数据信息用于分析。在进行WGCNA前，本研究对下载的矩阵信息进行了预处理：利用平台信息将探针名转化为基因名，并将对应多个探针的基因表达量取平均值作为该基因的最终表达量。最终获得行名为样本名，列名为基因名的矩阵用于后续共表达网络的构建。

1.2 WGCNA

1.2.1 共表达网络的构建和模块的识别 在R语言下运行WGCNA包，进行共表达网络的构建及模块的识别。为获取有效的共表达网络，本研究选取表达量方差大于所有方差四分位数的基因。通过剔除离群样本，保证网络构建的结果是稳定的。在剔除离群样本后，定义基因共表达的相关矩阵。通过选择合适的加权系数β使构建的网络更符合无尺度特征[5]，并应用此系数将相关矩阵转换为邻接矩阵后进一步转化为拓扑重叠矩阵（topological overlap matrix，TOM），并计算相异度。采用动态剪切树法层次聚类基因，其中将相异度作为距离测度，最小模块尺寸为30，进行模块识别并绘制基因树状图。为了进一步分析识别出来的模块，本研究计算出这些模块的模块特征值（module eigengene，ME），并采用层次聚类对ME进行聚类，绘制模块树状图，并设置高度0.25作为分割线来合并相似程度较高的模块。

1.2.2 模块的功能富集分析 为了解鉴别出来的模块具有的生物学功能，本研究将各模块内的基因映射至在线网站DAVID（https：//david-d.ncifcrf.gov/）中，进行基因本体（gene ontology，GO）及京都基因与基因组百 科 全 书（Kyoto Encyclopedia of Genes and Genomes，KEGG）通路分析[6]，并设置P＜0.05的条目富集程度具有显著性意义。

1.2.3 与临床信息相关的模块及枢纽基因的鉴别 本研究使用两种方法来鉴别与结直肠腺瘤样本临床信息相关的模块。第一种方法，计算各模块的ME与外部各个临床信息的皮尔森相关系数及其P值，通过皮尔森相关系数衡量不同模块与各个临床信息间的关系，并取相关系数最高的模块用于后续分析。第二种方法，计算每个基因的表达量与某一临床信息的皮尔森相关系数，称之为基因显著性（gene significance，GS），随后对模块内所有基因的GS取绝对值均值即为该模块的模块显著性（module significance，MS），MS值用以衡量某一模块与某一临床信息的整体相关程度。此外，在筛选出目的模块后，为了筛选该模块中的枢纽基因，本研究计算了目的模块中每一基因与该模块ME的相关系数，即模块隶属度（module membership，MM），用以衡量目的模块内的某一基因从属于该模块的程度。本研究定义某基因满足GS绝对值大于0.2且其归属模块MM绝对值大于0.8则为枢纽基因。

1.3 统计学处理 本研究采用的结直肠腺瘤的相关临床特征和各个共表达模块的ME间相关系数的计算是基于R语言平台Rstuio下实现。

2 结 果

2.1 结直肠腺瘤组织样本基因表达谱数据的预筛选本研究从GEO数据库中下载得到行名为51个样本，列名为41 076个核酸探针的矩阵。在利用平台文件GPL6480将探针名转换为基因名并对其中对应多个探针的基因表达量取平均值后，最终得到包含19 595个基因的结直肠腺瘤组织样本基因表达谱矩阵，并作为共表达网络构建的输入文件。

2.2 WGCNA

2.2.1 共表达网络的构建和模块的识别 如方法学部分所示，本文计算了各基因在所有样本中表达量的方差，并取方差在前25%的基因用于共表达网络的构建，共4 899个基因。通过构建51个结直肠腺瘤样本的4 899个基因的层次聚类树后，设置层高阈值为105，去除离群样本4个，编号分别为GSM1021177、GSM1021179、GSM1021137、GSM1021141。 然 后对 剩下的结直肠腺瘤样本再次进行层次聚类，同时绘制样本信息（图1）。

为选择合适的加权系数β用于无尺度网络的构建，同时需考虑对各基因节点（Node）平均连接度的适度保留，本研究最终选取β=7用于构建共表达网络（图2）。在β值确定后，本研究通过动态剪切树法进行模块初步识别并合并相似模块，去除灰色模块（此模块包含无法分配至任何一个模块的基因），最终得到15个模块，如表1及图3所示。

图1 47个结直肠腺瘤样本和临床信息的系统聚类

图2 加权基因共表达网络分析软阈值（β）的确定

表1 构建出的模块资料

图3 基因聚类树状图

2.2.2 模块的功能富集分析 利用DAVID在线工具对识别出来的各模块中基因分别进行GO和KEGG分析。主要结果如下：蓝色模块中基因主要富集在乏氧应答、细胞对肿瘤坏死因子的反应、胞外组织、氧化磷酸化等分子行为过程中，并主要富集于帕金森病等通路中；蓝绿色模块中基因主要富集在血管生成、骨形态生成蛋白信号通路负性调节、中心体等分子行为过程中，并富集于轴突导向等通路中；绿色模块中基因主要富集在细胞分裂等过程，并主要富集于细胞周期等通路；黄绿色模块主要富集于内皮细胞分化等，并富集于癌症通路等。浅灰色模块主要富集于主要组织相容性复合物Ⅱ蛋白质复合体等过程；淡青色模块主要富集于葡萄糖醛酸基转移酶活性等过程，并富集于T细胞受体等信号通路；淡绿色模块富集于细胞体等过程；淡黄色模块在信号传导子等过程及细胞因子受体相互作用等通路富集；品红色模块在金属肽酶活动等过程显著富集，并在病毒致癌等信号通路富集。午夜蓝色模块主要富集于信号转导等过程，并在信号转导及转录激活因子途径等通路富集；粉红色模块在类固醇代谢、外泌体等过程显著富集；紫色模块显著富集于膜电位调节等过程；橙红色模块主要富集于RNA聚合酶Ⅱ启动子转录的正调控、脱氧核糖核酸（deoxyribonucleic acid，DNA）模板的转录等过程；蓝绿色主要富集于细胞信号传导等；黄色模块主要富集于氧化还原等过程，并在代谢途径等通路富集。

2.2.3 与临床信息相关的模块及枢纽基因的鉴别 通过对结直肠腺瘤样本临床信息和各个模块进行关联分析，可选出与临床信息显著相关的目的模块。本文主要关注与上皮内瘤变级别显著相关的模块。如图4所示，共有8个模块与上皮内瘤变级别显著相关，其中淡绿色模块相关程度最高。为增加此结果的可信度，通过计算各模块MS的方式再次验证此结果，如图4b所示，淡绿色模块在所有模块中MS最大（P=4e-196）。此外，通过计算淡绿色模块中基因GS和MM相关系数（cor=0.65，P=4.5e-08）亦进一步验证此结果。基于上述结果可以得出，淡绿色模块和CRA样本上皮内瘤变级别的关联程度最大，本文后续对枢纽基因的筛选选择此模块。通过计算淡绿色内基因的GS和MM后，获得16个基因，排除其中5个没有注释的基因，最终得到11个枢纽基因，如表2所示。

3 讨 论

图4 结直肠腺瘤与不同临床特征相关模块的鉴定

表2 11个枢纽基因的主要信息

作为结直肠癌的癌前病变，CRA的发生发展机制较为复杂，涉及多种分子及通路。传统的仅针对单个或少数几个分子的研究无法对CRA发生发展机制进行全面探索。而WGCNA作为一种复杂的基因共表达网络构建方法，在处理多样本的复杂数据上有其独特优势[7]。

WGCNA通过在高通量数据中识别由功能相关的基因构成的模块，从生物功能整体考虑基因功能及其联系，弥补了传统方法的缺陷。通过此分析方法，研究者不仅可以发现某模块内基因间的联系，还可以关注该模块与其他模块内基因的相关性。此外，通过将临床信息与模块相关联，还可进一步获得与临床特征相关的基因，有助于建立疾病某临床特征研究基础。既往已有报道采用此方法对乳腺癌、肺癌等[8,9]恶性肿瘤进行分析，但尚缺乏应用WGCNA对结直肠腺瘤进行分析。本文通过采用WGCNA对GEO数据库中的51个CRA样本基因表达谱进行分析，最终获得了15个有生物学功能的模块。

通过对蓝色模块进行功能分析后发现，该模块内的基因在乏氧应答过程条目中显著富集。缺氧是肿瘤微环境的基本特征之一，过去较多研究显示缺氧与包括结直肠腺瘤及腺癌在内的各种肿瘤的发生发展有着密切的关系[10,11]。Giatromanolaki等[7]通过检测120个伴不同类型管状腺瘤样本的缺氧诱导因子1α、自噬基因LC3A和Beclin-1的表达水平后发现，这些腺瘤样本中缺氧引起的自噬与腺瘤发育程度和侵袭性有关。结合本研究的结果，笔者推断蓝色模块内的基因可能通过乏氧诱导参与腺瘤-癌发展过程。蓝绿色模块功能富集显示该模块基因在骨形态生成蛋白信号通路负性调节过程中显著富集。作为转化生长因子-β家族成员的一个亚组成员，骨形态生成蛋白作为肿瘤抑制因子在CRC形成的早期阶段起着重要作用[12]。Freeman等[13]发现骨形态生成蛋白信号通路在Smad4介导下抑制β-链蛋白的表达，继而抑制腺瘤发展到腺癌的过程。因此，对骨形态生成蛋白信号通路的各个靶点进行干预可能成为治疗早期腺癌的新思路，而富集于该通路条目的基因可作为下一步研究的重点。此外，蓝绿色模块中基因模块在中心体条目显著富集。既往研究显示，中心体数目和结构的改变在结直肠癌腺瘤-癌过程中十分常见，并早在低级别上皮内瘤变阶段已经出现[14]。结合本研究模块的富集信息，笔者推测中心体可能在腺瘤发展成腺癌的过程中发挥重要作用，而通过研究中心体在腺瘤-癌发展过程中的动态变化可能有助于深层次理解这种癌变过程的机制。

尿苷二磷酸葡萄糖醛酸转移酶作为一种重要的代谢酶，在部分体内内源性和外源性毒性物质代谢过程中起着关键的催化作用。Hoensch等[15]发现和腺瘤样本相比，结直肠癌样本中该解毒酶的代谢能力较低。本研究结果显示，葡萄糖醛酸基转移酶活性是淡青色模块基因重要富集条目之一。因此，葡萄糖醛酸转移酶可能在CRA发生发展过程中发挥重要作用，而提高尿苷二磷酸葡萄糖醛酸转移酶活性可能对于预防和治疗结直肠癌具有一定意义。黄色模块功能富集分析显示该模块内基因显著富集在脂质代谢过程条目。脂质代谢和腺瘤的发生具有密切关联。Wang等[16]发现低密度脂蛋白胆固醇通过激活活性氧水平和丝裂原活化蛋白激酶通路促进炎性反应和结直肠癌的发展。然而这种途径是否在腺瘤形成过程中产生作用尚未完全明确，关于脂质代谢在CRA形成过程中发挥的具体机制仍有待探究。

此外，通过与临床信息相关联，笔者所在团队鉴别出8个与CRA样本上皮内瘤变级别显著性相关的模块，其中淡绿色模块相关性最高，通过计算模块内基因的GS和MM后，笔者团队进一步识别出11个可能在上皮内瘤变级别改变过程中发挥关键作用的枢纽基因。其中部分基因在结直肠肿瘤领域的研究已有报道。例如，GBP6被认为在CRA发展过程中发挥重要作用[17]。GAGE7作为G抗原的重要成员，具有抗凋亡能力，在包括结直肠癌在内的众多恶性肿瘤中发挥重要作用[18,19]。本研究结果显示，淡绿色模块内基因与CRA样本上皮内瘤变级别呈显著负相关，而作为可能在该模块内发挥重要作用的枢纽基因GBP6及GAGE7等，则可能在腺瘤发展过程中起着关键的负性调控作用，这种负性调控作用的具体机制有待于进一步研究。此外，该模块中其他大多数枢纽基因在CRA及CRC研究中尚未见报道，其作用需要在未来细胞及动物实验中进一步验证。

总之，本研究通过WGCNA构建结直肠腺瘤的基因共表达网络，识别出15个共表达模和11个可能在CRA发展过程中发挥关键作用的枢纽基因，这些模块及基因对深层次理解CRA的分子机制可起到一定帮助。