knitr::opts_chunk$set(echo = TRUE, collapse = TRUE, warning = FALSE)

Condition_1 <- params$Condition_1
Condition_2 <- params$Condition_2

Top GO Analysis for hEGCLCs vs hiPSCs

For this functional we will focus only on ALL differentially methylated genomic regions of protein coding genes!

1. Environment Set Up

library(DT)
library(ggplot2)
library(AnnotationDbi)
library(org.Hs.eg.db)
library(gridExtra)
library(RColorBrewer)
library(viridis)
library(topGO)
library(data.table)
library(tidyr)
library(dplyr)

source("TopGO_helper_meth.R")
OutputFolder <- ifelse(is.null(params$OutputFolder), getwd(), params$OutputFolder) 


if (dir.exists(OutputFolder) == FALSE) {
  dir.create(OutputFolder, recursive=TRUE)
}

2. Data Load

DMRAnnotated <- readRDS(params$DMRAnnotated)
GeneUniverse_df <- readRDS(params$GeneUniverse)
GeneUniverse_df_protein_coding <- GeneUniverse_df[GeneUniverse_df$gene_biotype %in% "protein_coding",]
GeneUniverse <- unique(GeneUniverse_df_protein_coding$hgnc_symbol)

We have 18940 protein-coding genes in the universe.

Saving of selected DMRs

DMReg <- DMRAnnotated[[paste0(Condition_1,"vs",Condition_2)]]
DMReg <- DMReg %>% filter(gene_biotype == "protein_coding")

3. TOPGO for Gene Ontology Enrichment analysis

Gene ontology enrichment analysis is performed on the set of DMGs (Differentially Methylates Genes) using TopGO with Fisher statistics and weight01 algorithm. DMGs are now splitted in down- and up-methylated.

3.1 Selection of modulated genes and generation of gene vectors

I generate vectors for the gene universe, all modulated genes, up-meth genes and down-meth genes in the format required by TopGO.

GeneVectors <- topGOGeneVectors_meth(DMRAnnotated, Cond1 = Condition_1, Cond2 = Condition_2, genomic_type = "all", gene_type = "protein_coding", Universe = GeneUniverse)
## [1] "All the genes associated to DMRs are in the gene universe!"

Therefore:

  • universe genes: 18940 genes
  • modulated genes: 369 genes
  • down-methylated genes: 9 genes of interest
  • up-methylated genes: 360 genes of interest

N.B.: The sum of down- and up- methylated genes is > the number of modulated genes because some genes are shared between up and down lists. This is because for example in a promoter of a gene meth can be higher for a cell type and in the gene body of the same gene meth can be higher for the other cell type!

BpEval <- ifelse(length(grep('BP', params$TopGO))!=0, TRUE, FALSE)
MfEval <- ifelse(length(grep('MF', params$TopGO))!=0, TRUE, FALSE)
CcEval <- ifelse(length(grep('CC', params$TopGO))!=0, TRUE, FALSE)

3.2 TopGO analysis: Biological Process

On the basis of the analysis settings, the enrichment for Biological Process IS performed.

Biological Process Analysis for ALL modulated genes: 369 genes

# I generate a list that contains the association between each gene and the GO terms that are associated to it
BPann <- topGO::annFUN.org(whichOnto="BP", feasibleGenes=names(GeneVectors$DMGenes), 
                           mapping="org.Hs.eg.db", ID="symbol") %>% inverseList()

# Wrapper function for topGO analysis 
ResBPAll <- topGOResults_new(Genes=GeneVectors$DMGenes, gene2GO=BPann, ontology='BP', 
                         desc=NULL, nodeSize=15, algorithm='weight01', statistic='fisher', 
                         EnTh=params$GoEnTh, PvalTh=params$GoPvalTh, minTerms=12, geneTh=4,
                         saveRes=TRUE, outDir=paste0(OutputFolder), fileName='BPAll', maxAnn = 500)
## Gene vector contains levels: 0,1
## 
## Building most specific GOs .....
##  ( 12448 GO terms found. )
## 
## Build GO DAG topology ..........
##  ( 15834 GO terms and 35900 relations. )
## 
## Annotating nodes ...............
##  ( 16734 genes annotated to the GO terms. )
## 
##           -- Weight01 Algorithm -- 
## 
##       the algorithm is scoring 3458 nontrivial nodes
##       parameters: 
##           test statistic: fisher
## 
##   Level 16:  2 nodes to be scored    (0 eliminated genes)
## 
##   Level 15:  6 nodes to be scored    (0 eliminated genes)
## 
##   Level 14:  32 nodes to be scored   (49 eliminated genes)
## 
##   Level 13:  59 nodes to be scored   (185 eliminated genes)
## 
##   Level 12:  116 nodes to be scored  (1254 eliminated genes)
## 
##   Level 11:  208 nodes to be scored  (3519 eliminated genes)
## 
##   Level 10:  326 nodes to be scored  (5858 eliminated genes)
## 
##   Level 9:   459 nodes to be scored  (7499 eliminated genes)
## 
##   Level 8:   513 nodes to be scored  (9796 eliminated genes)
## 
##   Level 7:   588 nodes to be scored  (12148 eliminated genes)
## 
##   Level 6:   526 nodes to be scored  (14173 eliminated genes)
## 
##   Level 5:   334 nodes to be scored  (15345 eliminated genes)
## 
##   Level 4:   188 nodes to be scored  (16104 eliminated genes)
## 
##   Level 3:   83 nodes to be scored   (16375 eliminated genes)
## 
##   Level 2:   17 nodes to be scored   (16494 eliminated genes)
## 
##   Level 1:   1 nodes to be scored    (16556 eliminated genes)

Biological Process Analysis for DOWN-METHYLATED genes: 9 genes

# Wrapper function for topGO analysis 
ResBPDown <- topGOResults_new(Genes=GeneVectors$DMGenesDown, gene2GO=BPann, ontology='BP', 
                          desc=NULL, nodeSize=15, algorithm='weight01', statistic='fisher', 
                          EnTh=params$GoEnTh, PvalTh=params$GoPvalTh, minTerms=12, geneTh=4,
                          saveRes=TRUE, outDir=paste0(OutputFolder), fileName='BPDown', maxAnn = 500)
## Gene vector contains levels: 0,1
## 
## Building most specific GOs .....
##  ( 12448 GO terms found. )
## 
## Build GO DAG topology ..........
##  ( 15834 GO terms and 35900 relations. )
## 
## Annotating nodes ...............
##  ( 16734 genes annotated to the GO terms. )
## 
##           -- Weight01 Algorithm -- 
## 
##       the algorithm is scoring 681 nontrivial nodes
##       parameters: 
##           test statistic: fisher
## 
##   Level 14:  1 nodes to be scored    (0 eliminated genes)
## 
##   Level 13:  5 nodes to be scored    (0 eliminated genes)
## 
##   Level 12:  15 nodes to be scored   (20 eliminated genes)
## 
##   Level 11:  28 nodes to be scored   (1833 eliminated genes)
## 
##   Level 10:  39 nodes to be scored   (3209 eliminated genes)
## 
##   Level 9:   65 nodes to be scored   (3978 eliminated genes)
## 
##   Level 8:   72 nodes to be scored   (5322 eliminated genes)
## 
##   Level 7:   102 nodes to be scored  (7068 eliminated genes)
## 
##   Level 6:   118 nodes to be scored  (9032 eliminated genes)
## 
##   Level 5:   114 nodes to be scored  (11140 eliminated genes)
## 
##   Level 4:   70 nodes to be scored   (14221 eliminated genes)
## 
##   Level 3:   39 nodes to be scored   (15736 eliminated genes)
## 
##   Level 2:   12 nodes to be scored   (16306 eliminated genes)
## 
##   Level 1:   1 nodes to be scored    (16519 eliminated genes)

# Selection on enrichment of at least 2 is implemented (also to avoid depleted categories). Then categories are ranked by PVal and all the ones with Pval < th are selected. If the number is < minTerms, othter terms are included to reach the minimum number.

Biological Process Analysis for UP-METHYLATED genes: 360 genes

ResBPUp <- topGOResults_new(Genes=GeneVectors$DMGenesUp, gene2GO=BPann, ontology='BP', 
                        desc=NULL, nodeSize=15, algorithm='weight01', statistic='fisher', 
                        EnTh=params$GoEnTh, PvalTh=params$GoPvalTh, minTerms=12, geneTh=4,
                        saveRes=TRUE, outDir=OutputFolder, fileName='BPUp', maxAnn = 500)
## Gene vector contains levels: 0,1
## 
## Building most specific GOs .....
##  ( 12448 GO terms found. )
## 
## Build GO DAG topology ..........
##  ( 15834 GO terms and 35900 relations. )
## 
## Annotating nodes ...............
##  ( 16734 genes annotated to the GO terms. )
## 
##           -- Weight01 Algorithm -- 
## 
##       the algorithm is scoring 3379 nontrivial nodes
##       parameters: 
##           test statistic: fisher
## 
##   Level 16:  2 nodes to be scored    (0 eliminated genes)
## 
##   Level 15:  6 nodes to be scored    (0 eliminated genes)
## 
##   Level 14:  31 nodes to be scored   (49 eliminated genes)
## 
##   Level 13:  58 nodes to be scored   (185 eliminated genes)
## 
##   Level 12:  110 nodes to be scored  (1240 eliminated genes)
## 
##   Level 11:  199 nodes to be scored  (3512 eliminated genes)
## 
##   Level 10:  315 nodes to be scored  (5781 eliminated genes)
## 
##   Level 9:   446 nodes to be scored  (7431 eliminated genes)
## 
##   Level 8:   500 nodes to be scored  (9767 eliminated genes)
## 
##   Level 7:   574 nodes to be scored  (12136 eliminated genes)
## 
##   Level 6:   518 nodes to be scored  (14157 eliminated genes)
## 
##   Level 5:   333 nodes to be scored  (15341 eliminated genes)
## 
##   Level 4:   187 nodes to be scored  (16085 eliminated genes)
## 
##   Level 3:   82 nodes to be scored   (16375 eliminated genes)
## 
##   Level 2:   17 nodes to be scored   (16494 eliminated genes)
## 
##   Level 1:   1 nodes to be scored    (16556 eliminated genes)

Result visualization: Barplot

nterms=15
cols <- c(All='forestgreen', Down='blue', Up='red')
pals <- lapply(cols, {function(el) grDevices::colorRampPalette(c('gold', el))(nterms)})
if(nrow(ResBPAll$ResSel)>0){
  topGOBarplot_new(TopGORes=ResBPAll$ResSel, terms=nterms, pvalTh=0.01, plotTitle=paste0("All DMGenes - ", Condition_2, " vs ", Condition_1),
                             palette=pals$All, flip_x=FALSE)
}  

if(nrow(ResBPDown$ResSel)>0){
  topGOBarplot_new(TopGORes=ResBPDown$ResSel, terms=nterms, pvalTh=0.01, plotTitle=paste0("Down-methylated genes (associated to DMRs where methylation is lower in ", Condition_2, ")"),
                             palette=pals$Down, flip_x=FALSE)
}
if(nrow(ResBPUp$ResSel)>0){
  topGOBarplot_new(TopGORes=ResBPUp$ResSel, terms=nterms, pvalTh=0.01, plotTitle=paste0("Up-methylated genes (associated to DMRs where methylation is higher in ", Condition_2, ")"),
                             palette=pals$Up, flip_x=TRUE)
}

if(nrow(ResBPAll$ResSel)>0 && nrow(ResBPDown$ResSel)>0 && nrow(ResBPUp$ResSel)>0){
  topGOBarplotAll_new(TopGOResAll=ResBPAll$ResSel, TopGOResDown=ResBPDown$ResSel, TopGOResUp=ResBPUp$ResSel, 
                terms=nterms, pvalTh=0.01, plotTitle="TopGO analysis of DMGs - ")
}

Result visualization: Bubbleplot

if(nrow(ResBPAll$ResSel)>0){
  bubbleplot(ResBPAll$ResSel, terms=nterms, Ont = "BP", pvalTh=0.01, plotTitle="TopGO analysis of DMGs - All")
  }

if(nrow(ResBPDown$ResSel)>0){
  bubbleplot(ResBPDown$ResSel, terms=nterms, Ont = "BP", pvalTh=0.01, plotTitle="TopGO analysis of DMGs - Down", fillCol = "blue")
  }
if(nrow(ResBPUp$ResSel)>0){
  bubbleplot(ResBPUp$ResSel, terms=nterms, Ont = "BP", pvalTh=0.01, plotTitle="TopGO analysis of DMGs - Up", fillCol = "red")
  }

Result visualization: Top genes in top GO terms

if(nrow(ResBPAll$ResSel)>0){
  plotGenesInTerm_v2(ResBPAll$ResSel, ResBPAll$GOdata, DMRs = DMReg, nterms=nterms, ngenes=12, plotTitle="Top genes in top GO terms - All", Interactive=FALSE)
}

if(nrow(ResBPDown$ResSel)>0){
  plotGenesInTerm_v2(ResBPDown$ResSel, ResBPAll$GOdata, DMRs = DMReg, nterms=nterms, ngenes=12, plotTitle="Top genes in top GO terms - Down", Interactive=FALSE, fillCol='blue')
}  
if(nrow(ResBPUp$ResSel)>0){
  plotGenesInTerm_v2(ResBPUp$ResSel, ResBPAll$GOdata, DMRs = DMReg, nterms=nterms, ngenes=12, plotTitle="Top genes in top GO terms - Up", Interactive=FALSE, fillCol='red')
}

3.3 TopGO analysis: Molecular Function

On the basis of the analysis settings, the enrichment for Molecular Function IS performed.

Molecular Function Enrichment for ALL modulated genes: 369 genes

# I generate a list that contains the association between each gene and the GO terms that are associated to it
MFann <- topGO::annFUN.org(whichOnto='MF', feasibleGenes=names(GeneVectors$DMGenes), 
                           mapping='org.Hs.eg.db', ID='symbol') %>% inverseList()

# Wrapper function for topGO analysis 
ResMFAll <- topGOResults_new(Genes=GeneVectors$DMGenes, gene2GO=MFann, ontology='MF', 
                         desc=NULL, nodeSize=15, algorithm='weight01', statistic='fisher', 
                         EnTh=params$GoEnTh, PvalTh=params$GoPvalTh, minTerms=12, geneTh=4,
                         saveRes=TRUE, outDir=OutputFolder, fileName='MFAll', maxAnn = 500)
## Gene vector contains levels: 0,1
## 
## Building most specific GOs .....
##  ( 4495 GO terms found. )
## 
## Build GO DAG topology ..........
##  ( 4939 GO terms and 6405 relations. )
## 
## Annotating nodes ...............
##  ( 17198 genes annotated to the GO terms. )
## 
##           -- Weight01 Algorithm -- 
## 
##       the algorithm is scoring 572 nontrivial nodes
##       parameters: 
##           test statistic: fisher
## 
##   Level 11:  8 nodes to be scored    (0 eliminated genes)
## 
##   Level 10:  15 nodes to be scored   (0 eliminated genes)
## 
##   Level 9:   26 nodes to be scored   (200 eliminated genes)
## 
##   Level 8:   37 nodes to be scored   (1526 eliminated genes)
## 
##   Level 7:   77 nodes to be scored   (3773 eliminated genes)
## 
##   Level 6:   98 nodes to be scored   (4358 eliminated genes)
## 
##   Level 5:   124 nodes to be scored  (6255 eliminated genes)
## 
##   Level 4:   125 nodes to be scored  (9827 eliminated genes)
## 
##   Level 3:   48 nodes to be scored   (13123 eliminated genes)
## 
##   Level 2:   13 nodes to be scored   (14529 eliminated genes)
## 
##   Level 1:   1 nodes to be scored    (17012 eliminated genes)

Molecular Function Enrichment for DOWN-METHYLATED genes: 9 genes

ResMFDown <- topGOResults_new(Genes=GeneVectors$DMGenesDown, gene2GO=MFann, ontology='MF', 
                          desc=NULL, nodeSize=15, algorithm='weight01', statistic='fisher', 
                          EnTh=params$GoEnTh, PvalTh=params$GoPvalTh, minTerms=12, geneTh=4,
                          saveRes=TRUE, outDir=OutputFolder, fileName='MFDown', maxAnn = 500)
## Gene vector contains levels: 0,1
## 
## Building most specific GOs .....
##  ( 4495 GO terms found. )
## 
## Build GO DAG topology ..........
##  ( 4939 GO terms and 6405 relations. )
## 
## Annotating nodes ...............
##  ( 17198 genes annotated to the GO terms. )
## 
##           -- Weight01 Algorithm -- 
## 
##       the algorithm is scoring 107 nontrivial nodes
##       parameters: 
##           test statistic: fisher
## 
##   Level 11:  1 nodes to be scored    (0 eliminated genes)
## 
##   Level 10:  2 nodes to be scored    (0 eliminated genes)
## 
##   Level 9:   5 nodes to be scored    (15 eliminated genes)
## 
##   Level 8:   9 nodes to be scored    (1197 eliminated genes)
## 
##   Level 7:   12 nodes to be scored   (1477 eliminated genes)
## 
##   Level 6:   15 nodes to be scored   (1914 eliminated genes)
## 
##   Level 5:   20 nodes to be scored   (2981 eliminated genes)
## 
##   Level 4:   22 nodes to be scored   (4469 eliminated genes)
## 
##   Level 3:   13 nodes to be scored   (8867 eliminated genes)
## 
##   Level 2:   7 nodes to be scored    (11594 eliminated genes)
## 
##   Level 1:   1 nodes to be scored    (16546 eliminated genes)

Molecular Function Analysis for UP-METHYLATED genes: 360 genes

ResMFUp <- topGOResults_new(Genes=GeneVectors$DMGenesUp, gene2GO=MFann, ontology='MF', 
                        desc=NULL, nodeSize=15, algorithm='weight01', statistic='fisher', 
                        EnTh=params$GoEnTh, PvalTh=params$GoPvalTh, minTerms=12, geneTh=4,
                        saveRes=TRUE, outDir=OutputFolder, fileName='MFUp', maxAnn = 500)
## Gene vector contains levels: 0,1
## 
## Building most specific GOs .....
##  ( 4495 GO terms found. )
## 
## Build GO DAG topology ..........
##  ( 4939 GO terms and 6405 relations. )
## 
## Annotating nodes ...............
##  ( 17198 genes annotated to the GO terms. )
## 
##           -- Weight01 Algorithm -- 
## 
##       the algorithm is scoring 568 nontrivial nodes
##       parameters: 
##           test statistic: fisher
## 
##   Level 11:  7 nodes to be scored    (0 eliminated genes)
## 
##   Level 10:  15 nodes to be scored   (0 eliminated genes)
## 
##   Level 9:   26 nodes to be scored   (185 eliminated genes)
## 
##   Level 8:   37 nodes to be scored   (1526 eliminated genes)
## 
##   Level 7:   77 nodes to be scored   (3773 eliminated genes)
## 
##   Level 6:   97 nodes to be scored   (4358 eliminated genes)
## 
##   Level 5:   123 nodes to be scored  (6255 eliminated genes)
## 
##   Level 4:   125 nodes to be scored  (9825 eliminated genes)
## 
##   Level 3:   48 nodes to be scored   (13112 eliminated genes)
## 
##   Level 2:   12 nodes to be scored   (14529 eliminated genes)
## 
##   Level 1:   1 nodes to be scored    (17012 eliminated genes)

Result visualization: Barplot

if(nrow(ResMFAll$ResSel)>0){
  topGOBarplot_new(TopGORes=ResMFAll$ResSel, terms=nterms, pvalTh=0.01, plotTitle=paste0("All DMGenes - ", Condition_2, " vs ", Condition_1),
                             palette=pals$All, flip_x=FALSE)
} 

if(nrow(ResMFDown$ResSel)>0){
  topGOBarplot_new(TopGORes=ResMFDown$ResSel, terms=nterms, pvalTh=0.01, plotTitle=paste0("Down-methylated genes (associated to DMRs where methylation is lower in ", Condition_2, ")"),
                             palette=pals$Down, flip_x=FALSE)
}  
if(nrow(ResMFUp$ResSel)>0){
  topGOBarplot_new(TopGORes=ResMFUp$ResSel, terms=nterms, pvalTh=0.01, plotTitle=paste0("Up-methylated genes (associated to DMRs where methylation is higher in ", Condition_2, ")"),
                             palette=pals$Up, flip_x=TRUE)
}

if(nrow(ResMFAll$ResSel)>0 && nrow(ResMFDown$ResSel)>0 && nrow(ResMFUp$ResSel)>0){
  topGOBarplotAll_new(TopGOResAll=ResMFAll$ResSel, TopGOResDown=ResMFDown$ResSel, TopGOResUp=ResMFUp$ResSel, 
                terms=nterms, pvalTh=0.01, plotTitle= "TopGO analysis of DMGs - ")
}

Result visualization: Bubbleplot

if(nrow(ResMFAll$ResSel)>0){
  bubbleplot(ResMFAll$ResSel, terms=nterms, Ont = "MF", pvalTh=0.01, plotTitle="TopGO analysis of DMGs - All")
  }

if(nrow(ResMFDown$ResSel)>0){
  bubbleplot(ResMFDown$ResSel, terms=nterms, Ont = "MF", pvalTh=0.01, plotTitle="TopGO analysis of DMGs - Down", fillCol = "blue")
  }
if(nrow(ResMFUp$ResSel)>0){
  bubbleplot(ResMFUp$ResSel, terms=nterms, Ont = "MF", pvalTh=0.01, plotTitle="TopGO analysis of DMGs - Up", fillCol = "red")
  }

Result visualization: Top genes in top GO terms

if(nrow(ResMFAll$ResSel)>0){
  plotGenesInTerm_v2(ResMFAll$ResSel, ResMFAll$GOdata, DMRs = DMReg, nterms=nterms, ngenes=12, plotTitle="Top genes in top GO terms - All", Interactive=FALSE)
}

if(nrow(ResMFDown$ResSel)>0){
  plotGenesInTerm_v2(ResMFDown$ResSel, ResMFAll$GOdata, DMRs = DMReg, nterms=nterms, ngenes=12, plotTitle="Top genes in top GO terms - Down", Interactive=FALSE, fillCol='blue')
}  
if(nrow(ResMFUp$ResSel)>0){
  plotGenesInTerm_v2(ResMFUp$ResSel, ResMFAll$GOdata, DMRs = DMReg, nterms=nterms, ngenes=12, plotTitle="Top genes in top GO terms - Up", Interactive=FALSE, fillCol='red')
}

3.4 TopGO analysis: Cellular Component

On the basis of the analysis settings, the enrichment for Cellular Component IS performed.

Cellular Component Enrichment for ALL modulated genes: 369 genes

# I generate a list that contains the association between each gene and the GO terms that are associated to it
CCann <- topGO::annFUN.org(whichOnto='CC', feasibleGenes=names(GeneVectors$DMGenes), 
                           mapping='org.Hs.eg.db', ID='symbol') %>% inverseList()

# Wrapper function for topGO analysis 
ResCCAll <- topGOResults_new(Genes=GeneVectors$DMGenes, gene2GO=CCann, ontology='CC', 
                         desc=NULL, nodeSize=15, algorithm='weight01', statistic='fisher', 
                         EnTh=params$GoEnTh, PvalTh=params$GoPvalTh, minTerms=12, geneTh=4,
                         saveRes=TRUE, outDir=OutputFolder, fileName='CCAll', maxAnn = 500)
## Gene vector contains levels: 0,1
## 
## Building most specific GOs .....
##  ( 1818 GO terms found. )
## 
## Build GO DAG topology ..........
##  ( 2014 GO terms and 3390 relations. )
## 
## Annotating nodes ...............
##  ( 17642 genes annotated to the GO terms. )
## 
##           -- Weight01 Algorithm -- 
## 
##       the algorithm is scoring 469 nontrivial nodes
##       parameters: 
##           test statistic: fisher
## 
##   Level 13:  1 nodes to be scored    (0 eliminated genes)
## 
##   Level 12:  1 nodes to be scored    (0 eliminated genes)
## 
##   Level 11:  15 nodes to be scored   (30 eliminated genes)
## 
##   Level 10:  38 nodes to be scored   (42 eliminated genes)
## 
##   Level 9:   67 nodes to be scored   (764 eliminated genes)
## 
##   Level 8:   72 nodes to be scored   (2726 eliminated genes)
## 
##   Level 7:   70 nodes to be scored   (5607 eliminated genes)
## 
##   Level 6:   68 nodes to be scored   (9834 eliminated genes)
## 
##   Level 5:   60 nodes to be scored   (11986 eliminated genes)
## 
##   Level 4:   38 nodes to be scored   (14607 eliminated genes)
## 
##   Level 3:   36 nodes to be scored   (16577 eliminated genes)
## 
##   Level 2:   2 nodes to be scored    (17315 eliminated genes)
## 
##   Level 1:   1 nodes to be scored    (17557 eliminated genes)

Cellular Component Enrichment for DOWN-METHYLATED genes: 9 genes

# Wrapper function for topGO analysis 
ResCCDown <- topGOResults_new(Genes=GeneVectors$DMGenesDown, gene2GO=CCann, ontology='CC', 
                          desc=NULL, nodeSize=15, algorithm='weight01', statistic='fisher', 
                          EnTh=params$GoEnTh, PvalTh=params$GoPvalTh, minTerms=12, geneTh=4,
                          saveRes=TRUE, outDir=OutputFolder, fileName='CCDown', maxAnn = 500)
## Gene vector contains levels: 0,1
## 
## Building most specific GOs .....
##  ( 1818 GO terms found. )
## 
## Build GO DAG topology ..........
##  ( 2014 GO terms and 3390 relations. )
## 
## Annotating nodes ...............
##  ( 17642 genes annotated to the GO terms. )
## 
##           -- Weight01 Algorithm -- 
## 
##       the algorithm is scoring 102 nontrivial nodes
##       parameters: 
##           test statistic: fisher
## 
##   Level 11:  3 nodes to be scored    (0 eliminated genes)
## 
##   Level 10:  4 nodes to be scored    (0 eliminated genes)
## 
##   Level 9:   1 nodes to be scored    (145 eliminated genes)
## 
##   Level 8:   6 nodes to be scored    (488 eliminated genes)
## 
##   Level 7:   9 nodes to be scored    (857 eliminated genes)
## 
##   Level 6:   19 nodes to be scored   (5796 eliminated genes)
## 
##   Level 5:   19 nodes to be scored   (8091 eliminated genes)
## 
##   Level 4:   18 nodes to be scored   (12751 eliminated genes)
## 
##   Level 3:   20 nodes to be scored   (16364 eliminated genes)
## 
##   Level 2:   2 nodes to be scored    (17286 eliminated genes)
## 
##   Level 1:   1 nodes to be scored    (17527 eliminated genes)

Cellular Component Analysis for UP-METHYLATED genes: 360 genes

# Wrapper function for topGO analysis 
ResCCUp <- topGOResults_new(Genes=GeneVectors$DMGenesUp, gene2GO=CCann, ontology='CC', 
                        desc=NULL, nodeSize=15, algorithm='weight01', statistic='fisher', 
                        EnTh=params$GoEnTh, PvalTh=params$GoPvalTh, minTerms=12, geneTh=4,
                        saveRes=TRUE, outDir=OutputFolder, fileName='CCUp', maxAnn = 500)
## Gene vector contains levels: 0,1
## 
## Building most specific GOs .....
##  ( 1818 GO terms found. )
## 
## Build GO DAG topology ..........
##  ( 2014 GO terms and 3390 relations. )
## 
## Annotating nodes ...............
##  ( 17642 genes annotated to the GO terms. )
## 
##           -- Weight01 Algorithm -- 
## 
##       the algorithm is scoring 466 nontrivial nodes
##       parameters: 
##           test statistic: fisher
## 
##   Level 13:  1 nodes to be scored    (0 eliminated genes)
## 
##   Level 12:  1 nodes to be scored    (0 eliminated genes)
## 
##   Level 11:  15 nodes to be scored   (30 eliminated genes)
## 
##   Level 10:  38 nodes to be scored   (42 eliminated genes)
## 
##   Level 9:   67 nodes to be scored   (764 eliminated genes)
## 
##   Level 8:   72 nodes to be scored   (2726 eliminated genes)
## 
##   Level 7:   69 nodes to be scored   (5607 eliminated genes)
## 
##   Level 6:   68 nodes to be scored   (9834 eliminated genes)
## 
##   Level 5:   59 nodes to be scored   (11980 eliminated genes)
## 
##   Level 4:   37 nodes to be scored   (14607 eliminated genes)
## 
##   Level 3:   36 nodes to be scored   (16577 eliminated genes)
## 
##   Level 2:   2 nodes to be scored    (17315 eliminated genes)
## 
##   Level 1:   1 nodes to be scored    (17557 eliminated genes)

Result visualization: Barplot

if(nrow(ResCCAll$ResSel)>0){
  topGOBarplot_new(TopGORes=ResCCAll$ResSel, terms=nterms, pvalTh=0.01, plotTitle=paste0("All DMGenes - ", Condition_2, " vs ", Condition_1),
                             palette=pals$All, flip_x=FALSE)
}  

if(nrow(ResCCDown$ResSel)>0){
  topGOBarplot_new(TopGORes=ResCCDown$ResSel, terms=nterms, pvalTh=0.01, plotTitle=paste0("Down-methylated genes (associated to DMRs where methylation is lower in ", Condition_2, ")"),
                             palette=pals$Down, flip_x=FALSE)
}
if(nrow(ResCCUp$ResSel)>0){
  topGOBarplot_new(TopGORes=ResCCUp$ResSel, terms=nterms, pvalTh=0.01, plotTitle=paste0("Up-methylated genes (associated to DMRs where methylation is higher in ", Condition_2, ")"), palette=pals$Up, flip_x=TRUE)
}

if(nrow(ResCCAll$ResSel)>0 && nrow(ResCCDown$ResSel)>0 && nrow(ResCCUp$ResSel)>0){
  topGOBarplotAll_new(TopGOResAll=ResCCAll$ResSel, TopGOResDown=ResCCDown$ResSel, TopGOResUp=ResCCUp$ResSel, terms=nterms, pvalTh=0.01, plotTitle="TopGO analysis of DMGs - ")
}

Result visualization: Bubbleplot

if(nrow(ResCCAll$ResSel)>0){
  bubbleplot(ResCCAll$ResSel, terms=nterms, Ont = "CC", pvalTh=0.01, plotTitle="TopGO analysis of DMGs - All")
  }

if(nrow(ResCCDown$ResSel)>0){
  bubbleplot(ResCCDown$ResSel, terms=nterms, Ont = "CC", pvalTh=0.01, plotTitle="TopGO analysis of DMGs - Down", fillCol = "blue")
  }
if(nrow(ResCCUp$ResSel)>0){
  bubbleplot(ResCCUp$ResSel, terms=nterms, Ont = "CC", pvalTh=0.01, plotTitle="TopGO analysis of DMGs - Up", fillCol = "red")
  }

Result visualization: Top genes in top GO terms

if(nrow(ResCCAll$ResSel)>0){
  plotGenesInTerm_v2(ResCCAll$ResSel, ResCCAll$GOdata, DMRs = DMReg, nterms=nterms, ngenes=12, plotTitle="Top genes in top GO terms - All", Interactive=FALSE)
}

if(nrow(ResCCDown$ResSel)>0){
  plotGenesInTerm_v2(ResCCDown$ResSel, ResCCAll$GOdata, DMRs = DMReg, nterms=nterms, ngenes=12, plotTitle="Top genes in top GO terms - Down", Interactive=FALSE, fillCol='blue')
}  
if(nrow(ResCCUp$ResSel)>0){
  plotGenesInTerm_v2(ResCCUp$ResSel, ResCCAll$GOdata, DMRs = DMReg, nterms=nterms, ngenes=12, plotTitle="Top genes in top GO terms - Up", Interactive=FALSE, fillCol='red')
}

4. Savings

Most of the useful information has been saved during the analysis. Here I save the workspace and information about the session.

SessionInfo <- sessionInfo()
Date <- date()
Date
## [1] "Fri Jan 10 12:30:18 2025"
SessionInfo
## R version 4.2.1 (2022-06-23)
## Platform: x86_64-pc-linux-gnu (64-bit)
## Running under: Ubuntu 20.04.4 LTS
## 
## Matrix products: default
## BLAS:   /usr/lib/x86_64-linux-gnu/blas/libblas.so.3.9.0
## LAPACK: /usr/lib/x86_64-linux-gnu/lapack/liblapack.so.3.9.0
## 
## locale:
##  [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C               LC_TIME=en_US.UTF-8        LC_COLLATE=en_US.UTF-8     LC_MONETARY=en_US.UTF-8    LC_MESSAGES=en_US.UTF-8   
##  [7] LC_PAPER=en_US.UTF-8       LC_NAME=C                  LC_ADDRESS=C               LC_TELEPHONE=C             LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C       
## 
## attached base packages:
## [1] stats4    stats     graphics  grDevices datasets  utils     methods   base     
## 
## other attached packages:
##  [1] dplyr_1.1.2          tidyr_1.3.0          data.table_1.14.8    topGO_2.50.0         SparseM_1.81         GO.db_3.16.0         graph_1.76.0         viridis_0.6.2       
##  [9] viridisLite_0.4.2    RColorBrewer_1.1-3   gridExtra_2.3        org.Hs.eg.db_3.16.0  AnnotationDbi_1.60.2 IRanges_2.32.0       S4Vectors_0.36.2     Biobase_2.58.0      
## [17] BiocGenerics_0.44.0  ggplot2_3.4.2        DT_0.28             
## 
## loaded via a namespace (and not attached):
##   [1] utf8_1.2.3                              R.utils_2.12.2                          tidyselect_1.2.0                        RSQLite_2.3.1                          
##   [5] htmlwidgets_1.6.2                       grid_4.2.1                              BiocParallel_1.32.6                     methylSig_1.10.0                       
##   [9] scatterpie_0.1.8                        munsell_0.5.0                           codetools_0.2-19                        withr_2.5.0                            
##  [13] colorspace_2.1-0                        GOSemSim_2.24.0                         filelock_1.0.2                          highr_0.10                             
##  [17] knitr_1.43                              rstudioapi_0.15.0                       DOSE_3.24.2                             MatrixGenerics_1.10.0                  
##  [21] labeling_0.4.2                          GenomeInfoDbData_1.2.9                  polyclip_1.10-4                         bit64_4.0.5                            
##  [25] farver_2.1.1                            rhdf5_2.42.1                            vctrs_0.6.3                             treeio_1.23.1                          
##  [29] generics_0.1.3                          xfun_0.39                               BiocFileCache_2.6.1                     R6_2.5.1                               
##  [33] GenomeInfoDb_1.34.9                     graphlayouts_0.8.4                      locfit_1.5-9.7                          bitops_1.0-7                           
##  [37] rhdf5filters_1.10.1                     cachem_1.0.8                            fgsea_1.24.0                            gridGraphics_0.5-1                     
##  [41] DelayedArray_0.24.0                     BiocIO_1.8.0                            scales_1.2.1                            bsseq_1.34.0                           
##  [45] ggraph_2.1.0                            enrichplot_1.18.4                       gtable_0.3.3                            tidygraph_1.2.3                        
##  [49] rlang_1.1.1                             splines_4.2.1                           rtracklayer_1.58.0                      lazyeval_0.2.2                         
##  [53] DSS_2.46.0                              BiocManager_1.30.25                     yaml_2.3.7                              reshape2_1.4.4                         
##  [57] GenomicFeatures_1.50.4                  qvalue_2.30.0                           tools_4.2.1                             ggplotify_0.1.0                        
##  [61] gplots_3.1.3                            jquerylib_0.1.4                         Rcpp_1.0.11                             plyr_1.8.8                             
##  [65] sparseMatrixStats_1.10.0                progress_1.2.2                          zlibbioc_1.44.0                         purrr_1.0.1                            
##  [69] RCurl_1.98-1.12                         prettyunits_1.1.1                       cowplot_1.1.1                           SummarizedExperiment_1.28.0            
##  [73] ggrepel_0.9.3                           magrittr_2.0.3                          matrixStats_1.0.0                       hms_1.1.3                              
##  [77] patchwork_1.1.2                         evaluate_0.21                           HDO.db_0.99.1                           XML_3.99-0.14                          
##  [81] compiler_4.2.1                          biomaRt_2.54.1                          tibble_3.2.1                            KernSmooth_2.23-22                     
##  [85] crayon_1.5.2                            shadowtext_0.1.2                        R.oo_1.25.0                             htmltools_0.5.5                        
##  [89] ggfun_0.0.9                             aplot_0.1.10                            DBI_1.1.3                               tweenr_2.0.2                           
##  [93] ChIPseeker_1.34.1                       dbplyr_2.3.3                            MASS_7.3-60                             rappdirs_0.3.3                         
##  [97] boot_1.3-28.1                           Matrix_1.6-0                            permute_0.9-7                           cli_3.6.1                              
## [101] R.methodsS3_1.8.2                       parallel_4.2.1                          igraph_1.5.0                            GenomicRanges_1.50.2                   
## [105] pkgconfig_2.0.3                         TxDb.Hsapiens.UCSC.hg19.knownGene_3.2.2 GenomicAlignments_1.34.1                plotly_4.10.2                          
## [109] xml2_1.3.5                              ggtree_3.6.2                            bslib_0.5.0                             XVector_0.38.0                         
## [113] yulab.utils_0.0.6                       stringr_1.5.0                           digest_0.6.33                           Biostrings_2.66.0                      
## [117] rmarkdown_2.23                          fastmatch_1.1-3                         tidytree_0.4.2                          DelayedMatrixStats_1.20.0              
## [121] restfulr_0.0.15                         curl_5.0.1                              Rsamtools_2.14.0                        gtools_3.9.4                           
## [125] rjson_0.2.21                            lifecycle_1.0.3                         nlme_3.1-162                            jsonlite_1.8.7                         
## [129] Rhdf5lib_1.20.0                         limma_3.54.2                            BSgenome_1.66.3                         fansi_1.0.4                            
## [133] pillar_1.9.0                            lattice_0.21-8                          plotrix_3.8-2                           KEGGREST_1.38.0                        
## [137] fastmap_1.1.1                           httr_1.4.6                              glue_1.6.2                              png_0.1-8                              
## [141] bit_4.0.5                               ggforce_0.4.1                           stringi_1.7.12                          sass_0.4.7                             
## [145] HDF5Array_1.26.0                        blob_1.2.4                              caTools_1.18.2                          memoise_2.0.1                          
## [149] renv_1.0.9                              ape_5.7-1