-
Notifications
You must be signed in to change notification settings - Fork 2
/
Enrichment_missMethyl_GREAT.R
89 lines (72 loc) · 2.77 KB
/
Enrichment_missMethyl_GREAT.R
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
library(missMethyl)
library(methylGSA)
library(tidyr)
library(limma)
library(ddpcr)
library(dplyr)
library(tibble)
library(openxlsx)
#missMethy
asGO<-gometh(
AsCpGs, #interested CpGs
all.cpg = allCpGs, #Background CpGs
collection = "GO",
array.type = "450K",
plot.bias = FALSE,
prior.prob = TRUE,
anno = NULL,
equiv.cpg = TRUE,
fract.counts = TRUE
)
asGO1 <- asGO %>% filter(N > 5 & N < 2000)
asGO1 <- asGO1[order(asGO1$FDR),] %>% filter(FDR< 0.05)
#----------------------------------------------------------------------------------------------------------------------
#GREAT
library(rtracklayer)
gran<-as.data.frame(hm450.manifest)
subgran<-gran[rownames(gran) %in% IPS_all_CpGs,]
subgran <- rownames_to_column(subgran)
colnames(subgran)[1]<-"name"
subgranges<- subgran %>% dplyr::select(seqnames, start, end, name)
subgranges1<-GRanges(subgranges)
export(subgranges1, "background.bed", format = "bed")
subgran<-gran[rownames(gran) %in% AsCpGs,]
subgran <- rownames_to_column(subgran)
colnames(subgran)[1]<-"name"
subgranges<- subgran %>% dplyr::select(seqnames, start, end, name)
subgranges1<-GRanges(subgranges)
export(subgranges1, "Astrocytes.bed", format = "bed")
# get tables for each GO set ordered by variable of interest
get_GO_tables <- function(results, GO_table_to_get, order_by){
tb = getEnrichmentTables(results)
# restrict to which GP processes are desired
tb <- tb[[GO_table_to_get]]
# calculate FDR from raw P-values
tb$Hyper_FDR_Q_val <- p.adjust(tb$Hyper_Raw_PValue, method = "BH")
# restrict to GO terms w/ FDR<0.01
tb <- tb[tb$Hyper_FDR_Q_val<0.01,]
# order results by FDR or enrichment value
if(order_by == "qval")
tb <- tb[order(tb$Hyper_FDR_Q_val, decreasing = FALSE),]
if(order_by == "enrichment")
tb <- tb[order(tb$Hyper_Fold_Enrichment, decreasing = TRUE),]
# return processed table
tb
}
library(rGREAT)
library(genomation)
background <- readBed('background.bed')
Ast <- readBed('Astrocytes.bed')
# submit GREAT job (specifying background set as 450K probes that passed QC)
job = submitGreatJob(gr = Ast, bg = background, species = "hg19",request_interval = 0)
# GO Biological Process
tb_bp_qval_ast <- get_GO_tables(job, "GO Biological Process", "qval")
# GO Molecular Function
tb_mf_qval_ast <- get_GO_tables(job, "GO Molecular Function", "qval")
# GO Cellular Component
tb_cc_qval_ast <- get_GO_tables(job, "GO Cellular Component", "qval")
#combined
ast_GREAT<-rbind(tb_bp_qval_ast,tb_mf_qval_ast,tb_cc_qval_ast)
ast_GREAT <- ast_GREAT[order(ast_GREAT$Hyper_FDR_Q_val), ] %>% filter(Total_Genes_Annotated > 5 & Total_Genes_Annotated < 2000)
ast_GREAT$Category<-ifelse(ast_GREAT$ID %in% tb_bp_qval_ast$ID, "BP",ifelse(ast_GREAT$ID %in% tb_mf_qval_ast$ID,"MF","CC"))
ast_GREAT_search<-dplyr::filter(ast_GREAT , grepl("ast",name))