Enrichment after FC Scanner Analysis

1. load libraries

2. Enrichment Analysis_part1

# Load necessary libraries ------------------------------------


# Read the TSV file into R
Significant_Genes <- read_csv("/home/bioinfo/18-Enrichment_Analysis_final_Results/NewFiles/All_Results_Together_for_Enrichment/Cell_lines_vs_control/filtered_data_Cell_lines_vs_CD4Tcells.csv")

Rows: 746 Columns: 9── Column specification ────────────────────────────────────────────────────────────────────────────────────────────
Delimiter: ","
chr (1): gene
dbl (8): mean_All_cell_lines_cells, mean_PBMC_CD4T_cells, Relative_variance_All_cell_lines_cells, Relative_varia...
ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.

SS_up  <- subset(Significant_Genes,
                    Significant_Genes$log2FC > 3
                      & Significant_Genes$pValueBH <  0.05)
SS_up_genes <- SS_up$gene


SS_down  <- subset(Significant_Genes,
                    Significant_Genes$log2FC < -1 
                      & Significant_Genes$pValueBH <  0.05)
SS_down_genes <- SS_down$gene


SS_vs_norm_go <- clusterProfiler::enrichGO(SS_up_genes,
                                            "org.Hs.eg.db",
                                            keyType = "SYMBOL",
                                            ont = "BP",
                                            minGSSize = 50)

print(SS_vs_norm_go@result)

enr_go <- clusterProfiler::simplify(SS_vs_norm_go)

print(enr_go@result)

enrichplot::emapplot(enrichplot::pairwise_termsim(enr_go),
                     showCategory = 50)

SS_vs_norm_go <- clusterProfiler::enrichGO(SS_down_genes,
                                            "org.Hs.eg.db",
                                            keyType = "SYMBOL",
                                            ont = "BP",
                                            minGSSize = 50)

print(SS_vs_norm_go@result)

enr_go <- clusterProfiler::simplify(SS_vs_norm_go)

print(enr_go@result)

enrichplot::emapplot(enrichplot::pairwise_termsim(enr_go),
                     showCategory = 50)

NA
NA

3. Enrichment Analysis_part2

gmt <- msigdbr::msigdbr(species = "human", category = "H")

SS_vs_norm_enrich_up <- clusterProfiler::enricher(gene = SS_up_genes,
                                                pAdjustMethod = "BH",
                                                pvalueCutoff  = 0.05,
                                                qvalueCutoff  = 0.05,
                                                TERM2GENE = gmt[,c("gs_name", "gene_symbol")])


print(SS_vs_norm_enrich_up@result[SS_vs_norm_enrich_up@result$p.adjust < 0.05,])

SS_vs_norm_enrich_down <- clusterProfiler::enricher(gene = SS_down_genes,
                                                pAdjustMethod = "BH",
                                                pvalueCutoff  = 0.05,
                                                qvalueCutoff  = 0.05,
                                                TERM2GENE = gmt[,c("gs_name", "gene_symbol")])


print(SS_vs_norm_enrich_down@result[SS_vs_norm_enrich_down@result$p.adjust < 0.05,])

NA

4. Enrichment Analysis_part3

# We use the function `barplot` from package `enrichplot`
graphics::barplot(
  # We provide the variable pointing to enrichment results
  height = SS_vs_norm_enrich_up,
  # We display the best 15 results
  showCategory=15
)

enrichplot::dotplot(
  # We provide the variable pointing to enrichment results
  object = SS_vs_norm_enrich_up,
  # We display the best 15 results
  showCategory=15
)

# We use the function `barplot` from package `enrichplot`
graphics::barplot(
  # We provide the variable pointing to enrichment results
  height = SS_vs_norm_enrich_down,
  # We display the best 15 results
  showCategory=15
)

enrichplot::dotplot(
  # We provide the variable pointing to enrichment results
  object = SS_vs_norm_enrich_down,
  # We display the best 15 results
  showCategory=15
)

5. Enrichment Analysis_part5

# We use the function `heatplot` from `enrichplot` package
# We use the function `upsetplot` from `enrichplot` package
enrichplot::upsetplot(
  # We probide the variable pointing to GSEA results
  x = SS_vs_norm_enrich_up,
  # We show the 10 best results
  n = 10
)

enrichplot::upsetplot(
  # We probide the variable pointing to GSEA results
  x = SS_vs_norm_enrich_down,
  # We show the 10 best results
  n = 10
)

6. Enrichment Analysis_part6

# Load the packages
library(clusterProfiler)

clusterProfiler v4.14.4 Learn more at https://yulab-smu.top/contribution-knowledge-mining/

Please cite:

T Wu, E Hu, S Xu, M Chen, P Guo, Z Dai, T Feng, L Zhou, W Tang, L Zhan, X Fu, S Liu, X Bo, and G Yu.
clusterProfiler 4.0: A universal enrichment tool for interpreting omics data. The Innovation. 2021,
2(3):100141

Attaching package: ‘clusterProfiler’

The following object is masked from ‘package:stats’:

    filter

library(ReactomePA)

ReactomePA v1.50.0 Learn more at https://yulab-smu.top/contribution-knowledge-mining/

Please cite:

Guangchuang Yu, Qing-Yu He. ReactomePA: an R/Bioconductor package for reactome pathway analysis and
visualization. Molecular BioSystems. 2016, 12(2):477-479

library(org.Hs.eg.db)

Loading required package: AnnotationDbi
Loading required package: stats4
Loading required package: BiocGenerics

Attaching package: ‘BiocGenerics’

The following objects are masked from ‘package:dplyr’:

    combine, intersect, setdiff, union

The following object is masked from ‘package:SeuratObject’:

    intersect

The following objects are masked from ‘package:stats’:

    IQR, mad, sd, var, xtabs

The following objects are masked from ‘package:base’:

    anyDuplicated, aperm, append, as.data.frame, basename, cbind, colnames, dirname, do.call,
    duplicated, eval, evalq, Filter, Find, get, grep, grepl, intersect, is.unsorted, lapply, Map,
    mapply, match, mget, order, paste, pmax, pmax.int, pmin, pmin.int, Position, rank, rbind, Reduce,
    rownames, sapply, saveRDS, setdiff, table, tapply, union, unique, unsplit, which.max, which.min

Loading required package: Biobase
Welcome to Bioconductor

    Vignettes contain introductory material; view with 'browseVignettes()'. To cite Bioconductor, see
    'citation("Biobase")', and for packages 'citation("pkgname")'.

Loading required package: IRanges
Loading required package: S4Vectors

Attaching package: ‘S4Vectors’

The following object is masked from ‘package:clusterProfiler’:

    rename

The following objects are masked from ‘package:dplyr’:

    first, rename

The following object is masked from ‘package:utils’:

    findMatches

The following objects are masked from ‘package:base’:

    expand.grid, I, unname


Attaching package: ‘IRanges’

The following object is masked from ‘package:clusterProfiler’:

    slice

The following objects are masked from ‘package:dplyr’:

    collapse, desc, slice

The following object is masked from ‘package:sp’:

    %over%


Attaching package: ‘AnnotationDbi’

The following object is masked from ‘package:clusterProfiler’:

    select

The following object is masked from ‘package:dplyr’:

    select

library(pathview)

##############################################################################
Pathview is an open source software package distributed under GNU General
Public License version 3 (GPLv3). Details of GPLv3 is available at
http://www.gnu.org/licenses/gpl-3.0.html. Particullary, users are required to
formally cite the original Pathview paper (not just mention it) in publications
or products. For details, do citation("pathview") within R.

The pathview downloads and uses KEGG data. Non-academic uses may require a KEGG
license agreement (details at http://www.kegg.jp/kegg/legal.html).
##############################################################################

# Extract the gene symbols and log2FC values for SS_up
gene_list_up <- as.character(SS_up$gene)
log2fc_values_up <- SS_up$log2FC

# Get a list of valid gene symbols from org.Hs.eg.db
valid_symbols <- keys(org.Hs.eg.db, keytype = "SYMBOL")

# Check which gene symbols are not valid
invalid_symbols_up <- setdiff(gene_list_up, valid_symbols)

# Print invalid gene symbols
if (length(invalid_symbols_up) > 0) {
    cat("The following gene symbols are not valid and will be excluded from SS_up:\n")
    print(invalid_symbols_up)
}

The following gene symbols are not valid and will be excluded from SS_up:
 [1] "C19orf48"    "H2AFX"       "WARS"        "VARS"        "WDR34"       "HLA.DMA"     "PALM2.AKAP2" "FAM126A"    
 [9] "AC016747.1"  "CD3EAP"      "AL441992.1"  "C3orf14"     "HLA.DMB"     "ARNTL2"      "AC006064.4"  "AC011603.2" 
[17] "BHLHE40.AS1" "AC093865.1"  "HLA.DPB1"    "HIST1H2AL"   "MSC.AS1"     "HLA.DRB5"    "HLA.DQB1"    "AL590550.1" 
[25] "HIST1H3B"    "AC010967.1"  "HLA.DPA1"    "HIST1H1A"    "HLA.DRA"     "TRBV23.1"    "HLA.DQA1"    "HLA.DRB1"   
[33] "HIST1H1B"    "SLC7A11.AS1" "HLA.DQA2"    "TRAV38.2DV8"

# Filter out invalid gene symbols
gene_list_up <- intersect(gene_list_up, valid_symbols)

# Convert gene symbols to Entrez IDs
entrez_ids_up <- bitr(gene_list_up, fromType = "SYMBOL", toType = "ENTREZID", OrgDb = org.Hs.eg.db)

'select()' returned 1:1 mapping between keys and columns

# Check for successful conversion
if (nrow(entrez_ids_up) == 0) {
    stop("No gene symbols were converted to Entrez IDs for SS_up.")
}

# Merge Entrez IDs with log2FC values
entrez_ids_up <- merge(entrez_ids_up, SS_up, by.x = "SYMBOL", by.y = "gene")

# Create a named vector of log2FC values with Entrez IDs as names
gene_data_up <- setNames(entrez_ids_up$log2FC, entrez_ids_up$ENTREZID)

7. Enrichment Analysis_part7

# Perform GO enrichment analysis
ego <- enrichGO(gene = entrez_ids_up$ENTREZID,
                OrgDb = org.Hs.eg.db,
                keyType = "ENTREZID",
                ont = "BP",  # "BP" for Biological Process
                pAdjustMethod = "BH",
                pvalueCutoff = 0.05,
                qvalueCutoff = 0.2)

# Perform KEGG pathway enrichment analysis
ekegg <- enrichKEGG(gene = entrez_ids_up$ENTREZID,
                    organism = "hsa",
                    pvalueCutoff = 0.001)


# Perform Reactome pathway enrichment analysis
reactome <- enrichPathway(gene = entrez_ids_up$ENTREZID,
                          organism = "human",
                          pvalueCutoff = 0.05)


# Visualize GO enrichment results
barplot(ego, showCategory = 10)

dotplot(ego, showCategory = 10)

# Visualize KEGG pathway enrichment results
barplot(ekegg, showCategory = 30)

dotplot(ekegg, showCategory = 30)

# Visualize Reactome pathway enrichment results
barplot(reactome, showCategory = 10)

dotplot(reactome, showCategory = 10)

# Pathway visualization with pathview
# Example: Visualize KEGG pathway hsa04010 (replace with your pathway of interest)
pathway_up <- pathview(gene.data = as.numeric(entrez_ids_up$ENTREZID),
         pathway.id = "hsa05200", # Example KEGG pathway ID
         species = "hsa")

Warning: None of the genes or compounds mapped to the pathway!
Argument gene.idtype or cpd.idtype may be wrong.
'select()' returned 1:1 mapping between keys and columns
Info: Working in directory /home/bioinfo/17-SingleCellFCscanner
Info: Writing image file hsa05200.pathview.png

pathway_up

$plot.data.gene

$plot.data.cpd

NA

8. Enrichment Analysis_part8

# Load the packages
library(clusterProfiler)
library(ReactomePA)
library(org.Hs.eg.db)
library(pathview)

# Extract the gene symbols and log2FC values for SS_down
gene_list_down <- as.character(SS_down$gene)
log2fc_values_down <- SS_down$log2FC

# Check which gene symbols are not valid
invalid_symbols_down <- setdiff(gene_list_down, valid_symbols)

# Print invalid gene symbols
if (length(invalid_symbols_down) > 0) {
    cat("The following gene symbols are not valid and will be excluded from SS_down:\n")
    print(invalid_symbols_down)
}

The following gene symbols are not valid and will be excluded from SS_down:
[1] "TMEM8A"     "PCED1B.AS1" "AC243960.1" "AC139720.1" "AC119396.1"

# Filter out invalid gene symbols
gene_list_down <- intersect(gene_list_down, valid_symbols)

# Convert gene symbols to Entrez IDs
entrez_ids_down <- bitr(gene_list_down, fromType = "SYMBOL", toType = "ENTREZID", OrgDb = org.Hs.eg.db)

'select()' returned 1:1 mapping between keys and columns

# Check for successful conversion
if (nrow(entrez_ids_down) == 0) {
    stop("No gene symbols were converted to Entrez IDs for SS_down.")
}

# Merge Entrez IDs with log2FC values
entrez_ids_down <- merge(entrez_ids_down, SS_down, by.x = "SYMBOL", by.y = "gene")

# Create a named vector of log2FC values with Entrez IDs as names
gene_data_down <- setNames(entrez_ids_down$log2FC, entrez_ids_down$ENTREZID)

9. Enrichment Analysis_part9

# Install necessary packages if they are not already installed
if (!requireNamespace("BiocManager", quietly = TRUE)) install.packages("BiocManager")
BiocManager::install(c("pathview", "IRdisplay"))

'getOption("repos")' replaces Bioconductor standard repositories, see 'help("repositories", package =
"BiocManager")' for details.
Replacement repositories:
    CRAN: https://cloud.r-project.org
Bioconductor version 3.20 (BiocManager 1.30.25), R 4.4.2 (2024-10-31)
Warning: package(s) not installed when version(s) same as or greater than current; use `force = TRUE` to re-install:
  'pathview'Installing package(s) 'IRdisplay'
also installing the dependency ‘repr’

trying URL 'https://cloud.r-project.org/src/contrib/repr_1.1.7.tar.gz'
Content type 'application/x-gzip' length 32736 bytes (31 KB)
==================================================
downloaded 31 KB

trying URL 'https://cloud.r-project.org/src/contrib/IRdisplay_1.1.tar.gz'
Content type 'application/x-gzip' length 8286 bytes
==================================================
downloaded 8286 bytes

* installing *source* package ‘repr’ ...
** package ‘repr’ successfully unpacked and MD5 sums checked
** using staged installation
** R
** byte-compile and prepare package for lazy loading
** help
*** installing help indices
** building package indices
** testing if installed package can be loaded from temporary location
** testing if installed package can be loaded from final location
** testing if installed package keeps a record of temporary installation path
* DONE (repr)
* installing *source* package ‘IRdisplay’ ...
** package ‘IRdisplay’ successfully unpacked and MD5 sums checked
** using staged installation
** R
** byte-compile and prepare package for lazy loading
** help
*** installing help indices
** building package indices
** testing if installed package can be loaded from temporary location
** testing if installed package can be loaded from final location
** testing if installed package keeps a record of temporary installation path
* DONE (IRdisplay)

The downloaded source packages are in
    ‘/tmp/RtmpzBOl8G/downloaded_packages’
Installation paths not writeable, unable to update packages
  path: /usr/local/lib/R/site-library
  packages:
    BiocVersion
Old packages: 'argparse'

Update all/some/none? [a/s/n]: 

a

trying URL 'https://cloud.r-project.org/src/contrib/argparse_2.2.5.tar.gz'
Content type 'application/x-gzip' length 36066 bytes (35 KB)
==================================================
downloaded 35 KB

* installing *source* package ‘argparse’ ...
** package ‘argparse’ successfully unpacked and MD5 sums checked
** using staged installation
** R
** exec
** inst
** byte-compile and prepare package for lazy loading
** help
*** installing help indices
*** copying figures
** building package indices
** installing vignettes
** testing if installed package can be loaded from temporary location
** testing if installed package can be loaded from final location
** testing if installed package keeps a record of temporary installation path
* DONE (argparse)

The downloaded source packages are in
    ‘/tmp/RtmpzBOl8G/downloaded_packages’

# Load the required packages
library(pathview)
library(IRdisplay)


# Perform GO enrichment analysis
ego_down <- enrichGO(gene = entrez_ids_down$ENTREZID,
                     OrgDb = org.Hs.eg.db,
                     keyType = "ENTREZID",
                     ont = "BP",  # "BP" for Biological Process
                     pAdjustMethod = "BH",
                     pvalueCutoff = 0.05,
                     qvalueCutoff = 0.2)

ekegg_down <- enrichKEGG(
  gene         = entrez_ids_down$ENTREZID,  # Your list of gene IDs
  organism     = "hsa",                # Human pathways
  pvalueCutoff = 0.05,                 # Adjust if needed
  qvalueCutoff = 0.2
)

# Perform Reactome pathway enrichment analysis
reactome_down <- enrichPathway(gene = entrez_ids_down$ENTREZID,
                               organism = "human",
                               pvalueCutoff = 0.05)

# Visualize GO enrichment results
barplot(ego_down, showCategory = 10)

dotplot(ego_down, showCategory = 10)

# Plot if the object is valid
if (!is.null(ekegg_down) && nrow(as.data.frame(ekegg_down)) > 0) {
  barplot(ekegg_down, showCategory = 30)
} else {
  print("No enriched categories found.")
}

dotplot(ekegg_down, showCategory = 30)

if (!is.null(reactome_down) && nrow(as.data.frame(reactome_down)) > 0) {
  barplot(reactome_down, showCategory = 30)
} else {
  print("No enriched categories found.")
}

dotplot(reactome_down, showCategory = 10)

# Pathway visualization with pathview for downregulated genes
pathway_down <- pathview(gene.data = gene_data_down,
         pathway.id = "hsa05200", # Example KEGG pathway ID
         species = "hsa",
         kegg.native = TRUE) # Use KEGG's native pathway diagram

'select()' returned 1:1 mapping between keys and columns
Info: Working in directory /home/bioinfo/17-SingleCellFCscanner
Info: Writing image file hsa05200.pathview.png

pathway_down

$plot.data.gene

$plot.data.cpd

NA