se.rmd

---
title: "Quantitative Transcriptional Biomarkers of Xenobiotic Receptor Activation in Rat Liver"
author: "Lenny Kovac 3765442"
date: "`r Sys.Date()`"
output:
  html_document:
    theme: united
    fig_caption: true
    toc: yes
    toc_float:
      collapsed: true
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(
  echo       = TRUE,
  warning = FALSE,
  message = FALSE,
  fig.align  = "center",
  fig.width  = 7,
  fig.height = 5
)

# TODO: Move installs into Dockerfile
if (!requireNamespace("rtracklayer", quietly = TRUE))
  BiocManager::install("rtracklayer")
if (!requireNamespace("org.Rn.eg.db", quietly = TRUE))
  BiocManager::install("org.Rn.eg.db")
if (!requireNamespace("janitor", quietly = TRUE))
  BiocManager::install("janitor")

suppressPackageStartupMessages({
    library(SummarizedExperiment)
    library(DESeq2)
    library(vsn)
    library(ggplot2)
    library(ggrepel)
    library(ComplexHeatmap)
    library(RColorBrewer)
    library(hexbin)
    library(iSEE)
    library(tidyverse)
    library(rtracklayer)
    library(janitor)
    library(ExploreModelMatrix)
    library(cowplot)
    library(apeglm)
    library(org.Rn.eg.db)
    library(dplyr)
    library(clusterProfiler)
})

# =============================================================================
# 0. Load Feature count files from galxy each sample has its own file
# =============================================================================
data_dir <- "../data/"

# we have to mask certain functions cuz the name exists in other packages
select <- dplyr::select
filter <- dplyr::filter
rename <- dplyr::rename

# Find all the count files from our galxy pipeline
files <- list.files(
  path       = data_dir,
  pattern    = "\\.tabular",
  full.names = TRUE,
  recursive = TRUE
  )

# the experiment uses different ids for the sample. we want to annotate our samples later so we have to map to these
id_map <- tibble(
  gsm_id = c("GSM4283548","GSM4283552","GSM4283549","GSM4283553",
             "GSM4283550","GSM4283554","GSM4283551","GSM4283555"),
  us_id  = c("US-1420908","US-1420919","US-1420917","US-1420905",
             "US-1420872","US-1420907","US-1420937","US-1420893")
)

# Named vector: names = GSM, values = US
gsm_to_us <- setNames(id_map$us_id, id_map$gsm_id)

read_count_file <- function(path) {
  df <- read_tsv(path, show_col_types = FALSE)

  # The 2nd column name IS the sample id ("GSM4283555")
  sample_id <- colnames(df)[2]

  df %>%
    rename(count = 2) %>%
    mutate(
      sample_id = sample_id,
      us_id     = unname(gsm_to_us[sample_id]),   # lookup US-ID
    )
}

# Read all files and stack them into one long data frame
feature_counts_tbl <- files %>%
  map(read_count_file) %>%
  list_rbind()%>%
  rename(gene_id=Geneid)

# =============================================================================
# 1. Load sample annotation file, we have to filter for our samples
# =============================================================================
meta <- read_tsv("../data/GSE144219_Rx-TGx.Meta.txt", show_col_types = FALSE)
annotated_samples_tbl <- meta %>%
  filter(`Sample name` %in% id_map$us_id)

# =============================================================================
# 2. Load gene annotation file from ucsc rn7 and filter for genes we have
# =============================================================================
gtf <- import("../data/mRatBN7.gtf")
gtf_tbl <- as_tibble(gtf)

# Extract gene-level annotation, filtered to my feature_cpount genes
annotated_genes_tbl <- gtf_tbl %>%
  select(
    gene_id,
    type,
    chr       = seqnames,
    start, end, strand,
    source
  ) %>%
  distinct(gene_id, .keep_all = TRUE) %>%
  # Keep ONLY genes in my feature counts
  semi_join(feature_counts_tbl, by="gene_id")
```

# Experiment / Data-Aquisition

The raw RNA-seq data originates from DOI: 10.1093/toxsci/kfaa026 and was pre-processed using a Galaxy pipeline prior to analysis in R. The experiment consists of 8 rat liver samples divided into two groups: Methimazole-treated and Control (n=4).

The simplified upstream Galaxy workflow:

1. Quality Control (FastQC) — Raw reads were assessed for sequencing quality, adapter contamination, and GC-content biases to determine whether trimming was necessary.

2. Trimming (Trimmomatic) — Reads were trimmed based on the quality control report to remove low-quality bases and adapter sequences, improving downstream mapping accuracy.

3. Mapping (HISAT2) — Trimmed reads were aligned against the rat reference genome rn7 (mRatBN7.2), producing per-sample alignment files.

4. Feature Counting (featureCounts) — Mapped reads were summarized at the gene level by counting how many reads overlapped annotated genes in the rn7 GTF, producing the raw count matrices used as input for this analysis.


```{r Build Experiment}

# =============================================================================
# 0. Construct The Experiment Object
# =============================================================================
# We have 3 Tables: annotated_genes_tbl, feature_counts_tbl and annotated_samples_tbl

# ---- 0.1 Build the count matrix from long feature_counts
count_matrix <- feature_counts_tbl %>%
  select(gene_id, sample_id, count) %>%
  pivot_wider(names_from = sample_id, values_from = count) %>%
  column_to_rownames("gene_id") %>%
  as.matrix()

# ---- 0.2 Reorder GENE annotations to match matrix ROWS 
annotated_genes_ordered_tbl<- annotated_genes_tbl %>%
  filter(gene_id %in% rownames(count_matrix)) %>%
  # arrange exactly as the matrix rows 
  arrange(match(gene_id, rownames(count_matrix))) %>%
  column_to_rownames("gene_id")

# ---- 0.3 Reorder SAMPLE annotations to match matrix COLUMNS
annotated_samples_ordered_tbl <- annotated_samples_tbl %>%
  # join GSM onto the US-keyed annotations
  left_join(id_map, by = c("Sample name" = "us_id")) %>%
  # keep only samples present in the matrix
  filter(gsm_id %in% colnames(count_matrix)) %>%
  # order exactly as the matrix columns (which are GSM ids)
  arrange(match(gsm_id, colnames(count_matrix))) %>%
  # use GSM as rownames so it matches the matrix
  column_to_rownames("gsm_id")

# =============================================================================
# 1. Verify the correct Alignment, stop if its wrong!!
# =============================================================================
stopifnot(
  identical(rownames(count_matrix), rownames(annotated_genes_ordered_tbl)),
  identical(colnames(count_matrix), rownames(annotated_samples_ordered_tbl))
)

# =============================================================================
# 2. Build the acutal experiment
# =============================================================================
se <- SummarizedExperiment(
  assays  = list(counts = count_matrix),
  rowData = annotated_genes_ordered_tbl,
  colData = annotated_samples_ordered_tbl
)
```
The three annotation tables [gene metadata, sample metadata, and the count matrix] are combined into a single SummarizedExperiment object. This ensures that rows (genes) and columns (samples) remain consistently aligned throughout the analysis. The colData slot holds sample-level information such as treatment group, while rowData holds gene-level genomic coordinates.

## Some Stats pre R processing.

```{r Experiment Stats0}
nrow(se)
```
Total number of genes. But we want to avaoid the cellular "background noise". We can do that by filtering for genes with a higher count. In our case we will from now only observe genes with a count > 5.

```{r Experiment Stats1}
# =============================================================================
# 0. Extract genes with count > 5 and build lib_size for each sample.
# =============================================================================
se <- se[rowSums(assay(se, "counts")) > 5, ]
nrow(se)
```

# Explorative Daten Analyse (EDA)

Before testing for differential expression, we explore the structure of the
count data to assess data quality, identify technical biases, and confirm that
Methimazole vs. Control is the dominant source of variation.

## Library Size Variation

The library size (total number of mapped reads per sample) reflects sequencing
depth. Large differences in library size between samples are a technical
artifact rather than a biological signal, and must be accounted for through
normalization prior to differential expression analysis.

```{r EDA0, fig.cap="Libary Size for each sample."}

# =============================================================================
# 1. Add an adjust for lib_size.
# =============================================================================

#lib-size represents how many genes we have for a given sample
se$libSize <-  colSums(assay(se))
colnames(colData(se))[colnames(colData(se)) == "characteristics: treatment"] <- "treatment"

#simple bar plot for our lib sizes per sample looks good ;) only one outlier
colData(se) %>%
  as.data.frame() %>%
  ggplot(aes(x = `Sample.name`, y = libSize / 1e6, fill = treatment)) +
         geom_col(stat = "identity") + theme_bw() +
         labs(title = "Library size per sample",
              x = "Sample", y = "Total count in millions") +
         theme(axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1))

# init DESeq2 dataset object. ~ treatment tells model to test for differences in gene expression driven by the treatment variable ( methimazole vs. control)
dds <- DESeq2::DESeqDataSet(se, design = ~ treatment)

# compute sizeFactors for libary size (we only have one outlier but still account for it)
dds <- estimateSizeFactors(dds)

# sizefactors against absolute libarysize
ggplot(data.frame(libSize = colSums(assay(dds)),
                  sizeFactor = sizeFactors(dds),
                  treatment = dds$treatment),
       aes(x = libSize, y = sizeFactor, col = treatment)) +
    geom_point(size = 5) + theme_bw() +
    labs(title = "Size factors vs library size",
      x = "Library size", y = "Size factor")
```

```{r EDA1, fig.cap="The raw count data shows a mean-variance dependency"}
# standard abweichung against sizefacotrs (as aspected our one outlier)
meanSdPlot(assay(dds), ranks = FALSE)

```

```{r EDA2, fig.cap="Variance stabilizing transformation (VST) is applied to remove this dependency before distance-based analyses"}
# rename cuz deseq2 has some reserved names :(
rd <- rowData(dds)
colnames(rd)[colnames(rd) == "start"] <- "gene_start"
colnames(rd)[colnames(rd) == "end"] <- "gene_end"
colnames(rd)[colnames(rd) == "strand"] <- "strandedness"
rowData(dds) <- rd

vsd <- DESeq2::vst(dds, blind = TRUE)
meanSdPlot(assay(vsd), ranks = FALSE)
```

```{r EDA3, fig.cap="A clear block structure separating methimazole-treated from control samples indicates that treatment is the primary driver of transcriptional variation."}
# cluster samples with a heatmap
dst <- dist(t(assay(vsd)))
colors <- colorRampPalette(brewer.pal(9, "Blues"))(255)
ComplexHeatmap::Heatmap(
    as.matrix(dst),
    col = colors,
    name = "Euclidean\ndistance",
    cluster_rows = hclust(dst),
    cluster_columns = hclust(dst),
    bottom_annotation = columnAnnotation(
        Treatment = vsd$treatment,
        col = list(
            Treatment = c("methimazole" = "red", "-control-" = "blue")
        )
    )
)

```

```{r EDA4, fig.cap= "PC1 captures 64% of all transcriptional variation and perfectly separates the two groups.Methimazole treatment is by far the dominant source of variation in the dataset"}
# PCA
pcaData <- DESeq2::plotPCA(vsd, intgroup = c("treatment"),
                           returnData = TRUE)

percentVar <- round(100 * attr(pcaData, "percentVar"))

ggplot(pcaData, aes(x = PC1, y = PC2)) +
    geom_point(aes(color = treatment), size = 5) +
    theme_minimal() +
    xlab(paste0("PC1: ", percentVar[1], "% variance")) +
    ylab(paste0("PC2: ", percentVar[2], "% variance")) +
    coord_fixed() +
    scale_color_manual(values = c("-control-" = "blue", "methimazole" = "red")) +
    labs(title = "PCA of variance-stabilized counts")
```

# Differantial Gene Expression

Having confirmed sample quality and treatment-driven variation, we proceed to formal statistical testing for
differential expression between methimazole-treated and control samples using DESeq2.Shrinkage of log2 fold-change
estimates is applied to stabilize estimates for low-count genes.

```{r DGE0, fig.cap= "Dispersions decreases as mean expression increases, following the fitted trend line. Gene having similar counts across all samples or with low mean counts are concentrated around the low dispersion."}
# Dispersions
dds <- estimateDispersions(dds)
plotDispEsts(dds)
```

```{r DGE1 }
# For each gene, DESeq2 has fit a negative binomial model.
# Is the log2 fold-change between methimazole and control significantly different from zero?
dds <- nbinomWaldTest(dds)

# Results
res_groups <- results(dds, name = "treatment_methimazole_vs_.control.")
summary(res_groups)
```
20,525 genes that survived erlier filter (rowSums(counts) > 5). With adjusted p-value < 0.1

LFC > 0 (up) : 1266, 6.2% genes are significantly upregulated by methimazole (positive log2 fold-change).

LFC < 0 (down) : 1212, 5.9% genes are significantly downregulated by methimazole.

17 genes were flagged as outliers which is good. 35% too low counts were removed.

```{r DGE2, fig.cap= "Shows expression of top DE genes across all 8 samples:"}
#A gene with counts 2 vs 6 shows a huge fold-change but is unreliable.
#so  top hits are genes with both large AND trustworthy changes.
res_groupsLfc <- lfcShrink(dds, coef = "treatment_methimazole_vs_.control.", res = res_groups)

# Transform counts
vsd <- vst(dds, blind = TRUE)

# Get top DE genes
genes <- res_groupsLfc[order(res_groupsLfc$pvalue), ] %>%
    head(10) %>%
    rownames()

heatmapData <- assay(vsd)[genes, ]

# annotaion
heatmapColAnnot <- data.frame(colData(vsd)[, c("treatment")])
heatmapColAnnot <- HeatmapAnnotation(df = heatmapColAnnot)

# Plot as heatmap
ComplexHeatmap::Heatmap(heatmapData,
                        top_annotation = heatmapColAnnot,
                        cluster_rows = TRUE, cluster_columns = FALSE)
```

# Gene Enrichment Analysis

## Gene Ontology and KEGG Enrichment

I decided for two enrichment options. KEGG and Gene Ontology. GO is more a functional annotations whereas KEGG is a
database for Pathways. KEGG is way smaller thats why I have implemented both.

```{r Gene Enrichtent Analyis}
# ========================================================
# Significant genes (RefSeq IDs from DESeq results)
# Threshold: padj < 0.05 & |log2FC| > log2(1.5)
# ========================================================

res_groups_df   <- as.data.frame(subset(res_groups,
                     padj < 0.05 & abs(log2FoldChange) > log2(1.5)))

sig_refseq      <- gsub("\\..*$", "", rownames(res_groups_df))           # strip version .1/.2
universe_refseq <- gsub("\\..*$", "", rownames(as.data.frame(res_groups)))

# ---- 1. Map RefSeq -> Entrez (enrichGO/KEGG need Entrez IDs)
sig_entrez <- AnnotationDbi::mapIds(org.Rn.eg.db,
                  keys = sig_refseq, keytype = "REFSEQ",
                  column = "ENTREZID", multiVals = "first")

universe_entrez <- AnnotationDbi::mapIds(org.Rn.eg.db,
                  keys = universe_refseq, keytype = "REFSEQ",
                  column = "ENTREZID", multiVals = "first")

# ---- 2. Clean 
sig_entrez      <- unique(na.omit(sig_entrez))
universe_entrez <- unique(na.omit(universe_entrez))

# ---- 3. GO enrichment
resGroupGO <- enrichGO(
  gene          = sig_entrez,
  universe      = universe_entrez,
  OrgDb         = org.Rn.eg.db,
  keyType       = "ENTREZID",
  ont           = "BP",
  pAdjustMethod = "BH",
  pvalueCutoff  = 0.05,
  qvalueCutoff  = 0.05,
  readable      = TRUE
)

# ---- 4. KEGG enrichment (for rat, KEGG IDs == Entrez IDs)
# queries KEGG REST API -> requires internet connection to knit
resGroupKEGG <- enrichKEGG(
  gene          = sig_entrez,
  universe      = universe_entrez,
  organism      = "rno",
  keyType       = "kegg",
  pAdjustMethod = "BH",
  pvalueCutoff  = 0.05,
  qvalueCutoff  = 0.05
)
```

## Both enrichment tools compared.

```{r goo, fig.cap="Go Enrichment."}
dotplot(resGroupGO,   showCategory = 15) + ggtitle("GO Biological Process Enrichment")
```

```{r keggg, fig.cap="KEGG Enrichment"}
dotplot(resGroupKEGG, showCategory = 15) + ggtitle("KEGG Pathway Enrichment")

```

```{r bar, fig.cap="Best UP and DOWN  regulated pathways."}
# ========================================================
# 0. Helper: map RefSeq -> Entrez for a subset
# ========================================================
map_to_entrez <- function(refseq_ids) {
  na.omit(AnnotationDbi::mapIds(org.Rn.eg.db,
          keys = gsub("\\..*$", "", refseq_ids),
          keytype = "REFSEQ", column = "ENTREZID", multiVals = "first"))
}

# Helper: compute log2 enrichment factor for an enrichResult table
add_log2_enrichment <- function(enrichObj) {
  tbl  <- as.data.frame(enrichObj)
  n_11 <- tbl$Count
  n_10 <- length(intersect(enrichObj@gene, enrichObj@universe))
  n_01 <- as.numeric(gsub("/.*$", "", tbl$BgRatio))
  n    <- length(enrichObj@universe)
  tbl$log2_Enrichment <- log2((n_11 / n_10) / (n_01 / n))
  tbl
}

# ========================================================
# 1. obtain UP-regulated genes
# ========================================================
groupsDEup <- as.data.frame(subset(res_groups,
                 padj < 0.05 & log2FoldChange >  log2(1.5)))
up_entrez  <- map_to_entrez(rownames(groupsDEup))

resGroupGOup <- enrichGO(gene         = up_entrez,
                         keyType      = "ENTREZID",
                         ont          = "BP",
                         OrgDb        = org.Rn.eg.db,
                         universe     = universe_entrez,
                         pvalueCutoff = 1, qvalueCutoff = 1,   # keep all for plotting
                         readable     = TRUE)

resGroupGOupTable <- add_log2_enrichment(resGroupGOup)

# ========================================================
# 2. obtain DOWN-regulated genes
# ========================================================
groupsDEdown <- as.data.frame(subset(res_groups,
                   padj < 0.05 & log2FoldChange < -log2(1.5)))
down_entrez  <- map_to_entrez(rownames(groupsDEdown))

resGroupGOdown <- enrichGO(gene         = down_entrez,
                           keyType      = "ENTREZID",
                           ont          = "BP",
                           OrgDb        = org.Rn.eg.db,
                           universe     = universe_entrez,
                           pvalueCutoff = 1, qvalueCutoff = 1,
                           readable     = TRUE)

resGroupGOdownTable <- add_log2_enrichment(resGroupGOdown)

# 2.3. Wrap long descriptions
resGroupGOupTable$Description   <- str_wrap(resGroupGOupTable$Description,   width = 40)
resGroupGOdownTable$Description <- str_wrap(resGroupGOdownTable$Description, width = 40)

# ========================================================
# 3. Combined barplot: top 5 up + top 5 down
# ========================================================
upTop   <- head(resGroupGOupTable,   5)
downTop <- head(resGroupGOdownTable, 5)

plotDF <- rbind(upTop, downTop)
plotDF$direction <- c(rep("up", nrow(upTop)), rep("down", nrow(downTop)))

ggplot(plotDF,
       aes(x = log2_Enrichment,
           y = factor(Description, levels = rev(Description)),
           fill = direction)) +
    geom_bar(stat = "identity") +
    scale_fill_manual(values = c("up" = "red", "down" = "darkgreen")) +
    geom_text(aes(label = sprintf("%.2e", p.adjust)),
              hjust = 1, col = "white", size = 3) +
    ylab("") + xlab("log2 Enrichment") +
    ggtitle("Top GO:BP Terms — Up vs Down Regulated") +
    theme_bw()

```

```{r volcano, fig.cap="A nice volcano plot with annotated genes. Higher log2fold becuase we are NOT looking at pathways"}
# ========================================================
# 4. VOLCANO PLOT
# ========================================================
volDF <- as.data.frame(res_groups)
volDF$RefSeq <- rownames(volDF)
volDF <- volDF[!is.na(volDF$padj) & !is.na(volDF$log2FoldChange), ]

lfc_cut  <- log2(1.5)
padj_cut <- 0.05

volDF$Significance <- "Not Significant"
volDF$Significance[volDF$padj < padj_cut & volDF$log2FoldChange >  lfc_cut] <- "Up"
volDF$Significance[volDF$padj < padj_cut & volDF$log2FoldChange < -lfc_cut] <- "Down"
volDF$Significance <- factor(volDF$Significance,
                             levels = c("Not Significant", "Up", "Down"))

# ---- Pick top genes to label (most significant up + down)
upSet   <- volDF[volDF$Significance == "Up",   ]
downSet <- volDF[volDF$Significance == "Down", ]
topUp   <- head(upSet[order(upSet$padj),   ], 5)
topDown <- head(downSet[order(downSet$padj), ], 5)
labelDF <- rbind(topUp, topDown)

#  map RefSeq -> gene SYMBOL for readable labels
labelDF$Symbol <- AnnotationDbi::mapIds(org.Rn.eg.db,
                    keys = gsub("\\..*$", "", labelDF$RefSeq),
                    keytype = "REFSEQ", column = "SYMBOL", multiVals = "first")
# fall back to RefSeq if no symbol found
labelDF$Symbol[is.na(labelDF$Symbol)] <- labelDF$RefSeq[is.na(labelDF$Symbol)]

# nice plot 
ggplot(volDF, aes(x = log2FoldChange, y = -log10(padj), color = Significance)) +
    geom_point(alpha = 0.7, size = 1.8) +
    scale_color_manual(values = c("Not Significant" = "#cde7e3",
                                  "Up"              = "#e8666e",
                                  "Down"            = "#4f7faf")) +
    geom_vline(xintercept = c(-lfc_cut, lfc_cut), lty = 2, col = "grey50") +
    geom_hline(yintercept = -log10(padj_cut),     lty = 2, col = "grey50") +
    geom_text_repel(data = labelDF, aes(label = Symbol),
                    size = 3, color = "black", max.overlaps = Inf,
                    box.padding = 0.4, segment.color = "grey60") +
    labs(title = "Volcano Plot - Differential Expression",
         x = "log2 Fold Change", y = "-log10(Adjusted P-value)") +
    theme_bw() +
    theme(plot.title   = element_text(hjust = 0.5, face = "bold"),
          legend.title = element_text())
```