生物醫(yī)學(xué)大數(shù)據(jù)挖掘及生物信息學(xué)案例分析系列課程第三講 表觀遺傳學(xué)數(shù)據(jù)分析

第三講

表觀遺傳學(xué)數(shù)據(jù)分析PartI表觀遺傳學(xué)概念PartIIDNA甲基化數(shù)據(jù)分析PartIII腫瘤甲基化數(shù)據(jù)庫網(wǎng)絡(luò)工具介紹PartIV案例：基于數(shù)據(jù)挖掘的全基因組DNA甲基化數(shù)據(jù)分析PartV拷貝數(shù)變異數(shù)據(jù)分析的軟件實現(xiàn)-Genovar

PartI表觀遺傳學(xué)概念表觀遺傳學(xué)是研究不涉及DNA序列改變的情況下，DNA甲基化譜、染色質(zhì)結(jié)構(gòu)狀態(tài)和基因表達譜在細胞代間傳遞的遺傳現(xiàn)象的一門科學(xué)。計算表觀遺傳學(xué)是應(yīng)用及開發(fā)生物信息學(xué)方法（統(tǒng)計分析，模式識別等）解決生物醫(yī)學(xué)相關(guān)的表觀遺傳學(xué)問題。生物信息學(xué)構(gòu)架了基因組學(xué)與表觀基因組學(xué)的橋梁計算表觀遺傳學(xué)表觀遺傳學(xué)數(shù)據(jù)分類腫瘤神經(jīng)退行性疾病心血管疾病精神性疾病代謝性疾病表觀遺傳學(xué)在醫(yī)學(xué)中的應(yīng)用DNA甲基化和組蛋白修飾有潛在的臨床用途用于普遍臨床實踐抑癌基因高甲基化譜可用于癌癥病人預(yù)后特定基因的高甲基化可對治療反應(yīng)進行預(yù)測診斷工具患者預(yù)后治療反應(yīng)預(yù)測PartII基因組的DNA甲基化數(shù)據(jù)分析一、CpG島的DNA甲基化調(diào)控基因表達（一）DNA甲基化與CpG島DNA甲基化是一種發(fā)生在DNA序列上的化學(xué)修飾，可以在轉(zhuǎn)錄及細胞分裂前后被穩(wěn)定地遺傳。DNA甲基化是重要的表觀遺傳代碼。CpG是胞嘧啶（C）—磷酸（p）—鳥嘌呤（G）的縮寫。CpG島主要位于基因的啟動子和外顯子區(qū)域，富含CpG，長度為300—3000bp。CpG島甲基化和轉(zhuǎn)錄的關(guān)系DNA甲基化的發(fā)生機制在DNA甲基化轉(zhuǎn)移酶（DNMTs）的作用下使CpG二核苷酸5’-端的胞嘧啶轉(zhuǎn)變?yōu)?’-甲基胞嘧啶。DNA甲基化通常抑制基因表達，去甲基化則誘導(dǎo)了基因的重新活化和表達。DNA甲基化轉(zhuǎn)移酶（二）DNA甲基化對轉(zhuǎn)錄的調(diào)控1.DNA甲基化阻礙轉(zhuǎn)錄因子的結(jié)合2.DNA甲基化識別染色質(zhì)標(biāo)記3.DNA甲基化募集其他蛋白引起染色質(zhì)沉默4.DNA甲基化影響核小體定位二、基因組CpG島的預(yù)測算法1）實驗方法確定的基因組范圍CpG島圖譜②基于相鄰CpG二核苷酸距離的CpG島預(yù)測算法計算基因組范圍內(nèi)相鄰CpG二核苷酸之間的距離，估計距離的理論分布并確定閾值，根據(jù)域值確定CpG簇，通過合并CpG簇，得到CpG島。2）CpG島的計算預(yù)測方法通過設(shè)定窗口寬度大小，判斷窗口內(nèi)的序列片段是否滿足CpG島的定義（如長度、GC含量等）。①基于窗口滑動法的CpG島預(yù)測算法結(jié)合DNA結(jié)構(gòu)、組蛋白修飾、DNA甲基化狀態(tài)以及轉(zhuǎn)錄因子結(jié)合等信息定位CpG島等方法，可提高預(yù)測的準(zhǔn)確性。UCSC數(shù)據(jù)庫的截圖展示了三個CpG島三、DNA甲基化的檢測方法目前常用的DNA甲基化檢測方法是將待檢序列中甲基化的胞嘧啶轉(zhuǎn)化為其他堿基組成的變化。1.限制性內(nèi)切酶法2.重亞硫酸鈉法3.高通量測序的DNA甲基化檢測方法(Microarray）不同癌癥之間基因的超甲基化水平存在差異TCGA下載的DNA甲基化微陣列的數(shù)據(jù)格式甲基化矩陣數(shù)據(jù)注釋數(shù)據(jù)甲基化矩陣中的DNA甲基化水平用beta值來表示，beta表示的是甲基化等位基因密度與非甲基化等位基因密度的比值，該值為[0,1]之間的連續(xù)值，0表示未甲基化，1表示完全甲基化。DNA甲基化數(shù)據(jù)矩陣分析和基因表達矩陣分析類似差異甲基化區(qū)域篩選聚類火山圖亞型分析富集分析數(shù)據(jù)挖掘，特征分析……案例：應(yīng)用R的CancerSubtypes軟件包實現(xiàn)DNA甲基化數(shù)據(jù)分析（一）準(zhǔn)備2個數(shù)據(jù)文件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數(shù)據(jù)導(dǎo)入library(CancerSubtypes)（三）檢查數(shù)據(jù)的分布data.checkDistribution(data)(四）缺失值填補data=data.imputation(data,fun="median")（五）數(shù)據(jù)正態(tài)化data=data.normalization(data,type="feature_Median",log2=FALSE)類型：feature_Median中位正態(tài)化feature_Mean均值正態(tài)化feature_zscore按基因標(biāo)化sample_zscore按樣本標(biāo)化log2如果是TRUE,做以2為底的對數(shù)轉(zhuǎn)換（常用于RNAseq數(shù)據(jù)）（六）繪制雙向聚類圖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數(shù)據(jù)。write.csv(result,file="d:\\result.csv")差異甲基化區(qū)域：CpG的編號對應(yīng)原數(shù)據(jù)的順序編號。（八）用輸出結(jié)果繪制火山圖

數(shù)據(jù)格式如下：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)紅色字體表示設(shè)置散點和分界線的顏色和樣式?；鹕綀D對提取的差異表達基因，也可以應(yīng)用主成分分析繪制分類圖轉(zhuǎn)置存為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甲基化區(qū)域(基因）能夠?qū)颖痉蛛x開（九）腫瘤亞型和正常的比較1)應(yīng)用甲基化數(shù)據(jù)將腫瘤分為三個亞型。result1=ExecuteCC(clusterNum=3,d=T_data,maxK=10,clusterAlg="hc",distance="pearson")group=result1$group2)每一個亞型和正常數(shù)據(jù)的差異甲基化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差異甲基化區(qū)域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腫瘤甲基化數(shù)據(jù)庫腫瘤相關(guān)的DNA甲基化數(shù)據(jù)庫-MethHC網(wǎng)絡(luò)工具簡介MethHC網(wǎng)絡(luò)工具（http://MethHC..tw）整合了TCGA（TheCancerGenomeAtlas）數(shù)據(jù)庫中的DNA甲基化數(shù)據(jù)、基因表達數(shù)據(jù)、miRNA表達數(shù)據(jù)以及DNA甲基化和基因表達的相關(guān)性。數(shù)據(jù)庫中包含了18種人類癌癥6000個樣本和6548個微陣列及12,567個RNA測序數(shù)據(jù)。甲基化水平所涉及的基因區(qū)域包含promoter、enhancer、TSS1500、TSS200、5'UTR、firstexon、genebody和3'UTR、CpGislands、shelves和shores1）瀏覽高甲（低）甲基化基因MethHC網(wǎng)絡(luò)工具可以幫助獲得腫瘤相關(guān)基因的DNA甲基化水平。例如選擇膀胱癌（BladderUrothelialCarcinoma,BLCA），選擇基因區(qū)域為CpGisland，DNA甲基化水平選擇平均水平（Average），則獲得膀胱癌相關(guān)的250個高甲基化基因點擊“mostdifferentially250”可以獲得250個腫瘤樣本與正常樣本最顯著的差異甲基化基因。比如對于第一個基因DGKZ，點擊“detail”，圖中結(jié)果表示的是DNA甲基化水平與基因表達值的相關(guān)性。我們也可以打開相應(yīng)的數(shù)據(jù)進行查看，點擊“Downloaddata”進行下載：2）基于基因搜索的DNA甲基化水平分析在主頁菜單中點擊“GeneSearch”，并勾選膀胱癌（blca：BladderUrothelialCarcinoma）。在彈出的界面中，基因區(qū)域選擇“Promoter”，甲基化水平評價方法選擇“Average”，在基因列表的空白框內(nèi)輸入兩個基因KRT4和KRT6C，用分號分隔開：點擊“Search”按鈕后，輸出了KRT4和KRT6C在膀胱癌和正常樣本中的DNA甲基化水平的箱式圖。如果還想獲得KRT6C在正常樣本和膀胱癌樣本中的DNA甲基化水平（beta-value）的分布，點擊“Downloadsummary”，可以獲得兩組樣本DNA甲基化水平的描述信息，如p值、腫瘤樣本數(shù)、正常樣本數(shù)、最小值、中位數(shù)、第一和第三四分位數(shù)以及最大值等，如圖所示。3）DNA甲基化水平的統(tǒng)計學(xué)描述

PartIV案例：突眼癥的DNA全基因組甲基化數(shù)據(jù)分析案例：突眼癥的DNA甲基化數(shù)據(jù)處理1）識別差異甲基化區(qū)域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（設(shè)置4列：分別放入OR值，OR值置信區(qū)間的上下限，P值）一次Logistic回歸的結(jié)果運行結(jié)果如下：4列分別為第1列：OR值第2列：OR值置信區(qū)間的下限第3列：OR值置信區(qū)間的上限第4列：P值3）聚類分析CancerSubtypes軟件包4）主成分分析5）GO功能富集分析6）與表型的線性相關(guān)分析

PartV拷貝數(shù)變異數(shù)據(jù)分析的軟件實現(xiàn)DNA拷貝數(shù)變異（CopyNumberVariants,CNV）是由基因組發(fā)生重排而導(dǎo)致的，一般指長度1kb以上的基因組大片段的拷貝數(shù)增加或者減少，主要表現(xiàn)為亞顯微水平的缺失和重復(fù)：缺失1或0拷貝CNVloss正常2拷貝重復(fù)>2拷貝CNVgain

當(dāng)前，對于全基因組范圍內(nèi)的CNV檢測主要有三種技術(shù)：基于芯片的的比較基因組雜交技術(shù)（array-basedComparativeGenomicHybridization,aCGH）、SNP分型芯片技術(shù)和新一代測序技術(shù)。aCGH技術(shù)具有高敏感度和準(zhǔn)確度的特點，探針密度越高所得結(jié)果越準(zhǔn)確。相對于aCGH技術(shù)而言，SNP芯片僅使用單雜交就可以完成。

當(dāng)前的新一代測序技術(shù)就具有很多優(yōu)點了，如該技術(shù)不需要知道更多的背景知識和設(shè)計工作，應(yīng)用配對測序就可以鑒定出復(fù)雜的結(jié)構(gòu)變化，而且費用低廉。DNA拷貝數(shù)變異數(shù)據(jù)的分析軟件-Genovar

Genovar是一款基于Java開發(fā)的軟件，可從/上下載。該軟件不僅可探查CNV區(qū)域，還可以將探查到的CNV區(qū)域和基因組變異數(shù)據(jù)庫DGV（http://projects.tcag.ca/variation）以及dbSNP數(shù)據(jù)庫作比較。同時也可以對aCGH或序列格式的基因組數(shù)據(jù)進行可視化。該軟件采用的是Smith-WatermanArray（SW-ARRAY）算法，算法提供了探查CNV的動態(tài)程序。由于該算法依賴于參數(shù)的域值，因此結(jié)果對域值的變化會更加敏感。下面來看具體的操作步驟。Genovar界面如圖：該軟件可以打開四種格式的數(shù)據(jù)文件，分別為：aCGH格式的文件、基因位置文件、BAM格式的文件和UCSCfasta格式的文件。下面打開一個aCGH格式的數(shù)據(jù)，該數(shù)據(jù)中包括正常組織樣本7例，腫瘤組織樣本28例。主要信息包括探針名稱、染色體信息、起始及終止位置信息及樣本的ratio值。在軟件中將文件打開后，就出現(xiàn)圖所示的界面了，該界面將第一個樣本的全基因組aCGH數(shù)據(jù)進行可視化，其中綠色表示探針的ratio為正值，紅色表示探針的ratio為負值。從Sample對應(yīng)的下拉菜單中選擇任意一個樣本，就可以按照該樣本的的數(shù)值進行基因組繪圖了，如圖所示。如果想看某一個腫瘤樣本T2_1在6號染色體的分布情況，在Sample下拉菜單中選擇“T2_1”，在CHR下拉菜單中選擇“Chr6”，就可以獲得如圖所示的結(jié)果，其中左側(cè)是數(shù)據(jù)表，右側(cè)為繪制出的腫瘤樣本T2_1的探針在6號染色體上的分布情況。下面對于T2_1樣本，進行CNV區(qū)域的探查，如果想探查其他樣本，只需在Sample下拉菜單中選擇相應(yīng)的樣本就可以了，如圖所示。此時彈出的對話框如圖所示。軟件中提供了兩種方法：域值法ThresholdBasedMethod和Smith-Waterman