Introduction

Ewing sarcoma is a pediatric bone cancer driven by the EWSR1-FLI1 fusion oncogene. This analysis examines the transcriptomic changes following EWSR1-FLI1 knockdown to identify differentially expressed genes and affected biological pathways.

Experimental Design

Control group (shCTR): 3 samples
EWSR1-FLI1 knockdown (shEF1): 4 samples
Total samples: 7

Setup

# Data manipulation and visualization
library(tidyverse)
library(scales)
library(dplyr)
library(knitr)

# RNA-seq analysis
library(DESeq2)
library(SummarizedExperiment)

# Visualization
library(EnhancedVolcano)
library(pheatmap)
library(RColorBrewer)

# Enrichment analysis
library(clusterProfiler)
library(enrichplot)
library(msigdbr)
library(org.Hs.eg.db)
library(AnnotationDbi)
library(biomaRt)

# Table formatting
library(kableExtra)
library(DT)

# Clean Ensembl IDs (remove version numbers)
ensembl_ids <- rownames(rse)
ensembl_ids_clean <- gsub("\\..*", "", ensembl_ids)

cat("Starting gene annotation...\n\n")

## Starting gene annotation...

# Step 1: Quick annotation with org.Hs.eg.db
gene_symbols_orgdb <- mapIds(org.Hs.eg.db,
                              keys = ensembl_ids_clean,
                              column = "SYMBOL",
                              keytype = "ENSEMBL",
                              multiVals = "first")

orgdb_annotated <- sum(!is.na(gene_symbols_orgdb))
cat("  ✓", orgdb_annotated, "genes annotated\n\n")

##   ✓ 35480 genes annotated

ensembl <- useEnsembl(biomart = "genes", 
                      dataset = "hsapiens_gene_ensembl",
                      mirror = "useast")

# Query in chunks to avoid timeout
chunk_size <- 5000
n_chunks <- ceiling(length(ensembl_ids_clean) / chunk_size)
biomart_results <- list()

for (i in 1:n_chunks) {
  start_idx <- (i - 1) * chunk_size + 1
  end_idx <- min(i * chunk_size, length(ensembl_ids_clean))
  chunk_ids <- ensembl_ids_clean[start_idx:end_idx]
  
  cat("  Chunk", i, "/", n_chunks, "...")
  
  tryCatch({
    chunk_result <- getBM(
      attributes = c('ensembl_gene_id', 'external_gene_name', 'gene_biotype'),
      filters = 'ensembl_gene_id',
      values = chunk_ids,
      mart = ensembl
    )
    biomart_results[[i]] <- chunk_result
    cat(" ✓\n")
    Sys.sleep(0.5)
  }, error = function(e) {
    cat(" ERROR\n")
    biomart_results[[i]] <- NULL
  })
}

##   Chunk 1 / 12 ... ✓
##   Chunk 2 / 12 ... ✓
##   Chunk 3 / 12 ... ✓
##   Chunk 4 / 12 ... ✓
##   Chunk 5 / 12 ... ✓
##   Chunk 6 / 12 ... ✓
##   Chunk 7 / 12 ... ✓
##   Chunk 8 / 12 ... ✓
##   Chunk 9 / 12 ... ✓
##   Chunk 10 / 12 ... ✓
##   Chunk 11 / 12 ... ✓
##   Chunk 12 / 12 ... ✓

biomart_annotation <- bind_rows(biomart_results)
cat("\n  ✓", nrow(biomart_annotation), "genes retrieved from biomaRt\n\n")

## 
##   ✓ 56062 genes retrieved from biomaRt

# Merge annotations: org.Hs.eg.db, fill gaps with biomaRt
gene_annotation <- tibble(
  ensembl_id = ensembl_ids,
  ensembl_id_clean = ensembl_ids_clean,
  gene_symbol_orgdb = gene_symbols_orgdb
) %>%
  left_join(biomart_annotation, by = c("ensembl_id_clean" = "ensembl_gene_id")) %>%
  mutate(
    gene_symbol = coalesce(gene_symbol_orgdb, external_gene_name),
    display_name = coalesce(gene_symbol, ensembl_id_clean)
  ) %>%
  dplyr::select(ensembl_id, ensembl_id_clean, gene_symbol, gene_biotype, display_name)

# Summary
total_annotated <- sum(!is.na(gene_annotation$gene_symbol))
pct_annotated <- round(100 * total_annotated / nrow(gene_annotation), 1)
improvement <- total_annotated - orgdb_annotated

cat("=== ANNOTATION SUMMARY ===\n")

## === ANNOTATION SUMMARY ===

cat("  Total genes:", nrow(gene_annotation), "\n")

##   Total genes: 58037

cat("  Annotated:", total_annotated, "(", pct_annotated, "%)\n")

##   Annotated: 56261 ( 96.9 %)

cat("  Improvement:", improvement, "additional genes from biomaRt\n\n")

##   Improvement: 20781 additional genes from biomaRt

Exploratory Data Analysis

Sample Overview

sample_info %>%
   dplyr::select(sample_id, condition) %>%
  kable(
    caption = "Sample Information",
    col.names = c("Sample ID", "Condition")
  ) %>%
  kable_styling(
    bootstrap_options = c("striped", "hover", "condensed"),
    full_width = FALSE
  )

Sample Information
Sample ID	Condition
SRR579382	shCTR
SRR579383	shCTR
SRR579384	shCTR
SRR579385	shEF1
SRR579386	shEF1
SRR579387	shEF1
SRR579513	shEF1

We have 3 control samples and 4 EWSR1-FLI1 knockdown samples.

Principal Component Analysis

PCA reveals the overall structure of the data and potential batch effects.

# Create DESeq2 object
dds <- DESeqDataSet(rse, design = ~condition)

cat("Sample sizes:\n")

## Sample sizes:

table(sample_info$condition)

## 
## shCTR shEF1 
##     3     4

# Determine smallest group size
smallest_group_size <- min(table(sample_info$condition))
cat("\nSmallest group size:", smallest_group_size, "\n")

## 
## Smallest group size: 3

# Pre-filtering: remove genes with very low counts
# At least 10 counts in smallest group
smallest_group_size <- min(table(colData(dds)$condition))
keep <- rowSums(counts(dds) >= 10) >= smallest_group_size
dds <- dds[keep, ]
cat("Filtering: retained", sum(keep), "genes\n")

## Filtering: retained 31504 genes

# Update annotation to match filtered genes
gene_annotation_filtered <- gene_annotation %>%
  filter(ensembl_id %in% rownames(dds))

# Variance stabilizing transformation for PCA
vsd <- vst(dds, blind = TRUE)

# Perform PCA
pca_data <- plotPCA(vsd, intgroup = "condition", returnData = TRUE)
percent_var <- round(100 * attr(pca_data, "percentVar"), 1)

# Enhanced PCA plot
ggplot(pca_data, aes(x = PC1, y = PC2, color = condition, label = name)) +
  geom_point(size = 5, alpha = 0.8) +
  geom_text(hjust = 0, vjust = -0.5, size = 3.5, show.legend = FALSE) +
  scale_color_manual(
    values = c("shCTR" = "#2E86AB", "shEF1" = "#A23B72"),
    labels = c("shCTR" = "Control", "shEF1" = "EWSR1-FLI1 KD")
  ) +
  labs(
    title = "Principal Component Analysis",
    subtitle = "Samples cluster by EWSR1-FLI1 knockdown status",
    x = paste0("PC1: ", percent_var[1], "% variance"),
    y = paste0("PC2: ", percent_var[2], "% variance"),
    color = "Condition"
  ) +
  theme_minimal(base_size = 14) +
  theme(
    plot.title = element_text(face = "bold", size = 16),
    legend.position = "bottom"
  )

PCA plot showing sample clustering by condition

The first principal component explains 98.3% of the variance and clearly separates control from knockdown samples, indicating a strong transcriptomic response to EWSR1-FLI1 suppression.

Differential Gene Expression Analysis

DESeq2 Analysis

# Set reference level
dds$condition <- relevel(dds$condition, ref = "shCTR")

# Run DESeq2
dds <- DESeq(dds)

# Extract results
res_original <- results(dds, contrast = c("condition", "shEF1", "shCTR"))

# Create summary
cat("DESeq2 Results Summary - EF1 vs CTR:\n")

## DESeq2 Results Summary - EF1 vs CTR:

print(summary(res_original))

## 
## out of 31504 with nonzero total read count
## adjusted p-value < 0.1
## LFC > 0 (up)       : 8202, 26%
## LFC < 0 (down)     : 7130, 23%
## outliers [1]       : 250, 0.79%
## low counts [2]     : 0, 0%
## (mean count < 6)
## [1] see 'cooksCutoff' argument of ?results
## [2] see 'independentFiltering' argument of ?results
## 
## NULL

# Apply LFC shrinkage for more accurate effect sizes
library(apeglm)
res_shrink <- lfcShrink(dds, 
                       coef = "condition_shEF1_vs_shCTR",
                       type = "apeglm")

# Convert to tibble for tidyverse manipulation
res_original_tbl <- res_original %>%
  as.data.frame() %>%
  rownames_to_column("ensembl_id") %>%
  as_tibble() %>%
  left_join(gene_annotation_filtered, by = "ensembl_id") %>%  # Join annotation here!
  arrange(padj) %>%
  mutate(significance = case_when(
      is.na(padj) ~ "Not significant",
      padj < 0.05 & log2FoldChange > 1 ~ "Up-regulated (LFC>1)",
      padj < 0.05 & log2FoldChange > 0.5 ~ "Up-regulated (0.5<LFC<1)",
      padj < 0.05 & log2FoldChange < -1 ~ "Down-regulated (LFC<-1)",
      padj < 0.05 & log2FoldChange < -0.5 ~ "Down-regulated (-1<LFC<-0.5)",
      padj < 0.05 ~ "Significant (|LFC| < 0.5)",
      TRUE ~ "Not significant"
    )
  )
    
res_shrink_tbl <- res_shrink %>%
  as.data.frame() %>%
   rownames_to_column("ensembl_id") %>%
  as_tibble() %>%
  left_join(gene_annotation_filtered, by = "ensembl_id") %>%  # Join annotation here!
  arrange(padj) %>%
  mutate(significance = case_when(
      is.na(padj) ~ "Not significant",
      padj < 0.05 & log2FoldChange > 1 ~ "Up-regulated (LFC>1)",
      padj < 0.05 & log2FoldChange > 0.5 ~ "Up-regulated (0.5<LFC<1)",
      padj < 0.05 & log2FoldChange < -1 ~ "Down-regulated (LFC<-1)",
      padj < 0.05 & log2FoldChange < -0.5 ~ "Down-regulated (-1<LFC<-0.5)",
      padj < 0.05 ~ "Significant (|LFC| < 0.5)",
      TRUE ~ "Not significant"
    )
  )

# Summary statistics
n_up <- sum(res_original_tbl$padj < 0.05 & res_original_tbl$log2FoldChange > 1, na.rm = TRUE)
n_down <- sum(res_original_tbl$padj < 0.05 & res_original_tbl$log2FoldChange < -1, na.rm = TRUE)
n_sig <- sum(res_original_tbl$padj < 0.05, na.rm = TRUE)

# Summary statistics
n_up_shrink <- sum(res_shrink_tbl$padj < 0.05 & res_shrink_tbl$log2FoldChange > 1, na.rm = TRUE)
n_down_shrink <- sum(res_shrink_tbl$padj < 0.05 & res_shrink_tbl$log2FoldChange < -1, na.rm = TRUE)
n_sig_shrink <- sum(res_shrink_tbl$padj < 0.05, na.rm = TRUE)

Summary

Total genes tested: 31504
Significant genes (padj < 0.05): 13979
Up-regulated (LFC > 1, padj < 0.05): 4278
Down-regulated (LFC < -1, padj < 0.05): 1947

MA Plot

res_shrink_tbl %>%
  mutate(
    sig_color = case_when(
      significance == "Up-regulated (LFC>1)" ~ "Up-regulated",
      significance == "Down-regulated (LFC<-1)" ~ "Down-regulated",
      TRUE ~ "Not significant"
    )
  ) %>%
  ggplot(aes(x = baseMean, y = log2FoldChange)) +
  geom_point(aes(color = sig_color), alpha = 0.5, size = 1) +
  scale_x_log10(
    labels = label_comma(),
    breaks = c(1, 10, 100, 1000, 10000, 100000)
  ) +
  scale_color_manual(
    values = c(
      "Up-regulated" = "#D62828",
      "Down-regulated" = "#003049",
      "Not significant" = "grey70"
    )
  ) +
  geom_hline(yintercept = 0, linetype = "dashed", color = "black") +
  labs(
    title = "MA Plot: Mean Expression vs. Log2 Fold Change",
    subtitle = paste0(n_sig_shrink, " genes significantly altered by EWSR1-FLI1 knockdown"),
    x = "Mean of Normalized Counts",
    y = "Log2 Fold Change (shEF1 / shCTR)",
    color = "Regulation"
  ) +
  theme_minimal(base_size = 14) +
  theme(
    plot.title = element_text(face = "bold", size = 16),
    legend.position = "bottom"
  )

MA plot showing the relationship between mean expression and log2 fold change

Differentially Expressed Genes

res_original_tbl %>%
  dplyr::filter(padj < 0.05) %>%
  dplyr:: select(gene_symbol, ensembl_id_clean, gene_biotype, baseMean, log2FoldChange, padj) %>%
  dplyr::mutate(
    baseMean = round(baseMean, 2),
    log2FoldChange = round(log2FoldChange, 3),
    padj = format(padj, scientific = TRUE, digits = 5)
  ) %>%
  DT::datatable(
    caption = "Differentially Expressed Genes (FDR < 0.05)",
    colnames = c("Gene", "Ensembl ID", "Biotype", "Mean Expr", "Log2 FC", "Adj p-value"),
    filter = "top",
    options = list(pageLength = 15, scrollX = TRUE)
  )

Volcano Plot

# Get top genes for labeling
top_genes <- res_shrink_tbl %>%
  filter(!is.na(padj), !is.na(gene_symbol)) %>%
  arrange(padj) %>%
  slice_head(n = 15) %>%
  pull(gene_symbol)

# Prepare data with labels
res_volcano <- res_shrink_tbl %>%
 mutate(label = if_else(!is.na(gene_symbol) & gene_symbol %in% top_genes, gene_symbol, ""))

EnhancedVolcano(
  res_shrink_tbl,
  lab = res_volcano$label,
  x = "log2FoldChange",
  y = "padj",
  title = "EWSR1-FLI1 Knockdown vs. Control",
  subtitle = paste0(n_sig_shrink, " genes with FDR < 0.05"),
  pCutoff = 0.05,
  FCcutoff = 1,
  pointSize = 2.0,
  labSize = 4.5,
  selectLab = top_genes,
  col = c("grey30", "#2E86AB", "#A23B72", "#D62828"),
  colAlpha = 0.6,
  legendPosition = "right",
  legendLabSize = 12,
  legendIconSize = 4.0,
  drawConnectors = TRUE,
  widthConnectors = 0.5,
  maxoverlapsConnectors = 20
)

Volcano plot showing statistical significance vs. fold change

Expression Heatmap

# Get top 10 up and down regulated genes
top_genes_list <- res_original_tbl %>%
 filter(!is.na(padj), padj < 0.05, !is.na(gene_symbol)) %>%
  arrange(desc(log2FoldChange)) %>%
  slice_head(n = 10) %>%
  bind_rows(
    res_original_tbl %>%
     filter(!is.na(padj), padj < 0.05, !is.na(gene_symbol)) %>%
      arrange(log2FoldChange) %>%
      slice_head(n = 10)
  ) 

# Extract normalized counts
norm_counts <- counts(dds, normalized = TRUE)[top_genes_list$ensembl_id, ]
heatmap_data <- t(scale(t(norm_counts)))
# Rename rows to gene symbols
rownames(heatmap_data) <- top_genes_list$gene_symbol

# Annotation
annotation_col <- sample_info %>%
  dplyr::select(sample_id, condition) %>%
  tibble::column_to_rownames("sample_id")


annotation_colors <- list(
  condition = c(shCTR = "#2E86AB", shEF1 = "#A23B72")
)

# Create heatmap
pheatmap(
  heatmap_data,
  annotation_col = annotation_col,
  annotation_colors = annotation_colors,
  color = colorRampPalette(rev(brewer.pal(n = 11, name = "RdBu")))(100),
  breaks = seq(-2, 2, length.out = 100),
  cluster_cols = TRUE,
  cluster_rows = TRUE,
  show_rownames = TRUE,
  show_colnames = TRUE,
  fontsize = 10,
  fontsize_row = 9,
  main = "Top 20 Differentially Expressed Genes\n(10 Up-regulated + 10 Down-regulated)",
  border_color = NA
)

Heatmap of top 20 differentially expressed genes

Pathway Enrichment Analysis

# Prepare gene lists with multiple thresholds
genes_up_strict <- res_original_tbl %>%
  filter(padj < 0.05, log2FoldChange > 1, !is.na(gene_symbol)) %>%
  pull(gene_symbol)

genes_down_strict <- res_original_tbl %>%
  filter(padj < 0.05, log2FoldChange < -1, !is.na(gene_symbol)) %>%
  pull(gene_symbol)

genes_up_relaxed <- res_original_tbl %>%
  filter(padj < 0.05, log2FoldChange > 0.5, !is.na(gene_symbol)) %>%
  pull(gene_symbol)

genes_down_relaxed <- res_original_tbl %>%
  filter(padj < 0.05, log2FoldChange < -0.5, !is.na(gene_symbol)) %>%
  pull(gene_symbol)

# Prepare ranked gene list for GSEA (uses all genes)
ranked_genes <- res_original_tbl %>%
  filter(!is.na(padj), !is.na(log2FoldChange), !is.na(gene_symbol)) %>%
  mutate(rank_metric = -log10(padj) * sign(log2FoldChange)) %>%
  arrange(desc(rank_metric)) %>%
  dplyr::select(gene_symbol, rank_metric) %>%
  deframe()

# Get KEGG gene sets from msigdbr
kegg_gene_sets <- msigdbr(species = "Homo sapiens", category = "C2", subcategory = "CP:KEGG")
kegg_t2g <- kegg_gene_sets %>%
  dplyr::select(gs_name, gene_symbol) %>%
  distinct()

# Load gene set databases
hallmark_sets <- msigdbr(species = "Homo sapiens", category = "H")
hallmark_t2g <- hallmark_sets %>%
  dplyr::select(gs_name, gene_symbol) %>%
  distinct()

go_bp_sets <- msigdbr(species = "Homo sapiens", category = "C5", subcategory = "GO:BP")
go_bp_t2g <- go_bp_sets %>%
  dplyr::select(gs_name, gene_symbol) %>%
  distinct()

reactome_sets <- msigdbr(species = "Homo sapiens", category = "C2", subcategory = "CP:REACTOME")
reactome_t2g <- reactome_sets %>%
  dplyr::select(gs_name, gene_symbol) %>%
  distinct()

cat("Gene lists prepared:\n")

## Gene lists prepared:

cat("  Up-regulated (LFC > 1):", length(genes_up_strict), "genes\n")

##   Up-regulated (LFC > 1): 4172 genes

cat("  Up-regulated (LFC > 0.5):", length(genes_up_relaxed), "genes\n")

##   Up-regulated (LFC > 0.5): 6001 genes

cat("  Down-regulated (LFC < -1):", length(genes_down_strict), "genes\n")

##   Down-regulated (LFC < -1): 1905 genes

cat("  Down-regulated (LFC < -0.5):", length(genes_down_relaxed), "genes\n")

##   Down-regulated (LFC < -0.5): 4104 genes

Up-regulated Pathways KEGG

if (length(genes_up_strict) > 0) {
  enrich_up <- enricher(
    gene = genes_up_strict,
    TERM2GENE = kegg_t2g,
    pvalueCutoff = 0.05,
    qvalueCutoff = 0.2
  )
  
  if (!is.null(enrich_up) && nrow(as.data.frame(enrich_up)) > 0) {
    enrich_up_df <- enrich_up %>%
      as.data.frame() %>%
      as_tibble() %>%
      arrange(p.adjust) %>%
      slice_head(n = 10) %>%
      mutate(
        Description = str_remove(Description, "KEGG_"),
        Description = str_replace_all(Description, "_", " ")
      )
    
    # Table
    enrich_up_df %>%
      dplyr::select(Description, GeneRatio, p.adjust, Count) %>%
      mutate(p.adjust = format(p.adjust, scientific = TRUE, digits = 3)) %>%
      kable(
        caption = "Top Up-regulated KEGG Pathways",
        col.names = c("Pathway", "Gene Ratio", "Adjusted p-value", "Gene Count")
      ) %>%
      kable_styling(
        bootstrap_options = c("striped", "hover", "condensed"),
        full_width = FALSE
      )
    
    # Plot
    enrich_up_df %>%
      mutate(
        Description = fct_reorder(Description, -p.adjust),
        neg_log_p = -log10(p.adjust)
      ) %>%
      ggplot(aes(x = neg_log_p, y = Description, fill = Count)) +
      geom_col() +
      scale_fill_gradient(low = "#FFC09F", high = "#D62828") +
      labs(
        title = "Up-regulated KEGG Pathways",
        subtitle = "Top 10 enriched pathways in EWSR1-FLI1 knockdown",
        x = "-log10(Adjusted p-value)",
        y = NULL,
        fill = "Gene\nCount"
      ) +
      theme_minimal(base_size = 14) +
      theme(
        plot.title = element_text(face = "bold", size = 16),
        axis.text.y = element_text(size = 11)
      )
  } else {
    cat("No significantly enriched pathways found for up-regulated genes.\n")
  }
} else {
  cat("No up-regulated genes available for enrichment analysis.\n")
}

Down-regulated Pathways KEGG

if (length(genes_down_strict) > 0) {
  enrich_down <- enricher(
    gene = genes_down_strict,
    TERM2GENE = kegg_t2g,
    pvalueCutoff = 0.05,
    qvalueCutoff = 0.2
  )
  
  if (!is.null(enrich_down) && nrow(as.data.frame(enrich_down)) > 0) {
    enrich_down_df <- enrich_down %>%
      as.data.frame() %>%
      as_tibble() %>%
      arrange(p.adjust) %>%
      slice_head(n = 10) %>%
      mutate(
        Description = str_remove(Description, "KEGG_"),
        Description = str_replace_all(Description, "_", " ")
      )
    
    # Table
    enrich_down_df %>%
      dplyr::select(Description, GeneRatio, p.adjust, Count) %>%
      mutate(p.adjust = format(p.adjust, scientific = TRUE, digits = 3)) %>%
      kable(
        caption = "Top Down-regulated KEGG Pathways",
        col.names = c("Pathway", "Gene Ratio", "Adjusted p-value", "Gene Count")
      ) %>%
      kable_styling(
        bootstrap_options = c("striped", "hover", "condensed"),
        full_width = FALSE
      )
    
    # Plot
    enrich_down_df %>%
      mutate(
        Description = fct_reorder(Description, -p.adjust),
        neg_log_p = -log10(p.adjust)
      ) %>%
      ggplot(aes(x = neg_log_p, y = Description, fill = Count)) +
      geom_col() +
      scale_fill_gradient(low = "#A7C7E7", high = "#003049") +
      labs(
        title = "Down-regulated KEGG Pathways",
        subtitle = "Top 10 enriched pathways in EWSR1-FLI1 knockdown",
        x = "-log10(Adjusted p-value)",
        y = NULL,
        fill = "Gene\nCount"
      ) +
      theme_minimal(base_size = 14) +
      theme(
        plot.title = element_text(face = "bold", size = 16),
        axis.text.y = element_text(size = 11)
      )
  } else {
    cat("No significantly enriched pathways found for down-regulated genes.\n")
  }
} else {
  cat("No down-regulated genes available for enrichment analysis.\n")
}

Hallmark Pathways (GSEA)

Gene Set Enrichment Analysis uses all genes ranked by significance and effect size.

# Check for non-finite values
sum(!is.finite(ranked_genes))  # number of problematic entries

## [1] 12

# Remove non-finite values
ranked_genes_clean <- ranked_genes[is.finite(ranked_genes)]

# Optional: make sure names are still there (for GSEA)
ranked_genes_clean <- ranked_genes_clean[!is.na(names(ranked_genes_clean))]

gsea_hallmark <- GSEA(
  geneList = ranked_genes_clean,
  TERM2GENE = hallmark_t2g,
  pvalueCutoff = 0.25,
  pAdjustMethod = "BH",
  minGSSize = 10,
  maxGSSize = 500,
  verbose = FALSE
)

if (!is.null(gsea_hallmark) && nrow(as.data.frame(gsea_hallmark)) > 0) {
  gsea_df <- gsea_hallmark %>%
    as.data.frame() %>%
    arrange(pvalue) %>%
    mutate(
      Description = str_remove(Description, "HALLMARK_"),
      Description = str_replace_all(Description, "_", " ")
    )
  
  # Table
  gsea_df %>%
    dplyr::select(Description, NES, pvalue, p.adjust, setSize) %>%
    mutate(
      pvalue = format(pvalue, scientific = TRUE, digits = 3),
      p.adjust = format(p.adjust, scientific = TRUE, digits = 3),
      NES = round(NES, 3)
    ) %>%
    slice_head(n = 20) %>%
    kable(
      caption = "GSEA: Hallmark Pathways",
      col.names = c("Pathway", "NES", "p-value", "Adjusted p", "Set Size")
    ) %>%
    kable_styling(
      bootstrap_options = c("striped", "hover", "condensed"),
      full_width = FALSE
    )
  
  # Dotplot
  dotplot(gsea_hallmark, showCategory = 20, font.size = 10) +
    ggtitle("GSEA: Hallmark Pathways") +
    theme(plot.title = element_text(face = "bold", size = 14))
  
  # Enrichment plot for top pathways
  if (nrow(gsea_df) >= 3) {
    gseaplot2(gsea_hallmark, geneSetID = 1:3, 
              pvalue_table = TRUE,
              title = "Top 3 Hallmark Pathways")
  }
} else {
  cat("No significant Hallmark pathways found.\n")
}

ORA: Up-regulated Pathways

if (length(genes_up_strict) > 5) {
  ora_up <- enricher(gene = genes_up_strict, TERM2GENE = hallmark_t2g,
                     pvalueCutoff = 0.1, qvalueCutoff = 0.25, minGSSize = 5)
  
  if (!is.null(ora_up) && nrow(as.data.frame(ora_up)) > 0) {
    ora_df <- ora_up %>%
      as.data.frame() %>%
      mutate(Description = str_remove(Description, "HALLMARK_"),
             Description = str_replace_all(Description, "_", " ")) %>%
      arrange(p.adjust) %>%
      slice_head(n = 15)
    
    ora_df %>%
      dplyr::select(Description, GeneRatio, p.adjust, Count) %>%
      kable(caption = "Up-regulated Pathways") %>%
      kable_styling(bootstrap_options = c("striped", "hover"), full_width = FALSE)
    
    ora_df %>%
      mutate(Description = fct_reorder(Description, -p.adjust)) %>%
      ggplot(aes(x = -log10(p.adjust), y = Description, fill = Count)) +
      geom_col() +
      scale_fill_gradient(low = "#FFC09F", high = "#D62828") +
      labs(title = "Up-regulated Pathways", x = "-log10(Adj p)", y = NULL) +
      theme_minimal(base_size = 12)
  } else {
    cat("No significant up-regulated pathways.\n")
  }
}

ORA: Down-regulated Pathways

if (length(genes_down_strict) > 5) {
  ora_down <- enricher(gene = genes_down_strict, TERM2GENE = hallmark_t2g,
                       pvalueCutoff = 0.1, qvalueCutoff = 0.25, minGSSize = 5)
  
  if (!is.null(ora_down) && nrow(as.data.frame(ora_down)) > 0) {
    ora_df <- ora_down %>%
      as.data.frame() %>%
      mutate(Description = str_remove(Description, "HALLMARK_"),
             Description = str_replace_all(Description, "_", " ")) %>%
      arrange(p.adjust) %>%
      slice_head(n = 15)
    
    ora_df %>%
      dplyr::select(Description, GeneRatio, p.adjust, Count) %>%
      kable(caption = "Down-regulated Pathways") %>%
      kable_styling(bootstrap_options = c("striped", "hover"), full_width = FALSE)
    
    ora_df %>%
      mutate(Description = fct_reorder(Description, -p.adjust)) %>%
      ggplot(aes(x = -log10(p.adjust), y = Description, fill = Count)) +
      geom_col() +
      scale_fill_gradient(low = "#A7C7E7", high = "#003049") +
      labs(title = "Down-regulated Pathways", x = "-log10(Adj p)", y = NULL) +
      theme_minimal(base_size = 12)
  } else {
    cat("No significant down-regulated pathways.\n")
  }
}

Summary

Key Findings

This comprehensive differential gene expression analysis of EWSR1-FLI1 knockdown in Ewing sarcoma reveals:

Strong Transcriptomic Response: PCA analysis shows clear separation between control and knockdown samples, with PC1 explaining 98.3% of variance.
Widespread Gene Expression Changes:
- 13979 genes significantly altered (FDR < 0.05)
- 4278 genes up-regulated (LFC > 1)
- 1947 genes down-regulated (LFC < -1)
Biological Pathways: Pathway enrichment analysis identifies key biological processes affected by EWSR1-FLI1 suppression.

Results Summary Table

tibble(
  Metric = c(
    "Total Genes Tested",
    "Significantly Altered Genes",
    "Up-regulated Genes",
    "Down-regulated Genes",
    "PC1 Variance Explained"
  ),
  Value = c(
    format(nrow(res_original_tbl), big.mark = ","),
    format(n_sig, big.mark = ","),
    format(n_up, big.mark = ","),
    format(n_down, big.mark = ","),
    paste0(percent_var[1], "%")
  )
) %>%
  kable(
    caption = "Summary of Differential Gene Expression Analysis"
  ) %>%
  kable_styling(
    bootstrap_options = c("striped", "hover"),
    full_width = FALSE
  )

Summary of Differential Gene Expression Analysis
Metric	Value
Total Genes Tested	31,504
Significantly Altered Genes	13,979
Up-regulated Genes	4,278
Down-regulated Genes	1,947
PC1 Variance Explained	98.3%

Conclusions

EWSR1-FLI1 knockdown induces substantial transcriptomic changes in Ewing sarcoma cells, affecting thousands of genes and multiple biological pathways. These findings provide insight into the molecular mechanisms underlying EWSR1-FLI1-driven oncogenesis and may inform therapeutic strategies targeting this fusion oncogene.

Session Information

sessionInfo()

## R version 4.4.1 (2024-06-14 ucrt)
## Platform: x86_64-w64-mingw32/x64
## Running under: Windows 11 x64 (build 22631)
## 
## Matrix products: default
## 
## 
## locale:
## [1] LC_COLLATE=English_United States.utf8 
## [2] LC_CTYPE=English_United States.utf8   
## [3] LC_MONETARY=English_United States.utf8
## [4] LC_NUMERIC=C                          
## [5] LC_TIME=English_United States.utf8    
## 
## time zone: America/New_York
## tzcode source: internal
## 
## attached base packages:
## [1] stats4    stats     graphics  grDevices utils     datasets  methods  
## [8] base     
## 
## other attached packages:
##  [1] apeglm_1.26.1               DT_0.33                    
##  [3] kableExtra_1.4.0            biomaRt_2.58.2             
##  [5] org.Hs.eg.db_3.19.1         AnnotationDbi_1.66.0       
##  [7] msigdbr_7.5.1               enrichplot_1.24.4          
##  [9] clusterProfiler_4.12.6      RColorBrewer_1.1-3         
## [11] pheatmap_1.0.12             EnhancedVolcano_1.22.0     
## [13] ggrepel_0.9.6               DESeq2_1.44.0              
## [15] SummarizedExperiment_1.34.0 Biobase_2.64.0             
## [17] MatrixGenerics_1.16.0       matrixStats_1.4.1          
## [19] GenomicRanges_1.56.1        GenomeInfoDb_1.40.1        
## [21] IRanges_2.38.1              S4Vectors_0.42.1           
## [23] BiocGenerics_0.50.0         knitr_1.48                 
## [25] scales_1.4.0                lubridate_1.9.3            
## [27] forcats_1.0.0               stringr_1.5.1              
## [29] dplyr_1.1.4                 purrr_1.0.2                
## [31] readr_2.1.5                 tidyr_1.3.1                
## [33] tibble_3.2.1                ggplot2_4.0.0              
## [35] tidyverse_2.0.0            
## 
## loaded via a namespace (and not attached):
##   [1] splines_4.4.1           ggplotify_0.1.2         filelock_1.0.3         
##   [4] R.oo_1.26.0             polyclip_1.10-7         XML_3.99-0.17          
##   [7] lifecycle_1.0.4         httr2_1.0.5             lattice_0.22-6         
##  [10] MASS_7.3-61             crosstalk_1.2.1         magrittr_2.0.3         
##  [13] sass_0.4.9              rmarkdown_2.28          jquerylib_0.1.4        
##  [16] yaml_2.3.10             cowplot_1.1.3           DBI_1.2.3              
##  [19] abind_1.4-8             zlibbioc_1.50.0         R.utils_2.12.3         
##  [22] ggraph_2.2.1            yulab.utils_0.1.7       tweenr_2.0.3           
##  [25] rappdirs_0.3.3          GenomeInfoDbData_1.2.12 tidytree_0.4.6         
##  [28] svglite_2.1.3           codetools_0.2-20        DelayedArray_0.30.1    
##  [31] DOSE_3.30.5             xml2_1.3.6              ggforce_0.4.2          
##  [34] tidyselect_1.2.1        aplot_0.2.3             UCSC.utils_1.0.0       
##  [37] farver_2.1.2            viridis_0.6.5           BiocFileCache_2.12.0   
##  [40] jsonlite_1.8.9          tidygraph_1.3.1         systemfonts_1.3.1      
##  [43] bbmle_1.0.25.1          tools_4.4.1             progress_1.2.3         
##  [46] treeio_1.28.0           snow_0.4-4              Rcpp_1.1.0             
##  [49] glue_1.7.0              gridExtra_2.3           SparseArray_1.4.8      
##  [52] xfun_0.46               qvalue_2.36.0           withr_3.0.2            
##  [55] numDeriv_2016.8-1.1     fastmap_1.2.0           digest_0.6.36          
##  [58] timechange_0.3.0        R6_2.5.1                gridGraphics_0.5-1     
##  [61] colorspace_2.1-1        GO.db_3.19.1            dichromat_2.0-0.1      
##  [64] RSQLite_2.3.7           R.methodsS3_1.8.2       generics_0.1.3         
##  [67] data.table_1.16.0       prettyunits_1.2.0       graphlayouts_1.2.0     
##  [70] httr_1.4.7              htmlwidgets_1.6.4       S4Arrays_1.4.1         
##  [73] scatterpie_0.2.4        pkgconfig_2.0.3         gtable_0.3.6           
##  [76] blob_1.2.4              S7_0.2.0                XVector_0.44.0         
##  [79] shadowtext_0.1.4        htmltools_0.5.8.1       fgsea_1.30.0           
##  [82] png_0.1-8               ggfun_0.1.6             rstudioapi_0.16.0      
##  [85] tzdb_0.4.0              reshape2_1.4.4          coda_0.19-4.1          
##  [88] nlme_3.1-166            curl_5.2.3              bdsmatrix_1.3-7        
##  [91] cachem_1.1.0            parallel_4.4.1          pillar_1.10.1          
##  [94] grid_4.4.1              vctrs_0.6.5             dbplyr_2.5.0           
##  [97] evaluate_1.0.0          mvtnorm_1.3-1           cli_3.6.3              
## [100] locfit_1.5-9.10         compiler_4.4.1          rlang_1.1.4            
## [103] crayon_1.5.3            labeling_0.4.3          emdbook_1.3.13         
## [106] plyr_1.8.9              fs_1.6.4                stringi_1.8.4          
## [109] viridisLite_0.4.2       BiocParallel_1.38.0     babelgene_22.9         
## [112] Biostrings_2.72.1       lazyeval_0.2.2          GOSemSim_2.30.2        
## [115] Matrix_1.7-0            hms_1.1.3               patchwork_1.3.1        
## [118] bit64_4.5.2             KEGGREST_1.44.1         highr_0.11             
## [121] igraph_2.0.3            memoise_2.0.1           bslib_0.8.0            
## [124] ggtree_3.12.0           fastmatch_1.1-4         bit_4.5.0              
## [127] ape_5.8-1               gson_0.1.0

Differential Gene Expression Analysis of EWSR1-FLI1 Knockdown in Ewing Sarcoma

Ishani Ray

2025-10-19

Introduction

Experimental Design

Setup

Exploratory Data Analysis

Sample Overview

Principal Component Analysis

Differential Gene Expression Analysis

DESeq2 Analysis

Summary

MA Plot

Differentially Expressed Genes

Volcano Plot

Expression Heatmap

Pathway Enrichment Analysis

Up-regulated Pathways KEGG

Down-regulated Pathways KEGG

Hallmark Pathways (GSEA)

ORA: Up-regulated Pathways

ORA: Down-regulated Pathways

Summary

Key Findings

Results Summary Table

Conclusions

Session Information