Enrichment after FC Scanner Analysis-part2

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/2-filtered_All_cell_lines_cells_vs_CD4_Control.csv")

Rows: 3018 Columns: 9── Column specification ────────────────────────────────────────────────────────────────────────────────────────────────────
Delimiter: ","
chr (1): gene
dbl (8): mean.All_cell_lines_cells, mean..CD4_Control, Relative.variance.All_cell_lines_cells, Relative.variance.CD4_Con...
ℹ 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)
library(ReactomePA)
library(org.Hs.eg.db)
library(pathview)

# 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)

Reading KEGG annotation online: "https://rest.kegg.jp/link/hsa/pathway"...
Reading KEGG annotation online: "https://rest.kegg.jp/list/pathway/hsa"...

# 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/18-Enrichment_Analysis_final_Results
Info: Writing image file hsa05200.pathview.png

pathway_up

$plot.data.gene

$plot.data.cpd

# Assuming the image is saved with the name 'hsa05200.pathview.png'
image_path <- "hsa05200.hsa05200.png"

# Check if the image file exists and display it
if (file.exists(image_path)) {
    display_png(file = image_path)
} else {
    cat("Image file", image_path, "not found.")
}

Image file hsa05200.hsa05200.png not found.

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" "SATB1.AS1"  "PRKCQ.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"))

# 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.")
}

[1] "No enriched categories found."

dotplot(reactome_down, showCategory = 10)

Error in ans[ypos] <- rep(yes, length.out = len)[ypos] : 
  replacement has length zero

10. Enrichment Analysis_part10

11. Enrichment Analysis_part11

12. Enrichment Analysis_part12

13. Enrichment Analysis_part13