--- a
+++ b/9-25CpG clusterProfiler.R
@@ -0,0 +1,449 @@
+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)
+}