CN109872772B - 利用权重基因共表达网络挖掘结直肠癌放疗特异性基因的方法 - Google Patents

Abstract

本发明涉及权重基因共表达网络，特别涉及一种利用权重基因共表达网络挖掘结直肠癌放疗特异性基因的方法。本发明利用权重基因共表达网络挖掘结直肠癌放疗特异性响应基因的方法，通过权重基因共表达网络为挖掘结直肠癌放疗特异性响应基因提供了新途径，为预测结直肠癌患者生存预后提供了新的依据。针对本发明的方法发现的靶点基因的放疗可以在一定程度上提高结直肠癌患者的生存预后，降低结直肠癌患者死亡率，解决实际的临床问题，为广大患者提供了更佳的选择。

Description

利用权重基因共表达网络挖掘结直肠癌放疗特异性基因的 方法

技术领域

本发明涉及权重基因共表达网络，特别涉及一种利用权重基因共表达网络挖掘结直肠癌放疗特异性基因的方法。

背景技术

结直肠癌是中国及全球高发的恶性肿瘤，约占所有病例的10％，其死亡率高、易复发、预后差，是一种由多种癌基因和抗癌基因共同参与，多途径协同作用导致的基因病。全球疾病负担(Global Burden of Disease ，GBD)数据显示，2016年全球患有结直肠癌的人数达到632万，其中中国患病人数高达147万。2016年全球结直肠癌患者的死亡人数为83万，其中中国死亡患者数超过16万，占总死亡人数的1.73％。统计结果显示，从1990年到2016年全球结直肠癌患病率和死亡率持续增长，中国的患病率和死亡率增长较快，尤其是近几年中国的患病率和死亡率均超过全球平均水平。

随着全基因组水平检测技术，如GWAS、表观遗传学、第二代测序技术、高精度高通量质谱仪和代谢组学等的发展，检测成本不断降低，海量组学数据不断增多，传统的数据分析方法已不能满足分析需求。而权重基因共表达网络分析能对复杂的基因组水平的数据进行降维，将其简化为若干个模块。因此，权重基因共表达网络分析在多个研究领域得到了应用，例如：Iancu等首次将权重基因共表达网络分析应用于RNA-Seq数据，Shirasaki等首次将其应用于蛋白质组数据。

在对疾病的研究中，构建权重基因共表达网络是挖掘目的基因及其表达产物的功能的有效方法,是理解癌症生物学和肿瘤发展机制所必需的强有力的“数据驱动力”。目前已探索了许多与结直肠癌相关的基因及其蛋白质产物，但少有研究利用权重基因共表达网络探究癌症放疗前后的差异基因。而精确的短期放疗在直肠癌的治疗中起重要作用，挖掘癌症放疗的特异性响应基因可以大大提高放疗的效果。因此，利用权重基因共表达网络挖掘直肠癌放疗特异性响应基因，确定放疗特异性分子靶点，阻止癌症的发展，提高结直肠癌患者的生存预后显得尤为重要。

发明内容

鉴于现有技术存在的问题，本发明的目的在于提供一种利用权重基因共表达网络挖掘结直肠癌放疗特异性基因的方法，本发明构建权重基因共表达网络挖掘结直肠癌放疗特异性响应基因以更精确地指导临床治疗，阻止癌症的发展，提高结直肠癌患者的生存预后，降低结直肠癌患者的死亡率。

为了实现上述目的，本发明采用一下技术方案。

利用权重基因共表达网络挖掘5个结直肠癌放疗特异性响应基因，分别为CCNE1、CDT1、MCM6、DBF4、CENPK。

一种利用权重基因共表达网络挖掘结直肠癌放疗特异性响应基因的方法，包括以下步骤。

步骤1、下载GEO的基因表达谱数据，将数据标准化后，使用方差分析筛选出在放疗前后在癌细胞与癌旁细胞中表达存在差异的基因，使用WGCNA对这些基因构建加权共表达网络，挖掘共表达模块，并对每一个模块做KEGG富集分析，观察每个模块的功能。

步骤2、计算每个模块与癌症放疗的关系，筛选出与癌症放疗最相关的模块，提取出这些模块的共表达网络，整合人类蛋白质互作网络，构建共表达-互作子网，分析子网中的基因，最终筛选出最合适的目的基因。

步骤3、为证明结果的有效性，再次对目的基因进行KEGG通路分析以及利用gepia在线工具定制化分析关键基因在癌细胞与癌旁细胞中的表达变化。

进一步地，一种利用权重基因共表达网络挖掘结直肠癌放疗特异性响应基因的方法，所述步骤1具体包括以下步骤。

1）GEO数据下载和数据预处理：从GEO数据库下载GSE15781结直肠癌的芯片。为了标记芯片中的基因，需要插入相应的基因探针。首先将探针插入到对应基因上，去除掉空载的探针。然后下载数据集Quantile normalized signal intensity。当多个探针对应一个基因时取中位数作为该基因的表达水平。

2）筛选在各类样本中差异的基因：使用方差分析计算每个基因在各类样本中的表达情况以筛选出各类样本中差异的基因，选择p值小于0.05的基因作为后续分析的目标基因。

3）基因共表达网络构建：WGCNA是使用基因表达数据来构建无尺度网络的系统生物学方法。其基本思路如下：首先，构建基因表达相似性矩阵，即计算两两基因之间皮尔森相关系数的绝对值。使用公式1计算基因i和基因j之间的皮尔森相关系数，其中i和j分别是第i个基因和第j个基因的表达量。

然后，使用公式2将基因表达相似性矩阵转换成邻接矩阵，网络类型为signed.。其中β为软阈值，即每对基因的皮尔森相关系数的β次方。这一步能够从指数级别强化强相关性、减弱弱相关性。

下一步使用公式3将邻接矩阵转换成拓扑矩阵，拓扑重叠(topological overlapmeasure, TOM)可以用来描述基因之间的关联程度。

/>

1-TOM表示基因i和基因j之间的相异程度。先使用1-TOM作为距离对基因进行层次聚类，再使用动态剪切树的方法对模块进行识别。每个模块中最具有代表性的基因称为特征向量基因，简称ME，它代表了该模块内基因表达的整体水平，是每个模块中的第一主成分。使用公式4来计算ME，其中i表示模块q中的基因，l表示模块q中的芯片样本。

用某个基因在所有样本中的表达谱与特征向量基因ME表达谱的皮尔森相关性来衡量这个基因在该模块中的身份，即模块身份(module membership)，简称MM。使用公式5计算MM，其中

是第i个基因的表达谱，表示模块q的特征向量基因(ME)，/>

表示了基因i在模块q中的身份，当/>

= 0，则说明基因i不在模块q中，/>

越接近+1或−1，则说明基因i与模块q高度相关。正负号表示了基因i与模块q是正相关还是负相关。

在计算基因与模块之间的相关性时，还需要考虑基因与外部因素的相关性。基因显著性(gene significance)，简称GS，用来衡量基因与外部信息的关联程度，GS越高表示基因越具有生物学意义，GS = 0，说明这个基因不参与所研究的生物学问题。

将在各个样本中筛选出的表达有差异的基因的表达数据构建权重共表达网络，选择软阈值为9，筛选共表达模块。

4）共表达模块的富集分析：为观察各个共表达模块的功能，使用R软件包clusterProfiler对各个模块进行KEGG通路富集分析。

进一步地，一种利用权重基因共表达网络挖掘结直肠癌放疗特异性响应基因的方法，所述步骤2具体包括以下步骤。

1）各模块与癌症样本相关性分析：根据放疗样本（放疗，癌症样本放疗，正常样本放疗）类型，将对应属性的样本定义为1，其他样本定义为0，所有样本定义值形成一个表型矩阵，计算每个模块与表型矩阵中各个因素的相关性，筛选出与癌症放疗相关的模块。

2）癌症相关的模块的权重共表达网络构建：根据各个共表达模块中的基因的表达关系，选择共表达权重大于0.2的共表达模块对作为最终的共表达网络的边，导出癌症相关的模块的共表达网络。

3）共表达-互作调控网络构建：参与相关功能的基因在蛋白互作网络中通常为距离较近的节点，即在蛋白互作网络中表现出临近关系。为了识别出癌症预后显著相关的重要基因，从ps://www.ncbi.nlm.nih.gov/pmc/articles/PMC5210659/下载HIPPIE v2.0，将得到的癌症相关模块的共表达网络映射到人类蛋白互作网络中，构建共表达-互作特异性网络。

4）共表达-互作特异性网络度分布分析：生物学网络是复杂网络，往往具有自组织、自相似、吸引子、小世界、无尺度等特征，对共表达-互作特异性网络进行分析，观察共表达-互作特异性网络是否符合无尺度特征。

5）富集共表达基因：共表达基因往往和其周围临近基因共同参与相同的疾病通路或生物学过程，核心蛋白和其临近蛋白在表达量上也具有显著相关性，因此通过统计每个基因周围共表达的基因个数和连接的所有基因个数，以整个人类蛋白互作网络作为统计背景，构建统计学模型。对每个基因做Fisher富集检验，筛选显著富集FDR<0.01且共表达且互作的基因占所有与其互作的基因比例大于0.1的网络节点，最终得到目的基因（癌症放疗响应特异性基因）。

6）癌症放疗响应特异性基因的表达谱分析：为观察这些基因之间的表达关系，运用表达谱层次聚类分析的方法。

进一步地，一种利用权重基因共表达网络挖掘结直肠癌放疗特异性响应基因的方法，所述步骤3具体包括以下步骤。

1）癌症放疗响应特异性基因的富集分析：为观察这些基因的功能，对这些基因做KEGG Pathway富集分析。

2）癌症放疗响应特异性基因在外部数据集中的表达关系：通过gepia在线工具分析癌症特异性基因在TCGA数据集、GETx数据集中的表达关系。

与现有技术相比，本发明具有以下有益效果。

本发明利用权重基因共表达网络挖掘结直肠癌放疗特异性响应基因的方法，通过权重基因共表达网络为挖掘结直肠癌放疗特异性响应基因提供了新途径，为预测结直肠癌患者生存预后提供了新的依据。针对上述方法发现的靶点基因的放疗可以在一定程度上提高结直肠癌患者的生存预后，降低结直肠癌患者死亡率，解决实际的临床问题，为广大患者提供了更佳的选择。

附图说明

图1是模型构建的流程图。

图2是表达差异较大的基因的权重共表达网络构建过程中涉及到的图片。其中A、B图为各种软阈值功率的网络拓扑分析；C图为基因树状图和模块颜色；D图为模块-特征的相关性。

图3是各个模块之间的关系以及与放疗、肿瘤放疗、正常放疗三种方式之间的关系。

图4是7个模块富集到的KEGG Pathway的联系。

图5是共表达网络度分布。

图6是共表达-互作子网网络的度分布。

图7是放疗响应特异基因调控网络。

图8是肿瘤特异性基因的KEGG富集结果。

图9是5个肿瘤放疗特异基因的表达聚类分析结果，基因聚类使用欧氏距离。

图10是5个肿瘤放疗特异表达基因的差异表达分析。

具体实施方式

下面结合具体实施例，进一步阐述本发明。应理解，这些实施例仅用于说明本发明而不用于限制本发明的范围。此外应理解，在阅读了本发明讲授的内容之后，本领域技术人员可以对本发明作各种改动或修改，这些等价形式同样落于本申请所附权利要求书所限定的范围。

实施例 利用权重基因共表达网络挖掘结直肠癌放疗特异性响应基因的方法。

1. 数据预处理。

从GSE15781数据集中提取的直肠癌的芯片数据共有42个样本，其中包括癌症放疗样本9个、癌症未放疗样本13个、正常未放疗样本10个、正常接受放疗样本10个。插入基因探针，然后使用ABI Human Genome Survey Microarray Version 2芯片平台进行表达谱定量。为了明确基因筛选的流程，我们构建了分析流程图如图1所示。

2. 筛选差异基因。

根据筛选条件从16753个基因中共得到了7746个在各类样本中表达存在差异的基因，变化最大的前二十个基因，见表1。

3. 表达差异较大的基因的权重共表达网络构建。

在癌症样本和正常样本中存在许多表达有差异的基因，其中差异较大的基因是研究的关键，使用R语言中的WGCNA软件包对表达差异较大的基因进行基因网络的构建。为使构建的网络符合无尺度网络，出现连接度为k的节点的对数log(k)与该节点出现的概率的对数log(P(k))需要呈负相关，且相关系数要大于0.8。将表达矩阵转换成邻接矩阵，然后再将邻接矩阵转换成拓扑矩阵，设定β = 9（如图2A和2B所示)。然后使用average-linkage层次聚类法对基因进行聚类，按照混合动态剪切树的标准，设置每个基因网络模块最少的基因数目为30。在使用动态剪切法确定基因模块后，依次计算每个模块的特征向量值(eigengenes)，然后对模块进行聚类分析，将距离较近的模块合并成新的模块，设置height=0.25。共得到了10个模块，需要指出，grey模块是无法聚集到其它模块的基因集合(如图2C所示)。

各个模块中的基因统计如表2所示。7529个基因被分配到10个模块中且turquoise模块基因数显著高于其他模块。

计算每个模块的ME与放疗样本的皮尔森相关系数，相关系数越高说明这个模块越重要。如图2D所示，图2D中行表示每个模块的特征向量基因，列表示样本特征信息。从红色到绿色表示相关系数从高到低依次递减。每个小格子里的数字表示基因模块与相应特征的相关系数，括号中的数字表示P值。从图中可以看出与放疗显著正相关的模块为turquoise、brown、magenta，显著负相关的模块为green、blue、red、pink；与癌症样本放疗显著正相关的模块为turquoise、purple、magenta、yellow，显著负相关的模块为red、green、blue；与正常样本放疗并没有显著正相关的模块，显著负相关的模块为yellow、turquoise、pink、purple。

4. 模块及表型的关系。

为观察各个模块之间的相关性以及与放疗、肿瘤放疗、正常放疗三种方式的相关性，计算每个模块的特征向量，使用每个模块的特征向量聚类分析各个模块及放疗样本之间的关系。如图3所示，从图中可以看出这些模块之间大体上可以分三类，其中放疗和正常样本放疗与模块brown最为相近，癌症放疗与turquoise模块最相近(如图3所示)。

5. 富集分析。

如图4所示，不同的模块富集到不同的通路中，他们之间共同的通路很少。与癌症放疗最相关的turquoise模块中含有28个KEGG通路，其中包含多个与癌症相关的通路比如Pathways in cancer、Small cell lung cancer、Prostate cancer等，turquoise模块还与其他四个模块有共同通路。

6. turquoise模块的共表达网络构建。

根据各个共表达模块中的基因的表达关系，选择共表达权重大于0.2的基因对作为最终的共表达网络的边，导出癌症相关的模块的共表达网络。共表达网络中共包含563个基因。度分布如图5所示，从图中可以看出大部分基因的度都较低，只有少部分基因的度很高，这提示了这些度很高的基因在整个网络中具有重要的作用。

7. 构建共表达-互作调控网络。

从HIPPIE数据库下载人类蛋白互作数据构建人类蛋白质互作网络，网络中共包含17381个节点，平均邻居节点个数为19.6个。将375个基因全部映射到人类蛋白质互作网络中，结果显示具有相互作用且存在共表达的基因共有236个，平均邻居节点个数为3.26个。网络中的度分布如图6所示，从图中可以看出随着节点度越来越高，节点个数越来越少。

8. 放疗响应特异基因调控网络。

经过结合背景PPI网络和对网络的拓扑性质分析，可以计算共表达-互作子网中每个基因的共表达基因个数，构建统计学模型。然后计算网络中每个基因的共表达基因的富集情况，将其作为每个基因的网络属性，最终网络图如图7所示。从图中可以看出度最大且富集分数最高的前三个基因分别为CDH1、ATP1A1、ERBB2,这几个都已经被报道与癌症相关。

9. 放疗响应特异基因筛选。

根据网络中各基因的拓扑性质和共表达基因的富集显著性，最终确定筛选邻居节点中共表达基因的占比大于10%，并且在共表达-互作基因富集的结果中显著性p值的FDR小于0.01，同时邻居节点个数大于5的基因作为癌症放疗响应表达特异基因。最终结果见表3，从表中可以看出共得到5个基因。

10. 5个肿瘤放疗响应特异性基因的KEGG富集分析。

为观察这5个基因的功能，进一步分析这5个基因所参与的通路，使用R包clusterProfiler进行富集分析，最终结果如图8所示。从图中可以看出这5个基因富集到了2个KEGG通路中，其中有三个基因与Arrhythmogenic right ventricular cardiomyopathy(ARVC)通路有关，两个基因与蛋白质的消化和吸收有关。

11. 5个肿瘤放疗特异基因的表达谱聚类分析。

筛选出这5个基因的表达谱，做层次聚类分析结果如图9所示。从图中可以看出这5个基因之间的相关性较高，在癌症放疗样本中表达水平较低。

12. 5个肿瘤放疗特异基因的在TCGA数据中的表达分析。

为进一步分析这5个肿瘤放疗特异基因在结直肠癌中的表达变化，使用gepia（http://gepia.cancer-pku.cn/）在线工具分析在TCGA RNA-Seq数据集及GTEx数据集中的这5个基因的表达差异。如图10所示，从图中可以看出这5个基因的在癌症样本中都表现高表达。

通过上述具体的实施过程可以发现：本发明利用权重基因共表达网络挖掘了5个结直肠癌放疗特异性响应基因，针对这几个靶点基因的放疗可以在一定程度上提高结直肠癌患者的生存预后，降低结直肠癌患者死亡率，解决临床实际问题。

Claims (5)

1.利用权重基因共表达网络挖掘出CCNE1、CDT1、MCM6、DBF4、CENPK5个基因的组合在制备结直肠癌放疗特异性响应产品中的应用。

2.如权利要求1所述的应用，其特征在于，所述CCNE1、CDT1、MCM6、DBF4、CENPK5个基因的筛选方法具体包括以下步骤：

步骤1、下载GEO的基因表达谱数据，将数据标准化后，使用方差分析筛选出在放疗前后在癌细胞与癌旁细胞中表达存在差异的基因，使用WGCNA对这些基因构建加权共表达网络，挖掘共表达模块，并对每一个模块做KEGG富集分析，观察每个模块的功能；

步骤2、计算每个模块与癌症放疗的关系，筛选出与癌症放疗最相关的模块，提取出这些模块的共表达网络，整合人类蛋白质互作网络，构建共表达-互作子网，分析子网中的基因，最终筛选出最合适的目的基因；

步骤3、为证明结果的有效性，再次对目的基因进行KEGG通路分析以及利用gepia在线工具定制化分析关键基因在癌细胞与癌旁细胞中的表达变化；

筛选出5个结直肠癌放疗特异性响应基因分别为CCNE1、CDT1、MCM6、DBF4、CENPK。

3.如权利要求2所述的利用权重基因共表达网络挖掘结直肠癌放疗特异性响应基因的方法，其特征在于，所述步骤1具体包括以下步骤：

1）GEO数据下载和数据预处理：从GEO数据库下载GSE15781结直肠癌的芯片，为了标记芯片中的基因，需要插入相应的基因探针首先将探针插入到对应基因上，去除掉空载的探针；然后下载数据集Quantile normalized signal intensity；当多个探针对应一个基因时取中位数作为该基因的表达水平；

2）筛选在各类样本中差异的基因：使用方差分析计算每个基因在各类样本中的表达情况以筛选出各类样本中差异的基因，选择p值小于0.05的基因作为后续分析的目标基因；

3）基因共表达网络构建：WGCNA是使用基因表达数据来构建无尺度网络的系统生物学方法，首先构建基因表达相似性矩阵，即计算两两基因之间皮尔森相关系数的绝对值，使用公式1计算基因i和基因j之间的皮尔森相关系数，其中i和j分别是第i个基因和第j个基因的表达量；

然后使用公式2将基因表达相似性矩阵转换成邻接矩阵，网络类型为signed，其中β为软阈值，即每对基因的皮尔森相关系数的β次方，这一步能够从指数级别强化强相关性、减弱弱相关性；

下一步使用公式3将邻接矩阵转换成拓扑矩阵，拓扑重叠可以用来描述基因之间的关联程度；

公式3：

1-TOM表示基因i和基因j之间的相异程度，先使用1-TOM作为距离对基因进行层次聚类，再使用动态剪切树的方法对模块进行识别，每个模块中最具有代表性的基因称为特征向量基因，简称ME，它代表了该模块内基因表达的整体水平，是每个模块中的第一主成分，使用公式4来计算ME，其中i表示模块q中的基因，l表示模块q中的芯片样本；

公式4：

用某个基因在所有样本中的表达谱与特征向量基因ME表达谱的皮尔森相关性来衡量这个基因在该模块中的身份，即模块身份，简称MM；使用公式5计算MM，其中

是第i个基因的表达谱，表示模块q的特征向量基因(ME)，/>

表示了基因i在模块q中的身份，当

= 0，则说明基因i不在模块q中，/>

越接近+1或−1，则说明基因i与模块q高度相关，正负号表示了基因i与模块q是正相关还是负相关；

公式5：

在计算基因与模块之间的相关性时，还需要考虑基因与外部因素的相关性，基因显著性，简称GS，用来衡量基因与外部信息的关联程度，GS越高表示基因越具有生物学意义，GS= 0，说明这个基因不参与所研究的生物学问题，将在各个样本中筛选出的表达有差异的基因的表达数据构建权重共表达网络，选择软阈值为9，筛选共表达模块；

4）共表达模块的富集分析：为观察各个共表达模块的功能，使用R软件包clusterProfiler对各个模块进行KEGG通路富集分析。

4.如权利要求2所述的利用权重基因共表达网络挖掘结直肠癌放疗特异性响应基因的方法，其特征在于，所述步骤2具体包括以下步骤：

1）各模块与癌症样本相关性分析：根据放疗样本类型，将对应属性的样本定义为1，其他样本定义为0，所有样本定义值形成一个表型矩阵，计算每个模块与表型矩阵中各个因素的相关性，筛选出与癌症放疗相关的模块；

2）癌症相关的模块的权重共表达网络构建：根据各个共表达模块中的基因的表达关系，选择共表达权重大于0.2的共表达模块对作为最终的共表达网络的边，导出癌症相关的模块的共表达网络；

3）共表达-互作调控网络构建：参与相关功能的基因在蛋白互作网络中通常为距离较近的节点，即在蛋白互作网络中表现出临近关系，为了识别出癌症预后显著相关的重要基因，从ps://www.ncbi.nlm.nih.gov/pmc/articles/PMC5210659/ 下载HIPPIE v2.0，将得到的癌症相关模块的共表达网络映射到人类蛋白互作网络中，构建共表达-互作特异性网络；

4）共表达-互作特异性网络度分布分析：生物学网络是复杂网络，往往具有自组织、自相似、吸引子、小世界、无尺度等特征，对共表达-互作特异性网络进行分析，观察共表达-互作特异性网络是否符合无尺度特征；

5）富集共表达基因

共表达基因往往和其周围临近基因共同参与相同的疾病通路或生物学过程，核心蛋白和其临近蛋白在表达量上也具有显著相关性，因此通过统计每个基因周围共表达的基因个数和连接的所有基因个数，以整个人类蛋白互作网络作为统计背景，构建统计学模型，对每个基因做Fisher富集检验，筛选显著富集FDR<0.01且共表达且互作的基因占所有与其互作的基因比例大于0.1的网络节点，最终得到目的基因（癌症放疗响应特异性基因）；

6）癌症放疗响应特异性基因的表达谱分析：为观察这些基因之间的表达关系，运用表达谱层次聚类分析的方法。

5.如权利要求2所述的利用权重基因共表达网络挖掘结直肠癌放疗特异性响应基因的方法，其特征在于，所述步骤3具体包括以下步骤：

1）癌症放疗响应特异性基因的富集分析：为观察这些基因的功能，对这些基因做KEGGPathway富集分析；

2）癌症放疗响应特异性基因在外部数据集中的表达关系：通过gepia在线工具分析癌症特异性基因在TCGA数据集、GETx数据集中的表达关系。

Images