特异性分析

#设定工作路径
setwd("D:/DMP_mk_")
#加载需要的R包
#rm(list=ls())
library(reshape2)
library(dplyr)

#加载需要的数据
rm(list = ls())
#单个甲基化位点
meth <- read.table("meth_1.txt")#改动的地方
exprset <- read.table("exprset.txt")
#启动子区域平均甲基化
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")#改动的地方

#生成已经添加行名和列名的长数据
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)

#启动子甲基化与基因表达的关系
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"))

#组织特异性甲基化区域与组织特异性表达的关系
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)

#组织特异性甲基化位点与组织特异性表达的关系
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个组织特异性的启动子甲基化区域。（由12118个启动子甲基化区域得到）
不重复的基因只有1392个，有一些基因在不只一个组织中，也有的基因在组织中重复出现。可能是多个ucsc name对应到了一个gene symbol上。
12059个组织特异性基因。（由36598个基因得到）
overlap数目为24个。（0.01605351）
32312个组织特异性的甲基化位点。（由715765个位点得到）
overlap数目为5455个。（ 0.1688227）

1.相关性分析结果

> 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)
#得到染色体位置对应的探针，做表观关联分析
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.对组织特异性的位点进行描述性统计

高甲基化or低甲基化；

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

在基因的哪个区域；

 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转换
  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

不知道为什么循环运行之后图片就都显示不出来，没有找到具体的原因，所以只好手动。