特異性分析

#設(shè)定工作路徑
setwd("D:/DMP_mk_")
#加載需要的R包
#rm(list=ls())
library(reshape2)
library(dplyr)

#加載需要的數(shù)據(jù)
rm(list = ls())
#單個甲基化位點
meth <- read.table("meth_1.txt")#改動的地方
exprset <- read.table("exprset.txt")
#啟動子區(qū)域平均甲基化
meth_prom <- read.table("meth_prom_250updown.txt")#改動的地方
meth_spm <- read.table("meth_1.txt.spm") #改動的地方
#表達矩陣
exprset_spm <- read.table("exprset.txt.spm")
meth_prom_spm <- read.table("meth_prom_250updown.txt.spm")#改動的地方
#甲基化位點注釋信息
load(file = "anno_1.Rdata")#改動的地方

#生成已經(jīng)添加行名和列名的長數(shù)據(jù)
meth_prom_L <- cbind(rep(rownames(meth_prom),ncol(meth_prom)),melt(meth_prom))
colnames(meth_prom_L) <- c("gene","tissue","rate")
meth_prom_spm_L <- cbind(rep(rownames(meth_prom),ncol(meth_prom_spm)),melt(meth_prom_spm))
colnames(meth_prom_spm_L) <- c("gene","tissue","spm")
exprset_L <- cbind(rep(rownames(exprset),ncol(exprset)),melt(exprset))
colnames(exprset_L) <- c("gene","tissue","value")
exprset_spm_L <- cbind(rep(rownames(exprset),ncol(exprset_spm)),melt(exprset_spm))
colnames(exprset_spm_L) <- c("gene","tissue","spm")

meth_L <- melt(meth)
meth_spm_L <- melt(meth_spm)

#啟動子甲基化與基因表達的關(guān)系
vs <- read.csv(file="vs.csv")
table(meth_prom_L[,1] %in% vs[,1])
gene_symbol <- vs[match(meth_prom_L[,1],vs[,1]),2]
tmp <- transform(meth_prom_L,gene = gene_symbol)
merge_1 <- merge(exprset_L,tmp,by = c("gene","tissue"))
length(merge_1)
png(file="correlation.png", bg="transparent")
plot(merge_1$rate,merge_1$value)
dev.off()
cor.test(merge_1$value, merge_1$rate,method = c("pearson"))

#組織特異性甲基化區(qū)域與組織特異性表達的關(guān)系
meth_prom_plus <- merge(meth_prom_L,meth_prom_spm_L,by = c("gene","tissue"))
exprset_plus <- merge(exprset_L,exprset_spm_L,by = c("gene","tissue"))
methpromspec <- filter(meth_prom_plus,spm>0.6)
exprsetspec <- filter(exprset_plus,spm>0.9)
vs <- read.csv(file="vs.csv")
table(methpromspec[,1] %in% vs[,1])
gene_symbol_1 <- vs[match(methpromspec[,1],vs[,1]),2]
tmp_1 <- transform(methpromspec,gene = gene_symbol_1)
merge_2 <- merge(tmp_1,exprsetspec,by = c("gene","tissue"))
nrow(methpromspec)
nrow(exprsetspec)
nrow(merge_2)

#組織特異性甲基化位點與組織特異性表達的關(guān)系
meth_spm_L <- filter(meth_spm_L,variable != "DPM")
meth_plus <- cbind(meth_L,meth_spm_L$value)
colnames(meth_plus) <- c("tissue","rate","spm")
meth_plus <- cbind(anno_1,meth_plus)#改動的地方
methspec <- filter(meth_plus,spm>0.9)
merge_3 <- merge(methspec,exprsetspec,by.x = "gene_symbol",by.y = "gene")
nrow(methspec)
nrow(merge_3)

得到1495個組織特異性的啟動子甲基化區(qū)域。（由12118個啟動子甲基化區(qū)域得到）
不重復(fù)的基因只有1392個，有一些基因在不只一個組織中，也有的基因在組織中重復(fù)出現(xiàn)?？赡苁嵌鄠€ucsc name對應(yīng)到了一個gene symbol上。
12059個組織特異性基因。（由36598個基因得到）
overlap數(shù)目為24個。（0.01605351）
32312個組織特異性的甲基化位點。（由715765個位點得到）
overlap數(shù)目為5455個。（ 0.1688227）

1.相關(guān)性分析結(jié)果

> cor.test(merge_1$value, merge_1$rate,method = c("pearson"))

    Pearson's product-moment correlation

data:  merge_1$value and merge_1$rate
t = -5.1149, df = 143011, p-value = 3.142e-07
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
 -0.018705797 -0.008342175
sample estimates:
        cor 
-0.01352435 

2.epigenome association analysis

load(file = "united_data_1.Rdata")#改動的地方
coordinate <- data.frame(meth_1$chr,meth_1$start,meth_1$strand)#改動的地方
coordinate <- coordinate[rep(1:nrow(meth_1),times=ncol(meth)),]
meth_plus_ <- cbind(meth_plus,coordinate)
#篩選出特異性的甲基化位點
methspec <- filter(meth_plus_,spm>0.8)
#得到染色體位置對應(yīng)的探針，做表觀關(guān)聯(lián)分析
colnames(methspec)[11] <- "chr"
colnames(methspec)[12] <- "pos"
colnames(methspec)[13] <- "strand"
library(IlluminaHumanMethylation450kanno.ilmn12.hg19)
 data(Locations)#含有染色體位置和探針信息
Locations <- data.frame(transform(Locations,"probe"=rownames(Locations)))
specific_probe <- data.frame(merge(Locations,methspec,
                                   by=c("chr","pos","strand")))
for (i in 1:13)
{
  assign(paste0(unique(specific_probe$tissue)[i]),
         unique(specific_probe[specific_probe$tissue==
                                 unique(specific_probe$tissue)[i],4]))
  filename=paste0(unique(specific_probe$tissue)[i],".txt")
  write.table(get(paste0(unique(specific_probe$tissue)[i])),
                  file = filename,quote = F,row.names = F,col.names = F)
}
#驗證是否只在一個組織中特異
library(tidyfst)
methspec %>% count_dt(chr,pos,strand) -> count_methspec

使用工具：https://bigd.big.ac.cn/ewas/toolkit

3.對組織特異性的位點進行描述性統(tǒng)計

高甲基化or低甲基化；

> dim(methspec)
[1] 32312    10
> sum(methspec$rate<20)
[1] 26438
> sum(methspec$rate>70)
[1] 162

在基因的哪個區(qū)域；

 sum(methspec$prom==1)
[1] 14252
> sum(methspec$exon==1)
[1] 11419
> sum(methspec$intron==1)
[1] 14457
> sum(methspec$CpGi==1 | methspec$shores==1)
[1] 30196
> sum(methspec$CpGi==1)
[1] 28110
> sum(methspec$shores==1)
[1] 2086

在哪些組織；

4.組織特異性基因go和KEGG富集分析

setwd("D:/genes")
#加載包
library(clusterProfiler)
library(topGO)
library(DO.db)
library(pathview)
library(DOSE)
library(org.Hs.eg.db)
library(Rgraphviz)

i=5

  #得到gene_symbol
  unique(exprsetspec$tissue)[i]
  assign(paste0("x",i),unique(exprsetspec[exprsetspec$tissue==
                                         unique(exprsetspec$tissue)[i],1]))
  #id轉(zhuǎn)換
  eg <- bitr( get(paste0("x",i)), fromType="SYMBOL", toType=c("ENTREZID"), 
              OrgDb="org.Hs.eg.db")
  head(eg)
  assign(paste0("genelist",i),eg[,2])
  #GO富集分析
  go <- enrichGO( get(paste0("genelist",i)), OrgDb = org.Hs.eg.db, ont='ALL',
                  pAdjustMethod = 'BH',
                  pvalueCutoff = 0.05, qvalueCutoff = 0.2,keyType = 'ENTREZID')
  #查看主要富集到哪條通路
  dim(go)
  dim(go[go$ONTOLOGY=='BP',])
  dim(go[go$ONTOLOGY=='CC',])
  dim(go[go$ONTOLOGY=='MF',])
  #簡單的可視化
  filename=paste0("GObarplot_",unique(exprsetspec$tissue)[i],".png")
  png(file=filename, bg="transparent",height = 480,width = 700)
  barplot(go,showCategory=20,drop=T)
  dev.off()
  
  filename=paste0("GOdotplot_",unique(exprsetspec$tissue)[i],".png")
  png(file=filename, bg="transparent",height = 480,width = 700)
  dotplot(go,showCategory=20,orderBy = "x")
  dev.off()
  
  
  #KEGG富集分析
  if(F)
  {
    kegg <- enrichKEGG(
      get(paste0("genelist",i)),
      organism = "hsa",
      keyType = "kegg",
      pvalueCutoff = 0.05,
      pAdjustMethod = "BH",
      minGSSize = 10,
      maxGSSize = 500,
      qvalueCutoff = 0.2,
      use_internal_data = FALSE
    )
    head(kegg)
    dim(kegg)
    
    filename=paste0("kegg_",unique(exprsetspec$tissue)[i],".png")
    png(file=filename, bg="transparent")
    dotplot(kegg, showCategory=30,orderBy = "x")
    dev.off()
    
    
    
  }
  
  #疾病通路富集分析
  DO <- enrichDO(
    get(paste0("genelist",i)),
    ont = "DO",
    pvalueCutoff = 0.05,
    pAdjustMethod = "BH",
    minGSSize = 10,
    maxGSSize = 500,
    qvalueCutoff = 0.2,
    readable = FALSE
  )
  dim(DO)
  
  filename=paste0("DO_",unique(exprsetspec$tissue)[i],".png")
  png(file= filename, bg="transparent",height = 480,width = 480,)
  dotplot(DO, showCategory=20,orderBy = "x")
  dev.off()  
  i=i+1

不知道為什么循環(huán)運行之后圖片就都顯示不出來，沒有找到具體的原因，所以只好手動。