04_DMP_mFuzzClustering.Rmd

---
title: "Mfuzz Clustering"
format: html
author: "Dr Jamie Soul, Dr Euan McDonnell"
---

### Reset

```{r}

rm(list=ls(all.names=T)); gc()

```

### Load the libraries
```{r}
library(Mfuzz)
library(dplyr)
library(tidyr)
library(ggplot2)
library(cowplot)
library(RColorBrewer)
library(IlluminaHumanMethylationEPICanno.ilm10b4.hg19)
library(reshape2)
```

### Set up the directories

```{r}

# Defined main dir
scripts_dir <- getwd()
main_dir <- "methylation/"
data_dir <- "rawdata/"
res_dir <- "results/"
pca_dir <- "pca/"
dmr_dir <- "results/dmr/"
mfz_dir <- paste0(dmr_dir, "mfuzz/")
de_dir <- "differential_methylation/"
util_dir <- "utilities/"
goclust_dir <- paste0(mfz_dir, "go/")

# Create directories
sapply(c(mfz_dir, goclust_dir), dir.create)

```

### Source R scripts

```{r}

# Source script that sources files
source("install/source_rscripts.R")

# Source relevant scripts
SourceExternalScripts(paste0(util_dir, "R/"), "*.R$", ignore.case=FALSE)

```

## Parameters

```{r}

# Clustering
clust <- 3
seed_val <- 5
min_clust_membership <- 0.99
min_clust <- 2
max_clust <- 5

# DMRs
trt_log2fc_thr <- 0.1
maxgap_val <- 1000
min_dmr_size <- 1

# Workflow
dmin <- FALSE

```

## Load the dmrs and the data
```{r}

# Load in m-values
mvals <- read.delim(paste0(res_dir, "fun_methylation_matrix.mvals.tsv"), sep="\t", header=TRUE, row.names=1)

# Load in Stage CpGs
sigCpG_stage <- read.delim(paste0(de_dir, "differential_methylation.limma_results.continuous_covariates.treat_log2FC",trt_log2fc_thr,"_thr.tsv"), sep="\t", header=TRUE)
sigCpG_sex <- read.delim(paste0(de_dir, "differential_methylation.limma_results.continuous_covariates.ebayes_adjusted.tsv"), sep="\t", header=TRUE)
sigCpG_stage <- sigCpG_stage[sigCpG_stage$CONTRAST == "Z_STAGE_WEEKS",] 
sigCpG_sex <- sigCpG_sex[sigCpG_sex$CONTRAST == "sexmale - sexfemale" ,]
# sigCpG_stage <- sigCpG_stage[grepl("^cg",sigCpG_stage$Name),]
# sigCpG_sex <- sigCpG_sex[grepl("^cg",sigCpG_sex$Name),]

# Load in DMR CpGs
stage_dmr <- read.delim(paste0(dmr_dir, "dmrff_results.stage_as_continuous.treat_log2FC",trt_log2fc_thr,".cpgMvals.tsv"), sep="\t", header=TRUE)
sex_dmr <- read.delim(paste0(dmr_dir, "dmrff_results.sex.cpgMvals.tsv"), sep="\t", header=TRUE)
rownames(stage_dmr) <- stage_dmr$CpG
rownames(sex_dmr) <- sex_dmr$CpG

# Keep just sig CpGs
sigCpG_stage <- sigCpG_stage[sigCpG_stage$bonferroni  < 0.05,]
sigCpG_sex <- sigCpG_sex[sigCpG_sex$bonferroni < 0.05,]
sigCpG_stageDmr <- stage_dmr[stage_dmr$SigCpG != "ns",]
sigCpG_sexDmr <- sex_dmr[sex_dmr$SigCpG != "ns",]

# Get m-values
sigCpG_stage <- mvals[sigCpG_stage$Name,]
sigCpG_sex <- mvals[sigCpG_sex$Name,]

## Read in the metaedata
preProccessedData <- readRDS(paste0(data_dir,"preproccessed_methyl.RDS"))
metadata <- preProccessedData$metadata
metadata$Sample_Name <- paste0("X",metadata$Sample_Name)
metadata <- metadata[match(colnames(stage_dmr[,grepl("^X",colnames(stage_dmr))]),metadata$Sample_Name),]

# Load CpG annotations
data(IlluminaHumanMethylationEPICanno.ilm10b4.hg19)
annotdata <- getAnnotation(IlluminaHumanMethylationEPICanno.ilm10b4.hg19)

```

## Prepare the data for clustering

```{r}

# Generate expdata
sigCpGStage_expData <- GetExpData(sigCpG_stage, metadata)
sigCpGSex_expData <- GetExpData(sigCpG_sex, metadata)
stage_expData <- GetExpData(stage_dmr[,grepl("^X",colnames(stage_dmr))], metadata)
sex_expData <- GetExpData(sex_dmr[,grepl("^X",colnames(sex_dmr))], metadata)

```

## Define the cluster fuzziness

Based off of generating randomized data and testing the minimum m values at which no clusters are detected.

```{r}

set.seed(seed_val)
stage_m1 <- mestimate(stage_expData)
sex_m1 <- mestimate(sex_expData)
sigCpGStage_m1 <- mestimate(sigCpGStage_expData)
sigCpGSex_m1 <- mestimate(sigCpGSex_expData)

# Get list of expdatas
expDatas <- list(Stage=stage_expData, SigStage=sigCpGStage_expData, Sex=sex_expData, SigSex=sigCpGSex_expData)
m1s <- list(Stage=stage_m1, SigStage=sigCpGStage_m1, Sex=sex_m1, SigSex=sigCpGSex_m1)

```

## Defining N clusters

C selection a bit naff, but Dmin suggests 3.

```{r}

# Toggle if we want to run D min calculation for optimal number of clusters
if (dmin) {
  
  # Loop over expdatas
  for (name in names(expDatas)) {
    
    print(name)
  
    # D min
    clustCheck_dDim <- Mfuzz::Dmin(eset=expDatas[[name]], m=m1s[[name]], crange=seq(min_clust,15,1), repeats=3, visu=TRUE)
    
    # Save
    dir.create(paste0(mfz_dir, name, "/"), showWarnings=FALSE)
    plt <- ggplot(data.frame(Dmin=clustCheck_dDim, nClusts=min_clust:15), aes(x=nClusts, y=Dmin)) + 
      geom_point() + 
      geom_line() + 
      xlab("Number of Clusters") + ylab("D Min (Dissimilarity)") +
      scale_x_continuous(n.breaks=15) +
      theme_bw(base_size=16) + theme(axis.title=element_text(size=24)) 
    ggsave(paste0(mfz_dir, name, "/", name,"_dminElbowPlot_minClust",min_clust,".png"), plt, units="px", width=4000, height=2500)
    
  }
    
}

```

## Check optimal clusters visually

```{r}
#set seed results not dependent of running Dmin workflow or not

set.seed(seed_val)
# Loop over expdatas
for (name in names(expDatas)) {
  
  print(name)
  
  # Loop over clusters
  for (nclusts in min_clust:max_clust) {
    
    # Carry out clustering using the calculated cluster number and fuzziness parameters
    set.seed(seed_val)
    mfuzz <- mfuzz(expDatas[[name]], c=nclusts, m=m1s[[name]])
    
    # Post-process
    centroids_long <- PostProcessMFuzz(mfuzz)
    
    # Plot
    mfuzz_centroids_plot <- ggplot(centroids_long, aes(x=as.numeric(sample),y=value, group=cluster, colour=as.factor(cluster))) + 
      geom_line(lwd=1) +
      theme_cowplot() +
      scale_colour_brewer(palette = "Set1") +
      xlab("Weeks") +
      ylab("Expression") +
      labs(title=paste0("fuzzy c-means expression: ",nclusts," clusters"), color = "Cluster") 
    
    # Save
    save_plot(paste0(mfz_dir,name,"/",name,"_nClust",nclusts,"_mfuzzCentroidsPlot.png"),mfuzz_centroids_plot,base_height = 6,base_width = 6,  bg = "white")
    
  }

}

```

## Perform final clustering

```{r}

# Define number of clusters
chosenClusts <- list(Stage=3, SigStage=3, SigSex=5, Sex=5)
mfuzzClusts <- list(Stage=NULL, SigStage=NULL, SigSex=NULL, Sex=NULL)
all_centroids_long <- data.frame()

# Loop over expdatas
for (name in names(chosenClusts)) {

  # Carry out clustering using the calculated cluster number and fuzziness parameters
  set.seed(seed_val)
  mfuzz_cl <- mfuzz(expDatas[[name]],
                    c=chosenClusts[[name]],
                    m=m1s[[name]])

  # Convert centroid data into long dataframe
  centroids_long <- PostProcessMFuzz(mfuzz_cl)
  
  # Plot
  mfuzz_centroids_plot <- ggplot(centroids_long, aes(x=as.numeric(sample),y=value, group=cluster, colour=as.factor(cluster))) + 
      geom_line(lwd=1) +
      theme_cowplot() +
      scale_colour_brewer(palette = "Set1") +
      xlab("Weeks") +
      ylab("Expression") +
      labs(title=paste0("fuzzy c-means expression: ",chosenClusts[[name]]," clusters"), color = "Cluster") 
  
  # Save
  save_plot(paste0(mfz_dir,name,"/",name,"_mFuzzCentroidsPlot_finalVersion.png"),mfuzz_centroids_plot,base_height = 6,base_width = 6,  bg = "white")
  mfuzz_centroids_plot
  
  # Append to list
  mfuzzClusts[[name]] <- mfuzz_cl
  
  # Append to output dataframe
  centroids_long["ClustSet"] <- name
  all_centroids_long <- rbind(all_centroids_long, centroids_long)
  
}

```

## Plot the cluster membership
```{r}

all_cluster_plot_df <- data.frame()

# Loop over expadatas
for (name in names(chosenClusts)) {

  # Create copy of dataframe containing average z-score per time
  membership_df <- data.frame(exprs(expDatas[[name]]))
  
  # Add row with feature names
  membership_df$Feature <- row.names(membership_df)
  
  # Bind the cluster assignment 
  membership_df$cluster <- mfuzzClusts[[name]]$cluster
  
  # Fetch the membership for each feature in its top scoring cluster
  membership_df$membership <- sapply(1:length(membership_df$cluster),function(row){
    clust <- membership_df$cluster[row]
    mfuzzClusts[[name]]$membership[row,clust]
  })
  
  
  # Convert membership dataframe into long format
  # cluster_plot_df <- membership_df %>%
  #   dplyr::select(.,1:length(colnames(mfuzz_cl$cluster)), membership, Feature, cluster) %>%
  #   tidyr::gather(.,"sample",'value',1:length(colnames(mfuzz_cl$cluster)))
  cluster_plot_df <- reshape2::melt(membership_df, id.vars=c("Feature","cluster","membership"))
  cluster_plot_df$sample <- as.character(gsub("X", "", cluster_plot_df$variable))
  
  # Order the dataframe by membership score\
  cluster_plot_df <- cluster_plot_df[order(cluster_plot_df$membership),]
  
  # Convert the features into an ordered factor using forcats
  cluster_plot_df$Feature = forcats::fct_inorder(cluster_plot_df$Feature)
  
  # Filter only top features by membership
  cluster_plot_df <- cluster_plot_df[ cluster_plot_df$membership>=min_clust_membership,]
  cluster_plot_df["ClustSet"] <- name
  
  # Combine
  all_cluster_plot_df <- rbind(cluster_plot_df, all_cluster_plot_df)
  
}

```

```{r}

# Loop over expdatas
for (name in names(expDatas)[c(2,3)]) {
  
  mfuzz_plot_membership <- ggplot(all_cluster_plot_df[all_cluster_plot_df$ClustSet == name,], 
                                  aes(x=as.numeric(sample),y=value)) + 
          # Facet-wrap by cluster
          facet_wrap(~cluster, scales="free_y") +
          # Colour each panel by cluster
          # The oppacity of each feature is determined by membership
          geom_line(aes(colour=as.factor(cluster), alpha=membership, group =Feature), lwd=1) +
          scale_colour_brewer(palette = "Set1") +
          # Add the centroids as a black line
          geom_line(data=all_centroids_long[all_centroids_long$ClustSet == name,], 
                    aes(x=as.numeric(sample), y=value, group=cluster), lwd=1.2, color="black",inherit.aes=FALSE) +
          xlab(name) +
          ylab("Methylation") +
          labs(title="Per-cluster membership") + 
          theme_cowplot() + 
          theme(legend.position = "none")
  
  save_plot(paste0(mfz_dir,name,"/",name,"_mFuzzPlotMembership.png"),mfuzz_plot_membership,base_height = 6,base_width = 10,  bg = "white")
  
}

```
```{r}
mfuzz_plot_membership
```

### Assemble DMR outputs

```{r}

# Read in
StageDmr_df <- read.table(paste0(dmr_dir, "dmrff_results.stage_as_continuous.treat_log2FC",trt_log2fc_thr,".preClust_TallFormat.tsv"), sep="\t", header=TRUE)
SexDmr_df <- read.table(paste0(dmr_dir, "dmrff_results.sex.preClust_TallFormat.tsv"), sep="\t", header=TRUE)

# Combine with clustering
StageDmr_df["CLUSTER"] <- 0
SexDmr_df["CLUSTER"] <- 0
rownames(StageDmr_df) <- StageDmr_df$CpG
rownames(SexDmr_df) <- SexDmr_df$CpG
StageDmr_df[intersect(names(mfuzzClusts$SigStage$cluster), rownames(StageDmr_df)),]["CLUSTER"] <- unname(mfuzzClusts$SigStage$cluster[intersect(names(mfuzzClusts$SigStage$cluster), rownames(StageDmr_df))])
SexDmr_df[intersect(names(mfuzzClusts$SigSex$cluster), rownames(SexDmr_df)),]["CLUSTER"] <- unname(mfuzzClusts$SigSex$cluster[intersect(names(mfuzzClusts$SigSex$cluster), rownames(SexDmr_df))])
table(StageDmr_df$CLUSTER); intersect(StageDmr_df[StageDmr_df$CLUSTER == 0,]$CpG, rownames(sigCpG_stage))
table(SexDmr_df$CLUSTER); intersect(SexDmr_df[SexDmr_df$CLUSTER == 0,]$CpG, rownames(sigCpG_sex))

# Assemble to output
StageDmr_df <- StageDmr_df[,
  c("DMR_ID","SITE","chr","start","end","SYMBOL","DISTANCE","B","S","estimate","se","z","p.value","p.adjust","CONTRAST","MAXGAP","CpG","pos","CLUSTER","logFC","P.Value","adj.P.Val","probe_bonferroni","n","nSigCpg","nNonSigCpg","pSigCpg")]
colnames(StageDmr_df)[c(13,14,20,21,22,23)] <- c("DMR_PVAL","DMR_PADJ","CPG_ESTIMATE","CPG_PVAL","CPG_PADJ","CPG_BONFERRONI")
SexDmr_df <- SexDmr_df[,
  c("DMR_ID","SITE","chr","start","end","SYMBOL","DISTANCE","B","S","estimate","se","z","p.value","p.adjust","CONTRAST","MAXGAP","CpG","pos","CLUSTER","logFC","P.Value","adj.P.Val","probe_bonferroni","n","nSigCpg","nNonSigCpg","pSigCpg")]
colnames(SexDmr_df)[c(13,14,20,21,22,23)] <- c("DMR_PVAL","DMR_PADJ","CPG_ESTIMATE","CPG_PVAL","CPG_PADJ","CPG_BONFERRONI")

# Save
write.table(StageDmr_df[order(as.numeric(gsub("dDMR_","",StageDmr_df$DMR_ID)), decreasing=FALSE),], 
            file=paste0(dmr_dir, "dmrff_results.stage_as_continuous.treat_log2FC",trt_log2fc_thr,".withCpG_mValClusters_TallFormat.tsv"),
            col.names = TRUE,row.names = FALSE, sep = "\t", quote = FALSE)
write.table(SexDmr_df[order(as.numeric(gsub("sDMR_","",SexDmr_df$DMR_ID)), decreasing=FALSE),],
            file=paste0(dmr_dir, "dmrff_results.sex.withCpG_mValClusters_TallFormat.tsv"),
            col.names = TRUE,row.names = FALSE, sep = "\t", quote = FALSE)

```

### Assemble CpG outputs

```{r}

# Read in
CpGTrt_df <- read.table(paste0(de_dir, "differential_methylation.limma_results.continuous_covariates.treat_log2FC",trt_log2fc_thr,"_thr.tsv"), sep="\t", header=TRUE, quote = "")

CpGTrt_df <- CpGTrt_df[ CpGTrt_df$bonferroni <= 0.05,]

CpGEby_df <- read.table(paste0(de_dir, "differential_methylation.limma_results.continuous_covariates.ebayes_adjusted.padj_thr",0.05,".tsv"), sep="\t", header=TRUE, quote = "")
StageCpG_df <- CpGTrt_df[CpGTrt_df$CONTRAST == "Z_STAGE_WEEKS",]
SexCpG_df <- CpGEby_df[CpGEby_df$CONTRAST == "sexmale - sexfemale",]

# Drop non-ID CpGs
# StageCpG_df <- StageCpG_df[grepl("^cg",StageCpG_df$Name),]
# SexCpG_df <- SexCpG_df[grepl("^cg",SexCpG_df$Name),]

# Get sig
rownames(StageCpG_df) <- StageCpG_df$Name
rownames(SexCpG_df) <- SexCpG_df$Name

# Combine with clustering
StageCpG_df["CLUSTER"] <- 0
SexCpG_df["CLUSTER"] <- 0
StageCpG_df[intersect(names(mfuzzClusts$SigStage$cluster), rownames(StageCpG_df)),]["CLUSTER"] <- unname(mfuzzClusts$SigStage$cluster[intersect(names(mfuzzClusts$SigStage$cluster), rownames(StageCpG_df))])
SexCpG_df[intersect(names(mfuzzClusts$SigSex$cluster), rownames(SexCpG_df)),]["CLUSTER"] <- unname(mfuzzClusts$SigSex$cluster[intersect(names(mfuzzClusts$SigSex$cluster), rownames(SexCpG_df))])
table(StageCpG_df$CLUSTER); intersect(StageCpG_df[StageCpG_df$CLUSTER == 0,]$Name, rownames(StageCpG_df))
table(SexCpG_df$CLUSTER); intersect(SexCpG_df[SexCpG_df$CLUSTER == 0,]$Name, rownames(SexCpG_df))

# Make cluster non-numeric
StageCpG_df["CLUSTER"] <- paste0("Cluster_",StageCpG_df$CLUSTER)
SexCpG_df["CLUSTER"] <- paste0("Cluster_",SexCpG_df$CLUSTER)

# Save
write.table(StageCpG_df[order(as.numeric(gsub("cg|ch[.][0-9XYMT]*[.]|F|M|R","",StageCpG_df$Name)), decreasing=FALSE),], 
            file=paste0(dmr_dir, "deCpG.stage_as_continuous.treat_log2FC",trt_log2fc_thr,".withCpG_mValClusters_TallFormat.tsv"),
            col.names = TRUE,row.names = FALSE, sep = "\t", quote = FALSE)
write.table(SexCpG_df[order(as.numeric(gsub("cg|ch","",SexCpG_df$Name)), decreasing=FALSE),],
            file=paste0(dmr_dir, "deCpG.sex.withCpG_mValClusters_TallFormat.tsv"),
            col.names = TRUE,row.names = FALSE, sep = "\t", quote = FALSE)

```

### Save beds

```{r}

# Save bed
bedStageDmr_df <- StageDmr_df[order(as.numeric(gsub("dDMR_","",StageDmr_df$DMR_ID)), decreasing=FALSE),c("chr","start","end","DMR_ID")]
bedSexDmr_df <- SexDmr_df[order(as.numeric(gsub("sDMR_","",SexDmr_df$DMR_ID)), decreasing=FALSE),c("chr","start","end","DMR_ID")]
write.table(bedStageDmr_df[!duplicated(bedStageDmr_df$DMR_ID),], 
            file=paste0(dmr_dir, "dmrff_results.stage_as_continuous.treat_log2FC",trt_log2fc_thr,".bed"),
            col.names = FALSE,row.names = FALSE, sep = "\t", quote = FALSE)
write.table(bedSexDmr_df[!duplicated(bedSexDmr_df$DMR_ID),], 
            file=paste0(dmr_dir, "dmrff_results.sex.bed"),
            col.names = FALSE,row.names = FALSE, sep = "\t", quote = FALSE)


# Save bed - with gene IDs
bedStageDmr_df <- StageDmr_df[order(as.numeric(gsub("dDMR_","",StageDmr_df$DMR_ID)), decreasing=FALSE),c("chr","start","end","SYMBOL","DMR_ID")]
bedSexDmr_df <- SexDmr_df[order(as.numeric(gsub("sDMR_","",SexDmr_df$DMR_ID)), decreasing=FALSE),c("chr","start","end","SYMBOL","DMR_ID")]
write.table(bedStageDmr_df[!duplicated(bedStageDmr_df$DMR_ID),!grepl("DMR_ID",colnames(bedStageDmr_df))], 
            file=paste0(dmr_dir, "dmrff_results.stage_as_continuous.treat_log2FC",trt_log2fc_thr,".withGeneIDs.bed"),
            col.names = FALSE,row.names = FALSE, sep = "\t", quote = FALSE)
write.table(bedSexDmr_df[!duplicated(bedSexDmr_df$DMR_ID),!grepl("DMR_ID",colnames(bedSexDmr_df))], 
            file=paste0(dmr_dir, "dmrff_results.sex.withGeneIDs.bed"),
            col.names = FALSE,row.names = FALSE, sep = "\t", quote = FALSE)

```

### Enrichment

#### Read in background

```{r}

# Read in unfiltered background
backSet <- read.table(paste0(de_dir, "differential_methylation.limma_results.continuous_covariates.treat_log2FC",trt_log2fc_thr,"_thr.tsv"), sep="\t", header=TRUE, quote = "")

# Filter to keep only stage
backSet <- backSet[backSet$CONTRAST == "Z_STAGE_WEEKS",]

```

#### CpG

```{r}

# Get vector on which to filter significant results
sig_cpgs <- rep(1, length(unique(StageCpG_df$CLUSTER))); names(sig_cpgs) <- unique(StageCpG_df$CLUSTER)

# Just use all CpGs in the cluster
for (direction in c("all")) {
  
  # Get up/down-reg
  if (direction == "up") { 
    
    stageCpG <- StageCpG_df[StageCpG_df$logFC > 0,] 
    stageCpG$CLUSTER <- as.character(stageCpG$CLUSTER)
    
    # CpG-level
    print(paste0("RUNNING FOR :- ",direction))
    up_goCpGStage <- RunMethylEnrichment(stageCpG, groupCol="CLUSTER", cpg_df=backSet[,c("Name","CONTRAST")], groupIsBackGroundSubset=FALSE,
                                               sig_cpgs=sig_cpgs[unique(stageCpG$CLUSTER)], cpg_level=TRUE, array.type="EPIC",
                                               pval_col="bonferroni", dir=goclust_dir, CpGCol="Name", extra=paste0(".",direction,"CpGStage"))
     
  # ELse if down
  } else if (direction == "down") { 
    
    stageCpG <- StageCpG_df[StageCpG_df$logFC < 0,] 
    
    # CpG-level
    print(paste0("RUNNING FOR :- ",direction))
    down_goCpGStage <- RunMethylEnrichment(stageCpG, groupCol="CLUSTER", cpg_df=backSet[,c("Name","CONTRAST")], groupIsBackGroundSubset=FALSE,
                                                 sig_cpgs=sig_cpgs[unique(stageCpG$CLUSTER)], cpg_level=TRUE, array.type="EPIC",
                                                 pval_col="bonferroni", dir=goclust_dir, CpGCol="Name", extra=paste0(".",direction,"CpGStage"))
    
  # Else if all comparisons
  } else if (direction == "all") {
    
    stageCpG <- StageCpG_df
    
    # CpG-level
    print(paste0("RUNNING FOR :- ",direction))
    all_goCpGStage <- RunMethylEnrichment(stageCpG, groupCol="CLUSTER", cpg_df=backSet[,c("Name","CONTRAST")], groupIsBackGroundSubset=FALSE,
                                                sig_cpgs=sig_cpgs[unique(stageCpG$CLUSTER)], cpg_level=TRUE, array.type="EPIC",
                                                pval_col="bonferroni", dir=goclust_dir, CpGCol="Name", extra=paste0(".",direction,"CpGStage"))
    
  }
  
}

```

### Plot

#### CpGs

```{r }

# Plot
## Loop over clusters
for (clust in unique(all_goCpGStage$GROUP)) {
  
  ### All
  PlotTopEnrichTerms(df=all_goCpGStage[all_goCpGStage$ONTOLOGY == "BP" & all_goCpGStage$GROUP == clust,], 
                     sortCol="FDR", topN=10, outDir=paste0(goclust_dir, "/", clust, "/"), extraDetails=paste0("_",clust,"_allStageCpGs_mFuzzClusts_goBP"), width=2500, height=2250)
  PlotTopEnrichTerms(df=all_goCpGStage[all_goCpGStage$ONTOLOGY == "MF" & all_goCpGStage$GROUP == clust,], 
                     sortCol="FDR", topN=10, outDir=paste0(goclust_dir, "/", clust, "/"), extraDetails=paste0("_",clust,"_allStageCpGs_mFuzzClusts_goMF"), width=2500, height=2250)
  PlotTopEnrichTerms(df=all_goCpGStage[all_goCpGStage$ONTOLOGY == "CC" & all_goCpGStage$GROUP == clust,], 
                     sortCol="FDR", topN=10, outDir=paste0(goclust_dir, "/", clust, "/"), extraDetails=paste0("_",clust,"_allStageCpGs_mFuzzClusts_goCC"), width=2500, height=2250)

}

```