setwd("E:/workplace/mywork/methy/dbgap/chf/data_chf_contr/early_chf/c1_UMN/alldata")
#methylation <- read.csv("methylation_probe_info.csv",sep=",")
#save(methylation,file="methylation.Rdata")
load(file="methylation.Rdata")
load(file="E:/workplace/mywork/methy/dbgap/chf/data_chf_contr/early_chf/c1_UMN/alldata/methylation.Rdata")
setwd("E:\\workplace\\mywork\\methy\\dbgap\\chf\\data_chf_contr\\early_chf\\c1_UMN_JHU\\train_UMN_tset_JHU/1123_dataSummary")
id3 <- read.table("xgblasso_DMP.csv",sep=",",header = T)
head(id3)
id3 = as.character(id3$Feature)
venn_CpG = methylation[methylation$probe.id %in% id3[-c(2,6,17,20,21)],]
setdiff(id3[-c(2,6,17,20,21)],venn_CpG$probe.id)
#"cg06344265"
venn_CpG$probe.id = as.character(venn_CpG$probe.id)
venn_CpG$UCSC_RefGene_Name = as.character(venn_CpG$UCSC_RefGene_Name)
venn_CpG$Relation_to_UCSC_CpG_Island = as.character(venn_CpG$Relation_to_UCSC_CpG_Island)
add_cg = c("cg06344265",NA,NA,NA,NA,NA,NA,NA,NA,37,11,NA,NA,11,NA,NA,NA,NA,NA,NA,"GRIK4",NA,"TSS200",NA,"Open sea",NA,NA,NA,NA,NA,NA,NA)
venn_CpG = rbind(add_cg,venn_CpG)
write.table(venn_CpG,"venn_probe.csv",sep=",")
{
#all dmp
id3 = read.csv("new_champ/myDMP.txt",sep="\t")
id3 = rownames_to_column(id3,"probe.id")#23581
id3 <- as.character(id3$probe.id)
venn_CpG = methylation[methylation$probe.id %in% id3,]
setdiff(id3,venn_CpG$probe.id)
id3 = read.csv("new_champ/myDMP.txt",sep="\t")
id3 = rownames_to_column(id3,"probe.id")#23581
logFC_cutoff2 <- with( id3, mean( abs( logFC ) ) + 2 * sd( abs( logFC ) ) )
id3 <- id3[id3$adj.P.Val<0.05 & abs(id3$logFC) > logFC_cutoff2,]#318
venn_CpG = methylation[methylation$probe.id %in% id3$probe.id,]
setdiff(id3$probe.id,venn_CpG$probe.id)
}
#RefGene & Island
{
venn_CpG$UCSC_RefGene_Group = as.character(venn_CpG$UCSC_RefGene_Group)
venn_CpG$UCSC_RefGene_Name = as.character(venn_CpG$UCSC_RefGene_Name)
venn_CpG$UCSC_RefGene_Group = ifelse(venn_CpG$UCSC_RefGene_Group == "", "IGR",venn_CpG$UCSC_RefGene_Group)
venn_CpG$Relation_to_UCSC_CpG_Island = ifelse(venn_CpG$Relation_to_UCSC_CpG_Island == "", "Open sea",venn_CpG$Relation_to_UCSC_CpG_Island)
#UCSC_RefGene_Group
paste0(venn_CpG$UCSC_RefGene_Group,collapse=",")
A<-strsplit(as.character(venn_CpG$UCSC_RefGene_Group), ";")
tmp2 <- mapply( cbind, A )
df0 <- ldply (tmp2, data.frame)
prop.table(table(df0$X..i..))
#TSS200 IGR TSS1500 Body 1stExon 5'UTR
#0.050 0.200 0.250 0.275 0.050 0.175
#Relation_to_UCSC_CpG_Island
paste0(venn_CpG$Relation_to_UCSC_CpG_Island,collapse=",")
A<-strsplit(as.character(venn_CpG$Relation_to_UCSC_CpG_Island), ";")
tmp2 <- mapply( cbind, A )
df1 <- ldply (tmp2, data.frame)
prop.table(table(df1$X..i..))
#Open sea N_Shelf Island S_Shore N_Shore S_Shelf
#0.44 0.04 0.20 0.08 0.20 0.04
#all cpg--methylation
methylation$UCSC_RefGene_Group = as.character(methylation$UCSC_RefGene_Group)
methylation$Relation_to_UCSC_CpG_Island = as.character(methylation$Relation_to_UCSC_CpG_Island)
methylation$UCSC_RefGene_Group = ifelse(methylation$UCSC_RefGene_Group == "", "IGR",methylation$UCSC_RefGene_Group)
methylation$Relation_to_UCSC_CpG_Island = ifelse(methylation$Relation_to_UCSC_CpG_Island == "", "Open sea",methylation$Relation_to_UCSC_CpG_Island)
A<-strsplit(as.character(methylation$UCSC_RefGene_Group), ";")
tmp2 <- mapply( cbind, A )
df00 <- ldply (tmp2, data.frame)
head(df00)
prop.table(table(df00$X..i..))
#IGR Body 3'UTR TSS1500 5'UTR 1stExon TSS200
#0.13726348 0.35680909 0.03739231 0.15484313 0.13165007 0.06802518 0.11401675
A<-strsplit(as.character(methylation$Relation_to_UCSC_CpG_Island), ";")
tmp2 <- mapply( cbind, A )
df11 <- ldply (tmp2, data.frame)
head(df11)
prop.table(table(df11$X..i..))
#Open sea S_Shore N_Shelf S_Shelf N_Shore Island
#0.35281111 0.10314165 0.05029941 0.04472638 0.13191608 0.31710536
}
#The 25 CpGs associated with HFrisk model
{
library(ELMER)
nearGenes <- GetNearGenes(data = mae,
probes = id3[c(1,3:5,7:16,18:19,22:30)],
numFlankingGenes = 20)
}
#enrichKEGG
{
library(org.Hs.eg.db)
library(clusterProfiler)
library(plyr)
library(ggplot2)
library(forcats)
library(enrichplot)
library(ReactomePA)
library(DOSE)
library(org.Hs.eg.db)
library(clusterProfiler)
library(clusterProfiler)
library(topGO)
library(Rgraphviz)
library(RDAVIDWebService)
library(plyr)
library(stringr)
library(ggplot2)
#enrichKEGG
kegg_SYMBOL_hsa <- function(genes){
gene.df <- bitr(genes, fromType = "SYMBOL",
toType = c("SYMBOL", "ENTREZID"),
OrgDb = org.Hs.eg.db)
head(gene.df)
diff.kk <- enrichKEGG(gene = gene.df$ENTREZID,
organism = 'hsa',
pvalueCutoff = 0.99,
qvalueCutoff = 0.99
)
return( setReadable(diff.kk, OrgDb = org.Hs.eg.db,keyType = 'ENTREZID'))
}
paste0(venn_CpG$UCSC_RefGene_Name,collapse=",")
A<-strsplit(as.character(venn_CpG$UCSC_RefGene_Name), ";")
tmp2 <- mapply( cbind, A )
df0 <- ldply (tmp2, data.frame)
trunk_kk=kegg_SYMBOL_hsa(df0[,1])
trunk_df=trunk_kk@result
write.csv(trunk_df,file = 'enrichKEGG_result.csv')
#other
y <- mutate(trunk_df, richFactor = Count / as.numeric(sub("/\\d+", "", BgRatio)))
parse_ratio <- function(ratio) {
ratio <- sub("^\\s*", "", as.character(ratio))
ratio <- sub("\\s*$", "", ratio)
numerator <- as.numeric(sub("/\\d+$", "", ratio))
denominator <- as.numeric(sub("^\\d+/", "", ratio))
return(numerator/denominator)
}
y2 <- mutate(y, FoldEnrichment = parse_ratio(GeneRatio) / parse_ratio(BgRatio))
y3 <- mutate(y2, geneRatio2 = parse_ratio(GeneRatio)) %>%
arrange(desc(GeneRatio))
tmp <- y2[1:20,]
ggplot(tmp, showCategory = 20,aes(FoldEnrichment, fct_reorder(Description, FoldEnrichment))) +
geom_segment(aes(xend=0, yend = Description)) +
geom_point(aes(color=pvalue, size = Count)) +
scale_color_viridis_c(guide=guide_colorbar(reverse=TRUE)) +
scale_size_continuous(range=c(2, 10)) +
theme_minimal() +
xlab("Fold Enrichment") + #Rich Factor
ylab(NULL) +
ggtitle("Enriched Pathway Ontology")
write.csv(y2,file = 'enrichKEGG_result2.csv')
write.csv(y3,file = 'enrichKEGG_result3.csv')
}
#enrichGO
{
bit_gene <- function(genes){
gene.df <- bitr(genes, fromType = "SYMBOL",
toType = c("SYMBOL", "ENTREZID"),
OrgDb = org.Hs.eg.db)
return(gene.df$ENTREZID)
}
bp <- enrichGO(bit_gene( as.character(df0[,1])), ont="all",OrgDb = "org.Hs.eg.db")
bpp = bp@result
bp3 <- setReadable(bp, OrgDb = org.Hs.eg.db)
barplot(bp, split="ONTOLOGY")+ facet_grid(ONTOLOGY~., scale="free")
write.table(bp@result,"enrichGO_result.csv")
bp = as.data.frame(bp)
enrichMap(bp)
goplot(bp, showCategory = 10)
bp2 <-simplify(bp, cutoff=0.7,by="p.adjust", select_fun=min)
as.data.frame(bp2)
write.csv(bp2,file = 'enrichGO_simplify.csv')
}
#gometh: kegg and go;gsameth
{
library(limma)
library(minfi)
library(IlluminaHumanMethylation450kanno.ilmn12.hg19)
library(IlluminaHumanMethylation450kmanifest)
library(RColorBrewer)
library(missMethyl)
library(matrixStats)
library(minfiData)
setwd("E:\\workplace\\mywork\\methy\\dbgap\\chf\\data_chf_contr\\early_chf\\c1_UMN_JHU\\train_UMN_tset_JHU/1123_dataSummary")
id3 <- read.table("xgblasso_DMP.csv",sep=",",header = T)
head(id3)
id3 <- as.character(id3$Feature)
# Get all the CpG sites used in the analysis to form the background
load(file= "new_champ/CorrectedBeta.Rdata")
all <-rownames(CorrectedBeta)
# Total number of CpG sites tested
length(all)
#methylGSA
{
library(methylGSA)
library(IlluminaHumanMethylation450kanno.ilmn12.hg19)
library(IlluminaHumanMethylationEPICanno.ilm10b4.hg19)
DEG_filter = read.csv("new_champ/myDMP.txt",sep="\t")
DEG_filter = rownames_to_column(DEG_filter,"probe.id")#23581
cpg.pval <- as.data.frame(DEG_filter[DEG_filter$probe.id %in% c(id3[-c(2,6,17,20,21)]),c("probe.id","P.Value")])#514
c = paste(DEG_filter[DEG_filter$probe.id %in% c(id3[-c(2,6,17,20,21)]),"probe.id"],DEG_filter[DEG_filter$probe.id %in% c(id3[-c(2,6,17,20,21)]),"P.Value"],sep = "=",collapse = ",")
cpg.pval <- c(cg20051875=3.0845273397844e-09,cg24205914=9.48317106472548e-07,cg10083824=1.23286930624815e-06,cg03556243=4.95956762541123e-06,cg03233656=1.032434350337e-05,cg21429551=1.12835004011784e-05,cg16781992=1.33529811893581e-05,cg27401945=1.34687254097115e-05,cg08614290=3.97354136093555e-05,cg05363438=4.55923591304913e-05,cg08101977=5.07348656655275e-05,cg10556349=5.36252008917427e-05,cg13352914=9.88886388394135e-05,cg11853697=0.000116194512198325,cg00045910=0.000133863261077109,cg00495303=0.000141563110724823,cg00522231=0.000155960771094521,cg06344265=0.000258547818107834,cg21024264=0.000272396328465786,cg05481257=0.000326204702802846,cg23299445=0.000569220592642687,cg05845376=0.000637513022532287,cg17766026=0.000662632200376124,cg25755428=0.00083234067533089,cg07041999=0.000990852821362967)
head(cpg.pval)
library(org.Hs.eg.db)
res1 = methylglm(cpg.pval = cpg.pval, #Names should be CpG IDs.
minsize = 200,
maxsize = 500,
GS.type = "KEGG")#"GO", "KEGG", or "Reactome"
#2 gene sets
res6 = methylglm(cpg.pval = cpg.pval,
GS.type = "Reactome",
minsize = 100,
maxsize = 500)
#0 gene sets
res7 = methylglm(cpg.pval = cpg.pval,
GS.type = "GO",
minsize = 100,
maxsize = 500)
#449 gene sets
res2 = methylRRA(cpg.pval = cpg.pval,
method = "ORA",
minsize = 200,
maxsize = 500)
#259 gene sets
res3 = methylRRA(cpg.pval = cpg.pval,
method = "GSEA",
minsize = 200,
maxsize = 500)
#259 gene sets
res4 = methylgometh(cpg.pval = cpg.pval,
sig.cut = 0.05,
minsize = 200,
maxsize = 500)
data(CpG2Genetoy)
head(CpG2Gene)
FullAnnot = prepareAnnot(CpG2Gene)
data(GSlisttoy)
GS.list = GS.list[1:10]
res5 = methylRRA(cpg.pval = cpg.pval,
FullAnnot = FullAnnot,
method = "ORA",
GS.list = GS.list,
GS.idtype = "SYMBOL",
minsize = 100,
maxsize = 300)
#10 gene sets
}
#par(mfrow=c(1,1))
gst <- gometh(sig.cpg = id3[-c(2,6,17,20,21)], all.cpg=all, collection = "GO",plot.bias=TRUE)
GO_result = topGO(gst)
topGSA(gst, n=10)
write.table(GO_result,"gometh_toGO.csv",sep=",")
kegg <- gometh(sig.cpg = id3[-c(2,6,17,20,21)], all.cpg = all, collection = "KEGG", prior.prob=TRUE)
#kegg <- gometh(sig.cpg = set, all.cpg = all, collection = "KEGG", prior.prob=TRUE)
KEGG_result = topKEGG(kegg)
topGSA(kegg, n=10)
#gst.kegg.prom <- gometh(sig.cpg=id3[-c(2,6,17,20,21)], all.cpg=all, collection="KEGG",genomic.features = c("TSS200","TSS1500","1stExon"))
#topGSA(gst.kegg.prom, n=10)
#gst.kegg.body <- gometh(sig.cpg=id3[-c(2,6,17,20,21)], all.cpg=all, collection="KEGG", genomic.features = c("Body"))
#topGSA(gst.kegg.body, n=10)
write.table(KEGG_result,"gometh_toKEGG.csv",sep=",")
hallmark <- readRDS(url("http://bioinf.wehi.edu.au/MSigDB/v7.1/Hs.h.all.v7.1.entrez.rds"))
gsa <- gsameth(sig.cpg=id3[-c(2,6,17,20,21)], all.cpg=all, collection=hallmark)
GSA_result = topGSA(gsa)
write.table(GSA_result,"gsameth_topGSA.csv",sep=",")
}
#3.enrichPathway
{
x <- enrichPathway(gene = bit_gene(df0[,1]) ,pvalueCutoff=0.95, readable=T )
X = as.data.frame(x)
parse_ratio <- function(ratio) {
ratio <- sub("^\\s*", "", as.character(ratio))
ratio <- sub("\\s*$", "", ratio)
numerator <- as.numeric(sub("/\\d+$", "", ratio))
denominator <- as.numeric(sub("^\\d+/", "", ratio))
return(numerator/denominator)
}
y <- mutate(X, richFactor = Count / as.numeric(sub("/\\d+", "", BgRatio)))
y3 <- mutate(y, geneRatio1 = parse_ratio(GeneRatio))
y4 <- mutate(y3, geneRatio2 = parse_ratio(GeneRatio) / parse_ratio(BgRatio))
library(ggplot2)
library(forcats)
library(enrichplot)
tmp <- y4[1:20,]
ggplot(tmp, showCategory = 20,aes(geneRatio2, fct_reorder(Description, geneRatio2))) +
geom_segment(aes(xend=0, yend = Description)) +
geom_point(aes(color=pvalue, size = Count)) +
scale_color_viridis_c(guide=guide_colorbar(reverse=TRUE)) +
scale_size_continuous(range=c(2, 10)) +
theme_minimal() +
xlab("Fold Enrichment") +
ylab(NULL) +
ggtitle("Enriched Pathway Ontology")
tmp$log_pvalue <- -log10(tmp$pvalue)
ggplot(tmp, showCategory = 20,aes(log_pvalue, fct_reorder(Description, log_pvalue))) +
geom_segment(aes(xend=0, yend = Description)) +
geom_point(aes(size = Count)) +
scale_color_viridis_c(guide=guide_colorbar(reverse=TRUE)) +
scale_size_continuous(range=c(2, 10)) +
theme_minimal() +
xlab("Fold Enrichment") +
ylab(NULL) +
ggtitle("Enriched Pathway Ontology")
dotplot(x , color = "pvalue",showCategory = 20, font.size = 12, title = "")
barplot(x, showCategory=10)
emapplot(x,showCategory=20,color = "pvalue")
heatplot(x,showCategory=20)
#heatplot(x, foldChange=geneList)
cnetplot(x, categorySize="pvalue")#foldChange=geneList
upsetplot(x)
write.csv(y4,file = 'enrichPathway.csv')
}
#4.enrichDO
{
x <- enrichDO(gene = bit_gene( df0[,1]),
ont = "DO",
pvalueCutoff = 0.5,
pAdjustMethod = "BH",
universe = x@universe,
minGSSize = 5,
maxGSSize = 500,
qvalueCutoff = 0.5,
readable = FALSE)
x@result
write.csv(x@result,file = 'enrichDO.csv')
#plot
emapplot(x)
x <- setReadable(x, 'org.Hs.eg.db')
head(x)
cnetplot(x)
write.csv(x,file = 'enrichDO2.csv')
}
#5.enrichNCG
{
ncg <- enrichNCG(bit_gene(df0[,1]),pvalueCutoff = 0.95)
head(ncg,10)
dotplot(ncg, showCategory=30) + ggtitle("dotplot for GSEA")
ridgeplot(ncg)
gseaplot2(ncg, geneSetID = 1, title = ncg$Description[1])
gseaplot2(ncg, geneSetID = 1:3)
gseaplot2(ncg, geneSetID = 1:3, pvalue_table = TRUE,color = c("#E495A5", "#86B875", "#7DB0DD"), ES_geom = "dot")
write.csv(ncg,file = 'enrichNCG.csv')
}
#6. enrichDGN
{
dgn <- enrichDGN(bit_gene(df0[,1]),pvalueCutoff = 0.55)
head(dgn,20)
dgn <- data.frame(dgn)
y <- mutate(dgn, richFactor = Count / as.numeric(sub("/\\d+", "", BgRatio)))
library(ggplot2)
library(forcats)
library(enrichplot)
tmp <- y[1:20,]
ggplot(tmp, showCategory = 20,aes(richFactor, fct_reorder(Description, richFactor))) +
geom_segment(aes(xend=0, yend = Description)) +
geom_point(aes(color=pvalue, size = Count)) +
scale_color_viridis_c(guide=guide_colorbar(reverse=TRUE)) +
scale_size_continuous(range=c(2, 10)) +
theme_minimal() +
xlab("rich factor") +
ylab(NULL) +
ggtitle("Enriched Pathway Ontology")
write.csv(y,file = 'enrichDGN.csv')
edox <- setReadable(dgn, 'org.Hs.eg.db', 'ENTREZID')
cnetplot(edox)
emapplot(dgn)
}
#7.enricher
{
library(magrittr)
library(clusterProfiler)
#data(geneList, package="DOSE")
#gene <- names(geneList)[abs(geneList) > 2]
wpgmtfile <- system.file("extdata/wikipathways-20180810-gmt-Homo_sapiens.gmt", package="clusterProfiler")
#wpgmtfile <- system.file("extdata/c5.cc.v5.0.entrez.gmt", package="clusterProfiler")
wp2gene <- read.gmt(wpgmtfile)
head(wp2gene)
wp2gene <- wp2gene %>% tidyr::separate(ont, c("name","version","wpid","org"), "%")
wpid2gene <- wp2gene %>% dplyr::select(wpid, gene) #TERM2GENE
wpid2name <- wp2gene %>% dplyr::select(wpid, name) #TERM2NAME
ewp <- enricher(bit_gene(df0[,1]), TERM2GENE = wpid2gene, TERM2NAME = wpid2name)
head(ewp)
write.csv(ewp,file = 'enricher_gmt_Homo_sapiens.csv')
#======
gmtfile <- system.file("extdata", "c5.cc.v5.0.entrez.gmt", package="clusterProfiler")
c5 <- read.gmt(gmtfile)
egmt <- enricher(bit_gene(df0[,1]), TERM2GENE=c5)
head(egmt)
}
#8.msigdbr:enricher
{
library(msigdbr)
m_df <- msigdbr(species = "Homo sapiens")
head(m_df, 2) %>% as.data.frame
m_t2g <- msigdbr(species = "Homo sapiens", category = "C3") %>%
dplyr::select(gs_name, entrez_gene)
head(m_t2g)
em <- enricher(bit_gene(df0[,1]), TERM2GENE=m_t2g)
head(em)
}
#9.groupGO
{
ggo <- groupGO(gene = bit_gene(df0[,1]),
OrgDb = org.Hs.eg.db,
ont = "BP",
level = 3,
readable = TRUE)
ggo2 <- groupGO(gene = bit_gene(df0[,1]),
OrgDb = org.Hs.eg.db,
ont = "CC",
level = 3,
readable = TRUE)
ggo3 <- groupGO(gene = bit_gene(df0[,1]),
OrgDb = org.Hs.eg.db,
ont = "MF",
level = 3,
readable = TRUE)
head(ggo)
write.csv(ggo,file = 'groupGO_BP.csv')
write.csv(ggo2,file = 'groupGO_CC.csv')
write.csv(ggo3,file = 'groupGO_MF.csv')
}
#10.search:KEGG
{
library(clusterProfiler)
search_kegg_organism('Collagen', by='kegg_code')
ecoli <- search_kegg_organism('ece', by='scientific_name')
dim(ecoli)
}
#11.enrichMKEGG
{
mkk <- enrichMKEGG(gene = bit_gene(df0[,1]),organism = 'hsa')
}
#12.enrichMeSH
{
library(meshes)
x <- enrichMeSH(bit_gene(df0[,1]), MeSHDb = "MeSH.Hsa.eg.db", database='gendoo', category = 'C')
head(x)
}
Datasets
Models
