黄丽丽，杨 昆

(杭州电子科技大学计算机学院，浙江 杭州 310018)

年龄有关的差异甲基化且差异表达基因的特点

黄丽丽，杨 昆

(杭州电子科技大学计算机学院，浙江 杭州 310018)

当前DNA甲基化和基因表达之间的复杂关系还没有一个明确结论，年龄上差异甲基化基因的特点有待进一步研究.在年龄间距大的匹配的DNA甲基化和表达的数据上，基于区域分析方法识别年龄上差异甲基化区域，提取差异甲基化且差异表达基因，并集成多种数据鉴别它们的特点.分析结果表明，差异甲基化区域的变化模式在基因坐标区域内具有很高的一致性；Gene Body是DNA甲基化的重要调控区域，大部分差异甲基化区域分布于此区域，高甲基化区域和低甲基化区域分别更易于出现在CGI promoter基因和nonCGI promoter基因；在差异甲基化区域与基因表达水平的关系上，非差异表达/上调/下调的差异甲基化区域都同时存在正关联、负关联和无关联，在上述的复杂关系中，非差异表达的DMR偏向正相关，反之差异表达的DMR偏向负相关；对正关联和负关联的差异甲基化且差异表达基因利用蛋白质互作网络信息进行分类，同时对这两类基因进行GO富集性分析，结果显示，正关联和负关联的差异甲基化且差异表达基因与蛋白质互作的不同功能模块有关联性，提示两类基因很可能关联于不同的表达调控模式.

差异甲基化区域；差异表达；年龄；分类

0 引 言

DNA甲基化是一种重要的表观遗传修饰，具有重要的生物学功能，涉及众多的复杂生物过程，包括细胞分化和发育、X染色体失活、基因印记等[1].文献[2]对乳腺癌细胞的基因表达、CpG岛甲基化等数据进行了综合分析，发现甲基化的改变对基因表达的改变有直接的关联关系.文献[3]在T细胞介导的免疫应答和分化的关联基因上，确定了甲基化和表达水平之间的强负相关性.文献[4]发现CpG-rich区域的甲基化负关联于转录，CpG-poor区域的甲基化与基因表达既有正关联又有负关联.然而也有研究认为DNA甲基化对基因表达的影响作用很小.对转录开始位点200 bp区域和第一外显子的甲基化的分析，文献[5]发现大多数与年龄相关的DNA甲基化的变化只对顺式基因表达仅有一个很小的影响，主要表现在维持已有的基本表达水平.综上所述，DNA甲基化的变化和基因表达之间的复杂关系还没有一个明确结论，差异甲基化基因的特征还需要进一步研究.

本文利用DNA甲基化和基因表达的同源匹配数据来识别年龄有关的差异甲基化区域(Differentially Methylated Region，DMR)，并提取差异甲基化且差异表达基因，集成多种数据鉴别其特点.随着个体年龄的增长，个体免疫功能也会随之下降.与年龄有关的差异甲基化基因特点以及DNA甲基化改变对基因表达变化的影响的研究，对探索DNA甲基化功能以及人类与年龄相关的疾病研究具有重要意义.

1 实验数据与方法

1.1 实验数据的收集

本文实验数据来自文献[6]，包含140个志愿者外周血单核细胞的DNA甲基化数据和匹配的基因表达mRNA数据，分成2个样本集：1)年轻组21个样本，年龄范围19～30岁；2)老年组119个样本，年龄为90岁.具体实验数据通过基因表达综合数据库(GEO)下载，甲基化数据为Infinium 450 K芯片数据，芯片注释数据来自GEO，其他注释数据来自美国圣克鲁兹分校(The University of California Santa Cruz，UCSC)创建并维护的UCSC Genome Browser数据库.使用Official Symbol作为基因标识，本文提取并集成HPRD和BIOGRID这2个数据库中的Homo sapiens的相互作用数据，最终得到一个集成的蛋白质互作网络(PPIN).

1.2 差异表达基因的识别

由于同一基因有多个探针位点，因此计算多个探针位点的平均表达值作为该基因的表达值.以老年组为目标组，使用limma包计算相应参数，以错误发现率FDR<0.05为阈值识别差异表达基因(Differentially Expressed Gene，DEG)和非差异表达基因(Non-Differentially Expressed Gene，NoDEG).根据差异倍数(logFC)的正负，将差异表达基因细分为2个子类差异上调(Up)和差异下调(Down).

1.3 差异甲基化基因的识别

基于区域方法[7]，本文识别差异甲基化区域DMR，以基因为标识合并区域得到差异甲基化基因(Differentially Methylated Gene，DMG).具体处理流程：1)去除非唯一性的杂交探针；2)基于区域方法计算获得错误发现率FDR，β值等，其中，β值为甲基化程度的变化方向.以阈值FDR<0.05筛选得到差异甲基化区域，将β>0的区域记为高甲基化区域(Hyper)，β<0的区域记为低甲基化区域(Hypo)；3)删除冲突DMR(即属于一个基因的多个DMR其甲基化的变化方向不一致)，合并相同基因下的DMRs得到差异甲基化基因DMG，并细分为高甲基化基因(Hyper)和低甲基化基因(Hypo).

1.4 基因类别的划分

本文从多个视角把基因划分成不同的类别：1)根据基因表达角度分成非差异表达基因(NoDEG)、差异上调(Up)和差异下调(Down)；2)根据甲基化角度分成非差异甲基化基因(NoDMG)、高甲基化基因(Hyper)和低甲基化基因(Hypo)；3)根据启动子是否包含CpG岛(CGI)分成CGI promoter基因和nonCGI promoter基因；4)根据甲基化和基因表达水平之间的相关系数指标分成正关联(R+)和负关联(R-)基因.

1.5 网络拓扑指标的构建和分类

在集成的蛋白质互作网络中，本文计算节点的11种拓扑指标[8]：1)度中心性(DC)；2)聚类系数中心(CLUSC)；3)K-core；4)介数中心性(BETWC)；5)紧密度中心性(CLOSEC)；6)特征向量中心性(EC)；7)节点“重要性”PageRank；8)关于4个基因子类(Hyper_Down类、Hyper_Up类、Hypo_Down类和Hypo_Up类)的1N index[9]，并以它们构造分类特征，对正关联和负关联的差异甲基化且差异表达基因进行分类，鉴别它们的特点.在具体的分类中，采用支持向量机(Support Vector Machine，SVM)分类算法为训练分类器，在实际的分类计算中使用Matlab环境下的nSVMs[10]工具箱.

1.6 基因本体的分析

基因本体(Gene Ontology，GO)分为细胞组成、生物过程和分子功能3部分.为了对目标基因进行功能分析，本文以Infinium 450K芯片中所有的基因作为背景基因，使用GOrilla对所鉴别出的年龄相关的差异甲基化和差异表达基因(正关联和负关联)进行GO富集分析，其中以参数0.000 1为显著性过滤阈值.

2 实验结果与分析

2.1 差异甲基化和差异表达基因

本文的DNA甲基化和基因表达数据是匹配数据，以年轻组为对照组并以老年组为目标组，共识别出5 017个差异表达基因DEG，上调和下调分别标记为Up和Down，非差异表达基因标记为NoDEG类；识别出476个差异甲基化区域DMR，根据甲基化变化方向分别标记为Hyper(高甲基化)和Hypo(低甲基化)，排除其中的2个冲突DMR(即属于一个基因的多个DMR其甲基化变化方向不一致)，合并474个一致DMR得到450个差异甲基化基因DMG.分析发现20个基因包含多个DMR，并且仅有5%(1/20)的基因涉及冲突DMR.此结果显示单个基因的多个DMR同时高甲基化Hyper或同时低甲基化Hypo，表明DNA甲基化的变化模式在基因的坐标区域内一致性.

2.2 基因组注释特征的分析

近期对小鼠背根神经节的感觉神经元细胞的研究指出，基因启动子是否包含CpG岛(CGI)，可分为CGI promoter基因和nonCGI promoter基因，其甲基化对基因表达有不同的影响[11].根据UCSC的knownCanonical注释数据，计算得到444个差异甲基化基因(含468个DMR)TSS上游的1 kb启动子信息，结合注释数据cpgIslandExt计算CGI promoter基因和nonCGI promoter基因，结果如表1所示.费舍尔精确检验(Fisher’s exact test)结果显示高/低甲基化和基因启动子类别之间无关联的概率P值为9.2E-13，表明高甲基化Hyper和低甲基化Hypo分别更易于出现在CGI promoter基因和nonCGI promoter基因上.

表1 差异甲基化类型和基因启动子类型的列联表

为了鉴别的年龄相关的差异甲基化区域在基因不同区域上的分布规律，根据Infinium 450 k芯片的注释数据计算6类基因区域上的频率，区域的定义和详细信息参见文献[12].高/低甲基化区域的分布情况如图1所示.图1中，每个DMR区域可以涉及多个基因区域类别，结果显示差异表达基因/非差异表达基因的差异甲基化区域大多数分布于Gene Body区域.表明Gene Body区域的甲基化模式中容易出现包含多个相近CpG位点的区段性甲基化状态的改变，并且部分基因呈现功能变化，基因有差异表达.

图1 差异甲基化区域在基因功能区域上的分布

2.3 结合表达的分析

图2 三类差异甲基化区域的相关系数

由于缺少部分基因的表达数据，所识别的DMG中共有426个基因(涉及450个DMR)同时包含启动子信息和表达数据的，本文以此集合为研究对象分析DNA甲基化和基因表达的关系，并基于差异甲基化区域的CpG平均甲基化水平和对应基因表达水平之间的相关系数为指标来研究甲基化和基因表达间的关联关系.差异甲基化区域的相关系数按照3个类Up/Down/NoDEG分别统计，结果如图2所示.所有类别都同时存在3种关联性：正关联、负关联和无关联，表明DNA甲基化和基因表达水平的复杂关系，有别于已有的研究结果.

将DEG和NoDEG类按照差异甲基化区域的正负相关系数做统计，结果如表2所示.

表2 差异甲基化区域的正负相关系数和是否差异表达的列联表

在表2中，NoDEG类中相关系数为正和负的DMR分别是为193和142，而DEG类中分别为54和61，表明在上述复杂关系下非差异表达的DMR偏向正相关(193/335×100%=57.61%)；相反，差异表达的DMR偏向负相关(61/115×100%=53.04%)，通过计算得到卡方独立性检验的概率P值为0.047 6.

2.4 差异甲基化且差异表达基因的分类

图2表明正相关和负相关于基因表达的差异甲基化且差异表达基因同时存在.为了鉴别2类基因的特点，本文集成HPRD和BIOGRID，计算11种拓扑指标对正关联和负关联的差异表达且差异甲基化基因进行分类，结合10-fold交叉验证，以简单猜测结果Guess(即分类准确率为简单投票给样本个数最多的类别所占的比例)为参考结果(baseline)进行比较.基于单个特征的分类和基于特征子集(排序后前X个特征)的分类准确率如图3所示.基于单个特征分类时共有5种指标的分类结果优于参考结果，指标K-core，EC，Hypo_up’1N index最好；特征集合分类时前2个特征(CLUSC和K-core)即可达到最好的分类效果.结果表明，正关联和负关联的差异甲基化且差异表达基因与蛋白质的不同功能模块的有一定关联.

图3 正关联和负关联的差异甲基化且差异表达基因的分类准确率

2.5 GO富集分析

GO富集分析的分子功能和细胞成分以及生物过程的富集项个数和重复项个数如表3所示.分子功能、细胞成分和生物过程上，2个类别只有非常少的重复项.细胞成分上负相关类的6个term全部与突触膜和突触后膜等神经系统相关；生物过程上正相关类基因的富集term偏向负调控，涉及的负调控term数是正调控的2倍.以上结果表明，2个类别基因的富集GO term间存在极大差异，两类基因很可能关联于不同的表达调控模式，涉及不同的调控机制.

表3 正相关和负相关的差异甲基化且差异表达基因的GO富集项个数

3 结束语

本文针对目前DNA甲基化和基因表达之间的复杂关系以及年龄上差异甲基化基因的特点不明确等情况，根据DNA甲基化以及与之相配的基因表达数据来识别年龄有关的差异甲基化区域和差异表达基因，进而提取出差异甲基化且差异表达基因，利用多种注释数据挖掘其特点.为理解DNA甲基化的生物学功能提供了进一步参考，同时也为理解人类与年龄相关的疾病提供相关线索.今后将使用更多/更大的匹配的数据(包括相同的血液组织和其他组织的数据)来验证本文结果，进一步探索和挖掘正负关联的差异甲基化基因的特点.

[1]凡时财,张学工.DNA甲基化的生物信息学研究进展[J].生物化学与生物物理进展,2009,36(2):143-150.

[2]SUN Z, ASMANN Y W, KALARI K R, et al. Integrated analysis of gene expression, CpG island methylation, and gene copy number in breast cancer cells by deep sequencing[J]. PloS one, 2011,6(2):e17490.

[3]TSEREL L, KOLDE R, LIMBACH M, et al. Age-related profiling of DNA methylation in CD8+ T cells reveals changes in immune response and transcriptional regulator genes[J]. Scientific reports, 2015,5:13107.

[4]NG C W, YILDIRIM F, YAP Y S, et al. Extensive changes in DNA methylation are associated with expression of mutant huntingtin[J]. Proceedings of the National Academy of Sciences, 2013,110(6):2354-2359.

[5]YUAN T, JIAO Y, DE J S, et al. An Integrative Multi-scale Analysis of the Dynamic DNA Methylation Landscape in Aging[J]. Plos Genetics, 2015,11(2):e1004996.

[6]MARTTILA S, KANANEN L, HYRYNEN S, et al. Ageing-associated changes in the human DNA methylome: genomic locations and effects on gene expression[J]. BMC Genomics, 2015,16(1):1-17.

[7]MEILYN ONG J D H. Novel region discovery method for Infinium 450K DNA methylation data reveals changes associated with aging in muscle and neuronal pathways[J]. Aging Cell, 2014,13(1):142-155.

[8]FAISAL FE, MILENKOVIC T. Dynamic networks reveal key players in aging[J]. Bioinformatics, 2014,30(12):1721-1729.

[9]XU J, LI Y. Discovering disease-genes by topological features in human protein-protein interaction network[J]. Bioinformatics, 2006,22(22):2800-2805.

[10]FUNG G, MANGASARIAN O L. Finite Newton method for Lagrangian support vector machine classification[J]. Neurocomputing, 2003,55(1/2):39-55.

[11]HU Y. Simultaneous profiling of transcriptome and DNA methylome from a single cell[J]. Genome Biology, 2016,17(1):1-11.

[12]BIBIKOVA M, BARNES B, TSAN C, et al. High density DNA methylation array with single CpG site resolution[J]. Genomics, 2011,98(4):288-295.

The Characteristics of Age-related Differentially Methylated and Differentially Expressed Genes

HUANG Lili, YANG Kun

(SchoolofComputer,HangzhouDianziUniversity,HangzhouZhejiang310018,China)

Recently the complex relationship between DNA methylation and gene expression is not yet clear, and the characteristic of age-related differentially methylated genes is still vague. Further studies are required to investigate these questions. On the matched methylomic and transcriptomic data with the big age-gap of sample, based on the region analysis we identified age-related differentially methylated regions(DMR), differentially expressed and simultaneously differentially methylated genes, and then utilized multiple data to explore their characteristics. The results showed that the variation patterns of various DMRs of gene are highly consistent in its genomic location, gene body which contains the majority of differentially methylated regions is an important regulatory region of DNA methylation, and the hyper- and hypo-methylation regions tends to CGI promoter gene and nonCGI promoter gene, respectively. For the relationship between the DMR and its corresponding gene expression, the result showed three types of relevance (i.e. positive, negative and no correlation) exist at the same time. Furthermore, there still has a tendency: the majority of DMR of non-differential expression has a positive correlation while that of differential expression has negative correlation. The results of classification using the information of protein interaction networks and GO enrichment analysis indicated that the two kinds of differential methylation and differential expression genes (i.e. positive and negative correlation) could possibly relate to different function modules of protein interaction and be associated with different transcription regulation patterns.

differentially methylated region; differential expression; aging; classification

10.13954/j.cnki.hdu.2017.03.010

2016-11-18

国家自然科学基金资助项目(60903086)

黄丽丽(1991-)女，福建南平人，硕士研究生，生物信息学．通信作者：杨昆副教授，E-mail：yangkun@hdu.edu.cn．

Q811.4

A

1001-9146(2017)03-0045-06