生物医学大数据挖掘及生物信息学案例分析系列课程第三讲 表观遗传学数据分析

第三讲

表观遗传学数据分析PartI表观遗传学概念PartIIDNA甲基化数据分析PartIII肿瘤甲基化数据库网络工具介绍PartIV案例：基于数据挖掘的全基因组DNA甲基化数据分析PartV拷贝数变异数据分析的软件实现-Genovar

PartI表观遗传学概念表观遗传学是研究不涉及DNA序列改变的情况下，DNA甲基化谱、染色质结构状态和基因表达谱在细胞代间传递的遗传现象的一门科学。计算表观遗传学是应用及开发生物信息学方法（统计分析，模式识别等）解决生物医学相关的表观遗传学问题。生物信息学构架了基因组学与表观基因组学的桥梁计算表观遗传学表观遗传学数据分类肿瘤神经退行性疾病心血管疾病精神性疾病代谢性疾病表观遗传学在医学中的应用DNA甲基化和组蛋白修饰有潜在的临床用途用于普遍临床实践抑癌基因高甲基化谱可用于癌症病人预后特定基因的高甲基化可对治疗反应进行预测诊断工具患者预后治疗反应预测PartII基因组的DNA甲基化数据分析一、CpG岛的DNA甲基化调控基因表达（一）DNA甲基化与CpG岛DNA甲基化是一种发生在DNA序列上的化学修饰，可以在转录及细胞分裂前后被稳定地遗传。DNA甲基化是重要的表观遗传代码。CpG是胞嘧啶（C）—磷酸（p）—鸟嘌呤（G）的缩写。CpG岛主要位于基因的启动子和外显子区域，富含CpG，长度为300—3000bp。CpG岛甲基化和转录的关系DNA甲基化的发生机制在DNA甲基化转移酶（DNMTs）的作用下使CpG二核苷酸5’-端的胞嘧啶转变为5’-甲基胞嘧啶。DNA甲基化通常抑制基因表达，去甲基化则诱导了基因的重新活化和表达。DNA甲基化转移酶（二）DNA甲基化对转录的调控1.DNA甲基化阻碍转录因子的结合2.DNA甲基化识别染色质标记3.DNA甲基化募集其他蛋白引起染色质沉默4.DNA甲基化影响核小体定位二、基因组CpG岛的预测算法1）实验方法确定的基因组范围CpG岛图谱②基于相邻CpG二核苷酸距离的CpG岛预测算法计算基因组范围内相邻CpG二核苷酸之间的距离，估计距离的理论分布并确定阈值，根据域值确定CpG簇，通过合并CpG簇，得到CpG岛。2）CpG岛的计算预测方法通过设定窗口宽度大小，判断窗口内的序列片段是否满足CpG岛的定义（如长度、GC含量等）。①基于窗口滑动法的CpG岛预测算法结合DNA结构、组蛋白修饰、DNA甲基化状态以及转录因子结合等信息定位CpG岛等方法，可提高预测的准确性。UCSC数据库的截图展示了三个CpG岛三、DNA甲基化的检测方法目前常用的DNA甲基化检测方法是将待检序列中甲基化的胞嘧啶转化为其他碱基组成的变化。1.限制性内切酶法2.重亚硫酸钠法3.高通量测序的DNA甲基化检测方法(Microarray）不同癌症之间基因的超甲基化水平存在差异TCGA下载的DNA甲基化微阵列的数据格式甲基化矩阵数据注释数据甲基化矩阵中的DNA甲基化水平用beta值来表示，beta表示的是甲基化等位基因密度与非甲基化等位基因密度的比值，该值为[0,1]之间的连续值，0表示未甲基化，1表示完全甲基化。DNA甲基化数据矩阵分析和基因表达矩阵分析类似差异甲基化区域筛选聚类火山图亚型分析富集分析数据挖掘，特征分析……案例：应用R的CancerSubtypes软件包实现DNA甲基化数据分析（一）准备2个数据文件60个Tumor样本，8个normal样本dnamethy.csvdnamethynohead.csvread.table("d:\\dnamethy.csv",header=TRUE,sep=",")->aread.table("d:\\dnamethynohead.csv",header=FALSE,sep=",")->datadata<-as.matrix(data)rownames(data)<-a[,1]colnames(data)<-c(rep("Tumor",60),c(rep("Normal",8)))（二）将整理好的excel数据导入library(CancerSubtypes)（三）检查数据的分布data.checkDistribution(data)(四）缺失值填补data=data.imputation(data,fun="median")（五）数据正态化data=data.normalization(data,type="feature_Median",log2=FALSE)类型：feature_Median中位正态化feature_Mean均值正态化feature_zscore按基因标化sample_zscore按样本标化log2如果是TRUE,做以2为底的对数转换（常用于RNAseq数据）（六）绘制双向聚类图library(gplots)heatmap.2(data,col=redgreen)（七）求差异表达基因T_data<-data[,1:60]N_data<-data[,61:68]result=DiffExp.limma(Tumor_Data=T_data,Normal_Data=N_data,topk=NULL,RNAseq=FALSE)RNAseq=TRUE表示用RNAseq表达，RNAseq=FALSE表示microarrray数据。write.csv(result,file="d:\\result.csv")差异甲基化区域：CpG的编号对应原数据的顺序编号。（八）用输出结果绘制火山图

数据格式如下：library(ggplot2)read.table("d:\\volcano.csv",header=TRUE,sep=",")->datathreshold<-as.factor(abs(data$logFC)>=1.5&data$Pvalue<0.05)ggplot(data,aes(x=logFC,y=-log10(Pvalue),colour=threshold))+geom_point(alpha=0.4,size=1.75)+geom_vline(xintercept=c(-1,1),lty=4,col="grey",lwd=0.5)+geom_hline(yintercept=-log10(0.05),lty=4,col="grey",lwd=0.5)红色字体表示设置散点和分界线的颜色和样式。火山图对提取的差异表达基因，也可以应用主成分分析绘制分类图转置存为cpgpca.csv，存于d盘下。library(ggfortify)read.table("d:\\cpgpca.csv",header=TRUE,sep=",")->aautoplot(stats::prcomp(a[-1]),data=a,frame=TRUE,frame.type='t',frame.colour='group',colour='group')肿瘤样本正常样本说明差异的DNA甲基化区域(基因）能够将样本分离开（九）肿瘤亚型和正常的比较1)应用甲基化数据将肿瘤分为三个亚型。result1=ExecuteCC(clusterNum=3,d=T_data,maxK=10,clusterAlg="hc",distance="pearson")group=result1$group2)每一个亚型和正常数据的差异甲基化diff=DiffExp.limma(Tumor_Data=T_data,Normal_Data=N_data,group=group,topk=NULL,RNAseq=FALSE)write.csv(diff,file="d:\\subtype.csv")每一种亚型和正常比较的DNA差异甲基化区域3)绘制热图distanceMatrix=result1$distanceMatrixsilhouette=silhouette_SimilarityMatrix(group,distanceMatrix)drawHeatmap(T_data,group,silhouette=silhouette,scale="max_min",Title="DNAmethylation")sil=silhouette_SimilarityMatrix(result1$group,result1$distanceMatrix)plot(sil,col=c("red","green","blue"))4)绘制Silhouette图5)亚型之间的比较sigclustTest(T_data,group,nsim=500,nrep=1,icovest=1)PartIII肿瘤甲基化数据库肿瘤相关的DNA甲基化数据库-MethHC网络工具简介MethHC网络工具（http://MethHC..tw）整合了TCGA（TheCancerGenomeAtlas）数据库中的DNA甲基化数据、基因表达数据、miRNA表达数据以及DNA甲基化和基因表达的相关性。数据库中包含了18种人类癌症6000个样本和6548个微阵列及12,567个RNA测序数据。甲基化水平所涉及的基因区域包含promoter、enhancer、TSS1500、TSS200、5'UTR、firstexon、genebody和3'UTR、CpGislands、shelves和shores1）浏览高甲（低）甲基化基因MethHC网络工具可以帮助获得肿瘤相关基因的DNA甲基化水平。例如选择膀胱癌（BladderUrothelialCarcinoma,BLCA），选择基因区域为CpGisland，DNA甲基化水平选择平均水平（Average），则获得膀胱癌相关的250个高甲基化基因点击“mostdifferentially250”可以获得250个肿瘤样本与正常样本最显著的差异甲基化基因。比如对于第一个基因DGKZ，点击“detail”，图中结果表示的是DNA甲基化水平与基因表达值的相关性。我们也可以打开相应的数据进行查看，点击“Downloaddata”进行下载：2）基于基因搜索的DNA甲基化水平分析在主页菜单中点击“GeneSearch”，并勾选膀胱癌（blca：BladderUrothelialCarcinoma）。在弹出的界面中，基因区域选择“Promoter”，甲基化水平评价方法选择“Average”，在基因列表的空白框内输入两个基因KRT4和KRT6C，用分号分隔开：点击“Search”按钮后，输出了KRT4和KRT6C在膀胱癌和正常样本中的DNA甲基化水平的箱式图。如果还想获得KRT6C在正常样本和膀胱癌样本中的DNA甲基化水平（beta-value）的分布，点击“Downloadsummary”，可以获得两组样本DNA甲基化水平的描述信息，如p值、肿瘤样本数、正常样本数、最小值、中位数、第一和第三四分位数以及最大值等，如图所示。3）DNA甲基化水平的统计学描述

PartIV案例：突眼症的DNA全基因组甲基化数据分析案例：突眼症的DNA甲基化数据处理1）识别差异甲基化区域2）单变量的Logistic回归分析批量作Logistic回归read.table("d:\\logit.csv",header=TRUE,sep=",")->aresult<-array(0,dim=c(39,4))for(iin1:39){logr<-glm(group~sex+age+duration+a[,i+4],data=a,family="binomial")result[i,1]<-exp(summary(logr)$coefficients[5,1])result[i,2]<-exp(summary(logr)$coefficients[5,1]-1.96*summary(logr)$coefficients[5,2])result[i,3]<-exp(summary(logr)$coefficients[5,1]+1.96*summary(logr)$coefficients[5,2])result[i,4]<-summary(logr)$coefficients[5,4]}result（设置4列：分别放入OR值，OR值置信区间的上下限，P值）一次Logistic回归的结果运行结果如下：4列分别为第1列：OR值第2列：OR值置信区间的下限第3列：OR值置信区间的上限第4列：P值3）聚类分析CancerSubtypes软件包4）主成分分析5）GO功能富集分析6）与表型的线性相关分析

PartV拷贝数变异数据分析的软件实现DNA拷贝数变异（CopyNumberVariants,CNV）是由基因组发生重排而导致的，一般指长度1kb以上的基因组大片段的拷贝数增加或者减少，主要表现为亚显微水平的缺失和重复：缺失1或0拷贝CNVloss正常2拷贝重复>2拷贝CNVgain

当前，对于全基因组范围内的CNV检测主要有三种技术：基于芯片的的比较基因组杂交技术（array-basedComparativeGenomicHybridization,aCGH）、SNP分型芯片技术和新一代测序技术。aCGH技术具有高敏感度和准确度的特点，探针密度越高所得结果越准确。相对于aCGH技术而言，SNP芯片仅使用单杂交就可以完成。

当前的新一代测序技术就具有很多优点了，如该技术不需要知道更多的背景知识和设计工作，应用配对测序就可以鉴定出复杂的结构变化，而且费用低廉。DNA拷贝数变异数据的分析软件-Genovar

Genovar是一款基于Java开发的软件，可从/上下载。该软件不仅可探查CNV区域，还可以将探查到的CNV区域和基因组变异数据库DGV（http://projects.tcag.ca/variation）以及dbSNP数据库作比较。同时也可以对aCGH或序列格式的基因组数据进行可视化。该软件采用的是Smith-WatermanArray（SW-ARRAY）算法，算法提供了探查CNV的动态程序。由于该算法依赖于参数的域值，因此结果对域值的变化会更加敏感。下面来看具体的操作步骤。Genovar界面如图：该软件可以打开四种格式的数据文件，分别为：aCGH格式的文件、基因位置文件、BAM格式的文件和UCSCfasta格式的文件。下面打开一个aCGH格式的数据，该数据中包括正常组织样本7例，肿瘤组织样本28例。主要信息包括探针名称、染色体信息、起始及终止位置信息及样本的ratio值。在软件中将文件打开后，就出现图所示的界面了，该界面将第一个样本的全基因组aCGH数据进行可视化，其中绿色表示探针的ratio为正值，红色表示探针的ratio为负值。从Sample对应的下拉菜单中选择任意一个样本，就可以按照该样本的的数值进行基因组绘图了，如图所示。如果想看某一个肿瘤样本T2_1在6号染色体的分布情况，在Sample下拉菜单中选择“T2_1”，在CHR下拉菜单中选择“Chr6”，就可以获得如图所示的结果，其中左侧是数据表，右侧为绘制出的肿瘤样本T2_1的探针在6号染色体上的分布情况。下面对于T2_1样本，进行CNV区域的探查，如果想探查其他样本，只需在Sample下拉菜单中选择相应的样本就可以了，如图所示。此时弹出的对话框如图所示。软件中提供了两种方法：域值法ThresholdBasedMethod和Smith-Waterman