Differential Gene Expression Analysis
in Atopic Dermatitis
Using DESeq2

Introduction to Differential Gene Expression in Atopic Dermatitis (AD) Patients

• In this analysis, we perform Differential Gene Expression (DEG) analysis using the DESeq2 package on the GSE157194 dataset, which consists of RNA-seq data from patients with Atopic Dermatitis (AD). This dataset includes a total of 166 samples obtained from 57 patients with moderate to severe AD, prior to and after systemic treatment.

Dataset Overview:

• GSE157194: RNA-seq data from AD patients.
• 111 samples before treatment: Includes lesional (AL) and non-lesional (AN) skin biopsies. Non-lesional samples were taken at least 5 cm from the active lesion, ensuring that the comparison between these conditions captures the underlying gene expression differences in the disease state before any intervention.
• 40 samples post-treatment with Dupilumab (3 months).
• 15 samples post-treatment with Cyclosporine (3 months).

Objective of This Analysis:

• Our primary focus is to conduct a differential expression analysis between the lesional (AL) and non-lesional (AN) skin samples before treatment. This helps in identifying key differentially expressed genes (DEGs) and associated pathways that distinguish diseased skin from non-diseased skin. Insights gained from this comparison can:
• Reveal potential biomarkers for AD diagnosis and monitoring.
• Highlight novel therapeutic targets or pathways relevant to AD pathogenesis.
• Provide groundwork for drug discovery and personalized medicine by identifying pathways that could respond to specific therapeutic interventions.

Comprehensive Analysis Approach:

• This analysis incorporates various bioinformatics techniques and visualization tools (e.g., PCA, volcano plots, etc.) to ensure a holistic understanding of the gene expression profiles. While some aspects of the analysis may provide standalone insights, the goal is to create a comprehensive reference for DEG analysis using DESeq2. This can serve as a functional analysis resource for future studies.

In addition to the DEG analysis, we also perform pathway and gene enrichment analysis, including:

• KEGG and GO term analyses.
• Gene Set Enrichment Analysis (GSEA) and Overrepresentation Analysis (ORA) to explore the biological processes, molecular functions, and cellular components involved in AD.

1. Load Necessary Libraries

Before starting the analysis, we need to load the essential R libraries required for data manipulation, visualization, and performing differential gene expression (DEG) analysis.

2. Load and Prepare Data

2.1. Retrieving Data from GEO

• In this step, we download the raw gene expression data (count matrix) and the sample information (referred to as metadata or phenotype data) from the Gene Expression Omnibus (GEO) using GEO's accession ID GSE157194. The count matrix contains raw gene expression counts, while the phenotype data contains metadata for each sample (such as treatment groups and sample type).

• In this code (Code Block 2.1), we begin by downloading and preparing the dataset for analysis. Step 1 involves defining the GEO accession ID (gse_id <- "GSE157194") and using getGEOSuppFiles() from the GEOquery library to download the supplementary count matrix, which is then loaded into R using read.table() after decompressing the .gz file. In Step 2, the phenotype data (metadata) is obtained using getGEO() from GEOquery, and stored in metadata. The metadata provides additional descriptive information about the samples in the count matrix. Step 3 prints summary information ( Figure 2.1), including the number of genes and samples in the count matrix and the number of samples and features in the metadata, after downloading the data from GEO, for a quick overview of the dataset.

2.2. Preparing Data

• Effective data preparation is crucial for ensuring high-quality and accurate results in downstream analyses. This step includes identifying and removing potential outlier samples to avoid biases and pre-filtering low-count genes to enhance the overall data quality.

2.2.1. Identify and Remove Outlier Samples

• Outlier samples, particularly those with significantly low or high library sizes, can introduce biases into the analysis and skew the results. Therefore, it is crucial to identify and remove such samples before conducting further data manipulations to ensure that the remaining data is representative of high-quality samples, allowing for reliable downstream analyses.
• Traditionally, removing samples below the lower threshold (calculated as Q1 - 1.5 * IQR) is a common practice, as these samples often indicate technical issues, such as failed sequencing or low quality. Samples above the upper threshold (Q3 + 1.5 * IQR) are typically retained, as they are often biologically meaningful (e.g., having high sequencing depth), and removing them could risk discarding valuable information. However, if it can be demonstrated that such samples introduce bias or are clearly problematic, it might be necessary to remove them to maintain data integrity.
• Based on the dot plot of library sizes, we identified a total of six outlier samples from 166 original samples: three with significantly low counts (Patient_12_AN_m3, Patient_24_AN_m3, Patient_31_AL_m0) and three with significantly high counts (Patient_41_AL_m3, Patient_42_AL_m0, Patient_48_AL_m3). Although the high-count samples might be biologically distinct, their presence could mask potential biomarkers within the majority of the population by over-influencing the analysis . Therefore, to ensure that we focus on identifying biomarkers that are representative of the majority, we have decided to remove both lower and upper threshold outliers.
• The provided code is written to handle either scenario, allowing users to remove only the lower threshold outliers, or both lower and upper threshold outliers. Comments within the code explain which lines to uncomment or comment out based on the user's chosen approach. In this case, both thresholds are applied to remove the six outlier samples, resulting in a filtered dataset that better represents the core population.
• Not all the plots presented in the code are necessary to generate; they are provided here to demonstrate different ways of visualizing the effect of removing outlier samples.
• In this code (Code Block 2.2.1), we begin with Step 1, which involves calculating the library sizes for each sample using colSums(counts_matrix). This operation returns the total count of all reads for each sample. The summary statistics (summary(library_sizes)) is then printed to provide an overview of the data distribution (Min= 9,628,756, 1st Qu=20,647,659, Median=22,896,256, Mean=23,218,539, 3rd Qu= 25,719,325, Max= 43,076,805). Next, in Step 1.1, we set two thresholds for outlier identification, using quantiles (Q) and the interquartile range (IQR). The first quartile (quantile(library_sizes, 0.25)) and the third quartile (quantile(library_sizes, 0.75)) represent the 25th and 75th percentiles of the data, respectively. The IQR (IQR(library_sizes)) indicates the spread of the central 50% of the data and helps to assess data variability. We then calculate the lower threshold as Q1 - 1.5 * IQR, which targets low outliers that may indicate poor-quality or problematic samples. Optionally, the upper threshold is calculated as Q3 + 1.5 * IQR to identify samples with unusually high counts, which may also be removed if deemed necessary. In Step 1.2, depending on the user's choice, either or both thresholds can be used to define outlier samples. These outliers are identified using logical conditions, and their indices are used to remove them from both "counts_matrix" and "metadata", resulting in the creation of filtered versions ("counts_matrix_filtered" and "metadata_filtered") (Step 3). After identifying the outliers, in Step 2, we store the original dimensions of both "counts_matrix" and "metadata" before filtering. This allows us to compare them against the dimensions after outlier removal, which we store in Step 3.3. In Step 3, we remove the outlier samples from the "counts_matrix" (Step 3.1) and "metadata" (Step 3.2), ensuring that the metadata matches the "filtered counts_matrix". After filtering, we reassign the filtered versions back to "counts_matrix" and "metadata" to ensure consistency for subsequent analysis steps (Step 6). In Step 4, following the removal of outliers, we create various visualizations (dot plot (Figure 2.2.1.1), bar plot (Figure 2.2.1.2), boxplot (Figure 2.2.1.3), density plot (Figure 2.2.1.4), and violin plot (Figure 2.2.1.5)) to illustrate the changes in library sizes before and after filtering. These visualizations provide insights into the effects of outlier removal on the dataset. Specifically:

• Boxplots and Violin plots are used to compare the distributions of library sizes before and after filtering.
• Dot plots annotate specific outlier samples that were removed.
• Bar plots allow for a direct comparison of individual sample library sizes before and after filtering.
• Density plots provide an overview of the entire distribution of library sizes.

• In Step 5, we generate visualizations to confirm the changes in the dimensions of the "counts_matrix" and "metadata" before and after outlier removal using dimensions comparison plot (Figure 2.2.1.6) as well as Summary Information (Figure 2.2.1.7). To achieve this, we create data frames representing the count matrix dimensions and metadata dimensions, both before and after filtering (Step 5.1 and Step 5.2). These data frames are then used to generate four bar plots, which illustrate the number of genes, samples, and metadata features before and after filtering. These four plots are combined into a single figure using the grid.arrange() function from the gridExtra library (Step 5.3). This combined visualization provides a clear overview of how filtering affected the dataset dimensions.

2.2.2. Pre-filtering Low-Count Genes to Improve Data Quality

• Low-count genes, which are expressed at insufficient levels, can contribute noise and increase the likelihood of false-positive results. Removing these lowly-expressed genes ensures that the analysis focuses on biologically meaningful signals and helps to reduce computational complexity and improve the robustness of downstream analyses.
• Although DESeq2 performs independent filtering internally during the differential expression step, performing pre-filtering at this early stage is widely accepted in analyses to improve computational efficiency, reduce memory usage, and remove genes that are unlikely to contribute meaningful biological insights.
• In this code (Code Block 2.2.2), we first (Step 1) create a logical matrix to identify genes with sufficient expression (counts ≥ 10). Specifically, we retain only those genes that have counts of at least 10 in at least half of the samples (ncol(counts_matrix) / 2). This step ensures that the filtered dataset retains genes with meaningful levels of expression across a significant number of samples. After identifying the genes to keep, we filter the counts_matrix to retain only those genes, reassign the filtered matrix back to counts_matrix (Step 2), and then check the dimensions of the filtered dataset (Step 3) to confirm the number of genes and samples that remain (Figure 2.2.2).

3. Exploratory Data Analysis (EDA)- Before the "DESeqDataSet" object (dds) Creation

• EDA is essential before performing DEA (Differential Expression Analysis) to explore the relationship between samples and to gain insights into the structure and distribution of the data. This step involves exploring both the metadata and the count matrix to ensure that the dataset is well-prepared for downstream analysis.

3.1. Metadata Exploration

3.1.1. Metadata Overview

• We examine the metadata to identify key columns for analysis, such as sample ID, experimental conditions, tissue type, time points, and treatments. Selected columns are then viewed using View() function for a quick overview of the samples (Code Block 3.1.1.; Figure 3.1.1).

3.1.2. Check Distribution of Samples Across Key Metadata Variables (using table() and kable())

3.1.2.1. Summary of Sample Distribution by Patient Group and Timepoints

• By exploring key metadata variables such as sample type (e.g., lesional vs non-lesional or treated vs untreated), timepoints (e.g., baseline vs post-treatment), and therapy groups (e.g., Cyclosporine and Dupilumab), we gain insights into the distribution of the samples across different conditions. This exploration also helps to ensure that the dataset is well-balanced and appropriately structured for the upcoming differential expression analysis.
• For example, after filtering out sample outliers, the source_name_ch1 column in metadata shows the distribution of samples across various experimental conditions (Code Block 3.1.2; Figure 3.1.2.1 a), such as baseline samples (m0) for both lesional (AL) (m0_AL: 55 samples) and non-lesional (AN) (m0_AN: 54 samples) skin, which are key to our current analysis, as well as samples taken post-treatment with Cyclosporine (8 lesional & 6 non-lesional) and Dupilumab (19 lesional & 18 non-lesional samples). This data can be useful for exploring treatment effects in future experiments.
• In addition, characteristics_ch1.2 (Code Block 3.1.2.1; Figure 3.1.2.1 b) provides information on the timepoints (e.g., baseline or post-treatment). This is crucial for analyzing how gene expression changes over time, particularly in response to therapies like Cyclosporine and Dupilumab. In this dataset, 51 samples were collected 3 months post-treatment, and 109 samples were collected at baseline.

3.1.2.2. Summary of Sample Distribution Across Multiple Metadata Variables

• Other metadata columns can also be explored (Code Block 3.1.2.2):
• The characteristics_ch1.3 column (Figure 3.1.2.2 a) provides information about the treatment each patient received. It shows that 14 samples were treated with Cyclosporine and 37 samples with Dupilumab, while the remaining 109 samples represent baseline conditions (m0_AL and m0_AN combined). This is important for distinguishing between treatment groups when analyzing differential expression across conditions.
• The patient id:ch1 column (Figure 3.1.2.2 b) helps track individual patients across different conditions. Some patients may have samples from multiple conditions (e.g., lesional, non-lesional, and post-treatment with both Cyclosporine and Dupilumab), allowing for longitudinal analyses within the same individuals. For researchers, this provides an opportunity to focus on patients who have undergone multiple treatments, reducing variability and enhancing the accuracy of comparative analyses.
• In addition to these, other metadata columns can provide insights into the experimental setup and data processing. For example, platform_id shows all utilized platforms (GPL21290 in this experiment), extract_protocol_ch1 and extract_protocol_ch1.1 describe how the RNA was extracted and processed in the lab, and how the RNA libraries were prepared for sequencing. Furthermore, data_processing and its related fields provide detailed information on how the sequencing data was processed, from basecalling using Illumina software to trimming primers with cutadapt and mapping reads to the human reference genome with Tophat2 and Bowtie2. The use of Samtools to clean and sort the reads, and HTSeq to count reads per gene, ensures the data is accurately prepared for downstream analysis. The genome build used for mapping is GRCh38, and the raw count data is stored in a tab-delimited format.

3.1.3. Check Distribution of Samples Across Key Metadata Variables (using bar charts)

• Visualizing the distribution of samples with bar charts provides a more intuitive and visual approach compared to table-based summaries. It helps to quickly identify patterns, imbalances, or trends in the data, making it easier to understand the distribution across different experimental conditions or metadata variables.

3.1.3.1. Visualizing Sample Distribution for a Single Metadata Variable Using Bar Charts

• This code (Code Block 3.1.3.1) focuses on visualizing a single variable, such as source_name_ch1, to display the sample distribution by experimental condition (Figure 3.1.3.1). However, the variable can be easily substituted with others if needed. (Step 1): The ggplot2 library is loaded for creating the plot. (Step 2): A bar plot is constructed using geom_bar() to create the bars and geom_text() to add count labels, where ..count.. refers to the computed count statistic. The plot is styled using theme_minimal() and color is applied with scale_fill_brewer(). (Step 3): Finally, the ggsave() function saves the plot to a specified location.

3.1.3.2. Automated Bar Chart Generation for Multiple Metadata Variables

• This approach extends the visualization to multiple metadata variables (such as "source_name_ch1" (Figure 3.1.3.1), "characteristics_ch1.2" (Figure 3.1.3.2 a), "characteristics_ch1.3" (Figure 3.1.3.2 b), "patient id:ch1" (Figure 3.1.3.2 c)), allowing for the exploration of different aspects of the dataset. It automates the creation of bar charts, making it easy to visualize sample distributions across various metadata fields. (Step 1): The ggplot2 library is loaded for plotting and viridis library for using a color palette that supports many unique colors. (Step 2): A list of columns to visualize is defined. (Step 3): The code loops through each column, converting column names into symbols using sym() and sanitizing the names with gsub() for valid file naming. (Step 4): The plots are created using ggplot(), styled with theme_minimal(), and colored using scale_fill_viridis_d(option = "C"). We used scale_fill_viridis_d(option = "C") instead of scale_fill_brewer() because scale_fill_brewer() has a limited number of distinct colors to handle all categories, whereas scale_fill_viridis_d() provides more colors and ensures better contrast, making it suitable for larger datasets. Each plot is displayed and saved to a file using ggsave().

3.2. Count Matrix: Exploration and Visualization

3.2.1. Count Matrix Overview (Inspecting the Count Matrix using View())

• We examine the count matrix, which contains raw gene expression counts for each sample. This step helps us understand the structure of the data, including whether it is normalized or not, and allows us to check the column names (samples) and row names (Ensembl Gene Identifiers) to ensure the data is formatted correctly for the next steps of differential expression analysis (Code Block 3.2.1; Figure 3.2.1).

3.2.2. Boxplot of Raw Counts

• In RNA-seq analysis, boxplots of raw counts provide a crucial view of the distribution of gene expression across samples. By plotting the log-transformed counts, we can detect outliers, assess sequencing depth variations, and verify whether the data is ready for downstream analysis. Boxplots help identify potential issues, such as batch effects or technical biases, that may affect the accuracy of differential expression analysis (DEA) and the identification of biomarkers.

3.2.2.1. Boxplot of Raw Counts (Uncolored)

• This simple boxplot Figure 3.2.2.1 (after the removal of low-count genes and outlier samples) provides an overview of the log2-transformed gene expression counts across all samples without any group-based color distinction. It allows us to quickly assess the range of counts for each sample and identify potential outliers. This uncolored version gives a broad view of data distribution and helps to highlight sample-level variations after the removal of low-count genes and outlier samples, ensuring the data is cleaned and ready for downstream analyses.

• In this code (Code Block 3.2.2.1), we create a boxplot for the log2-transformed counts of the count matrix. (Step 1): We load the ggplot2 library to generate the plot. (Step 2): The count matrix is log2-transformed and converted into a long format for visualization. (Step 3): A basic boxplot is created with minimal styling using geom_boxplot() to show the count distribution across samples. (Step 4): The plot is saved to a specific location using the ggsave() function.

3.2.2.2. Boxplot of Raw Counts Colored by Experimental Condition

• This enhanced version of the boxplot introduces color coding based on experimental conditions. By grouping the samples according to metadata (e.g., different treatments), we can more easily visualize how the distribution of counts changes across different conditions. This plot is particularly useful for identifying condition-specific biases or batch effects that could affect downstream analyses, such as DEA for biomarker discovery. Figure 3.2.2.2 a & Figure 3.2.2.2 b show the distribution of raw counts before and after removing low-count genes and outlier samples, respectively, allowing us to assess how preprocessing affects the consistency of counts across experimental conditions.

• In this code (Code Block 3.2.2.2), we create a boxplot colored by experimental condition. (Step 1): The metadata and count matrix are aligned based on sample IDs. (Step 2): The log2-transformed count matrix is prepared for plotting and merged with metadata information. (Step 3): A boxplot is created with colors representing the experimental conditions using geom_boxplot() and a color palette. (Step 4): The total sample count is added to the legend, and (Step 5): The final plot is saved to the desired directory using the ggsave() function.

Code Block 3.2.2.2. Boxplot of Raw Counts Colored by Experimental Condition

# ============================
# Load necessary libraries
# ============================
library(ggplot2)
library(reshape2)

# ============================
# Step 1: Align Sample IDs and Metadata
# ============================

# Set the row names of metadata to the sample names (stored in metadata$title)
rownames(metadata) <- metadata$title

# Extract the sample IDs from counts_matrix (columns represent samples)
sample_ids <- colnames(counts_matrix)

# Filter metadata to contain only rows that match the sample IDs in the count matrix
metadata_filtered <- metadata[metadata$title %in% sample_ids, ]

# Reorder metadata to match the order of sample IDs in the count matrix
metadata_filtered <- metadata_filtered[match(sample_ids, rownames(metadata_filtered)), ]


# ============================
# Step 2: Prepare Data for Plotting
# ============================

# Transpose the counts matrix to make samples the rows
counts_matrix_t <- t(counts_matrix)

# Convert the transposed counts matrix to a data frame and log-transform the counts
# Add 1 to avoid log(0), which is undefined
df <- as.data.frame(log2(counts_matrix_t + 1))

# Add a sample column (rownames will become a new column for plotting)
df$Sample <- rownames(df)

# Add the group information from metadata (e.g., source_name_ch1)
df$group <- metadata_filtered$source_name_ch1[match(df$Sample, rownames(metadata_filtered))]

# Ensure `group` is a factor for proper grouping in the plot
df$group <- as.factor(df$group)


# ============================
# Step 3: Reshape the Data for ggplot2 (Long Format)
# ============================

# Melt the data frame to create a long format for ggplot2
df_long <- melt(df, id.vars = c("Sample", "group"))


# ============================
# Step 4: Define Colors for Each Unique Group and Add Sample Counts to the Labels
# ============================

# Calculate the number of samples for each group
group_counts <- table(df$group)

# Calculate the total number of samples
total_samples <- sum(group_counts)

# Create the color palette for the unique groups
unique_groups <- unique(df$group)
color_palette <- c("#FF6347", "#4682B4", "#32CD32", "#FFD700", "#8A2BE2", "#00CED1")

# Ensure the color palette has enough colors for all groups
group_colors <- setNames(color_palette[1:length(unique_groups)], unique_groups)

# ============================
# Step 5: Create the Plot with Updated Group Labels
# ============================

p <- ggplot(df_long, aes(x = factor(Sample, levels = sample_ids), y = value, fill = group)) +
  geom_boxplot() +
  theme_minimal() +
  labs(title = "Log2 Raw Counts Across Samples by Source Name", x = "Samples", y = "Log2 Counts", 
       fill = paste0("group (Total n=", total_samples, ")")) +  # Add total number to the legend title
  theme(axis.text.x = element_text(angle = 30, hjust = 1, size = 4), # Rotate and adjust x-axis labels
        plot.background = element_rect(fill = "white", colour = NA),  # Set white background
        panel.background = element_rect(fill = "white", colour = NA)) + # Set white background for the panel
  scale_fill_manual(values = group_colors, labels = paste0(levels(df$group), " (n=", group_counts[levels(df$group)], ")"))  # Use the custom color vector and add counts in the legend

# ============================
# Step 6: Save the Plot
# ============================

ggsave("Figures/3.2.2.Boxplot_Raw_Counts/3.2.2.2.Colored_Boxplot_Raw_Counts.png", 
       plot = p, width = 20, height = 6)

    

• Interpretation: In Figure 3.2.2.2 a, before removing low-count genes and outlier samples, the log2-transformed raw counts display a broad range of variability. There are several extreme outliers (black dots) far above the whiskers of the boxplots, indicating that some samples have unusually high counts. This suggests the presence of technical or biological variability that could skew downstream analysis. The boxplots themselves are wider and less consistent across different groups, reflecting the noise introduced by these low-count genes and outliers.
• After filtering, in Figure 3.2.2.2 b, the data shows a much tighter distribution. While there are still outliers (black dots), these are less extreme, and they appear closer to the whiskers, representing smaller deviations from the group median. The overall range of counts has become more consistent across groups, with fewer extreme variations, allowing for a clearer comparison between the experimental conditions. This filtering process has effectively removed noise from the data, enhancing its quality and reliability for subsequent analysis by reducing the impact of low-count genes and extreme outliers.

3.2.3. Histogram of Total Gene Counts Per Sample

• The histogram of total gene counts per sample is essential for assessing the overall sequencing depth and library size across samples in RNA-seq data. By visualizing the distribution of total gene counts (log-transformed), we can identify under-sequenced or over-sequenced samples that may need further normalization or removal during preprocessing. This visualization helps to ensure that sequencing depth is consistent across samples before differential expression analysis (DEA) and biomarker discovery. Figure 3.2.3 a & Figure 3.2.3 b display the distribution of gene counts before and after removing low-count genes and outlier samples, respectively, which helps to assess the impact of filtering on sequencing depth.

• In this code (Code Block 3.2.3), we create a histogram that visualizes the distribution of log10-transformed gene counts per sample. (Step 1): We define the output file path and specify where to save the generated plot. (Step 2): The PNG device is opened to capture the plot. (Step 3): We calculate the log10-transformed gene counts and create a histogram using the hist() function. (Step 4): To enhance the visualization, the gene count values are added as labels above each bar. Finally, (Step 5): The PNG device is closed, and the plot is saved.

Code Block 3.2.3. Histogram of Total Gene Counts Per Sample

# ============================
# Step 1: Define Output File Path
# ============================

# Define the file path where you want to save the image
output_file <- "C:/Users/Desktop/RNA_Seq_Task1/3/3_Count_Matrix_Exploration/3.2.3.a.histogram_Of_total_GeneCounts_Per_Sample_Before_Removing_LowCountGenes_AND_Outlier_Samples.png"

# ============================
# Step 2: Open PNG Device
# ============================

# Open the PNG device to save the plot
png(filename = output_file, width = 1200, height = 1000)

# ============================
# Step 3: Create Histogram of Log-Transformed Gene Counts
# ============================

# Calculate log10 of total counts per gene
log_gene_counts <- log10(rowSums(counts_matrix))

# Create the histogram and get the counts for each bin
hist_data <- hist(log_gene_counts,  
                  breaks = 50,  # Specify number of bins in the histogram
                  main = "Gene Count Distribution",  # Title of the plot
                  xlab = "Log10 Counts",  # Label for the x-axis
                  col = "steelblue")  # Add color to the bars

# Update the y-axis limit to make room for labels
ylim <- c(0, max(hist_data$counts) * 1.1)
plot(hist_data, col = "steelblue", main = "Gene Count Distribution", xlab = "Log10 Counts", ylim = ylim) # cex: to change the font size of the labels on the bars

# ============================
# Step 4: Add Numbers Above Bars
# ============================

# Add the text labels above the histogram bars
text(hist_data$mids, hist_data$counts, labels = hist_data$counts, pos = 3, cex = 1, srt = 45)

# ============================
# Step 5: Close the PNG Device
# ============================

# Close the PNG device
dev.off()
	

• Interpretation: The histograms visualize the distribution of gene counts before and after filtering low-count genes and outlier samples. Before filtering (Figure 3.2.3 a), the data exhibits a large proportion of genes with low expression levels, as shown by the sharp peak near the lower end of the Log10 scale. This suggests a significant amount of noise in the dataset, which can obscure meaningful biological signals. After filtering (Figure 3.2.3 b), the distribution shifts, with the removal of low-count genes leading to a more balanced representation of moderately to highly expressed genes. The refined dataset, with fewer extreme values and noise, is now more suitable for downstream analyses, offering a clearer and more reliable view of gene expression patterns. This cleaning process enhances data quality, improving the robustness and accuracy of subsequent biological interpretations.

# =========== Step 1: Load Required Libraries =============
library(matrixStats)  # Or DelayedMatrixStats for variance calculations
library(ggplot2)      # For enhanced plotting capabilities

# =========== Step 2: Log2 Transformation =============
# Convert counts_matrix to a matrix and apply log2 transformation
# Log transformation reduces the dynamic range of the data, making it easier to interpret
counts_matrix_log <- log2(as.matrix(counts_matrix) + 1) 

# =========== Step 3: Calculate Gene Variances =============
# Use the rowVars() function to calculate the variance for each gene across samples
# Gene variance is used to identify the most variable genes for further analysis
gene_variances <- rowVars(counts_matrix_log)

# =========== Step 4: Calculate Cumulative Variance =============
# Sort the variances in descending order and calculate the cumulative sum
cum_variance <- cumsum(sort(gene_variances, decreasing = TRUE)) / sum(gene_variances)

# =========== Step 5: Plot Cumulative Variance using ggplot2 for enhanced plotting =============
# Plot the cumulative proportion of variance explained as the number of genes increases
# This plot will help visualize how many top genes are needed to explain a high proportion of variance
ggplot(data.frame(Genes = 1:length(cum_variance), CumulativeVariance = cum_variance), aes(x = Genes, y = CumulativeVariance)) +
  geom_line() +
  labs(title = "Cumulative Variance Explained by Top Genes", 
       x = "Number of Genes", 
       y = "Cumulative Proportion of Variance Explained") +
  theme_minimal()
        

Code Block 3.2.5.1. Unsupervised Heatmap for All Samples (160 samples)

# ============ Load Required Libraries ============
# Load the required libraries for matrix operations and heatmap generation
library(matrixStats)  # Ensure matrixStats is loaded for row variance calculation
library(pheatmap)     # For generating heatmaps

# ============ Step 1: Convert counts_matrix to Matrix ============
# Convert counts_matrix to a matrix if it is currently a data frame
counts_matrix <- as.matrix(counts_matrix)  # Ensures compatibility with downstream functions

# ============ Step 2: Log2 Transformation of Counts ============
# Apply a log2 transformation to reduce the dynamic range of the count data
counts_matrix_log <- log2(counts_matrix + 1)  # Adding 1 to avoid log(0) issues

# ============ Step 3: Identify Top Variable Genes ============
# Recalculate the most variable genes based on the log-transformed counts
# We select the top 300 most variable genes using row variance
top_var_genes_log <- head(order(rowVars(counts_matrix_log), decreasing = TRUE), 300)  

# ============ Step 4: Generate Heatmap and Save ============
# Create the heatmap for all 166 samples based on the top 300 most variable genes
# Adjust clustering and visual parameters such as label angles and font sizes
pheatmap(counts_matrix_log[top_var_genes_log, ], 
         cluster_rows = TRUE, cluster_cols = TRUE, 
         show_rownames = FALSE, show_colnames = TRUE,  
         angle_col = 90,  # Rotate x-axis labels for better readability
         fontsize_col = 4,  # Reduce font size for the column labels
         filename = "Figures/3.2.5.Exploratory_Heatmaps/3.2.5.1.Unsupervised_Heatmap_160_Samples.png")  # Save the heatmap as an image
    

Code Block 3.2.5.2. Unsupervised Heatmap for Baseline Lesional and Non-lesional Samples (109 samples)

# ============ Load Required Libraries ============
# Load the required libraries for matrix operations and heatmap generation
library(matrixStats)  # Ensure matrixStats is loaded for row variance calculation
library(pheatmap)     # For generating heatmaps

# ============ Step 1: Subset Metadata for Lesional and Non-lesional Samples ============
# Subset the metadata to include only the baseline lesional (AL_m0) and non-lesional (AN_m0) samples
lesional_samples <- metadata[grepl("AL_m0", metadata$title), ]  # Filter metadata for lesional samples
non_lesional_samples <- metadata[grepl("AN_m0", metadata$title), ]  # Filter metadata for non-lesional samples

# Combine lesional and non-lesional samples into one dataset for comparison
selected_samples <- rbind(lesional_samples, non_lesional_samples)

# ============ Step 2: Subset the Count Matrix ============
# Subset the count matrix to include only the lesional and non-lesional samples
# Assuming the sample names (columns) in counts_matrix match those in metadata$title
selected_sample_names <- selected_samples$title
counts_matrix_filtered <- counts_matrix[, selected_sample_names]  # Filter the count matrix for the selected samples

# ============ Step 3: Log2 Transformation of Filtered Counts ============
# Apply a log2 transformation to the filtered count matrix to reduce the dynamic range
counts_matrix_log <- log2(as.matrix(counts_matrix_filtered) + 1)  # Adding 1 to avoid log(0) issues

# ============ Step 4: Identify Top Variable Genes ============
# Recalculate the most variable genes based on the log-transformed counts
# Select the top 300 most variable genes
top_var_genes_log <- head(order(rowVars(counts_matrix_log), decreasing = TRUE), 300)  

# ============ Step 5: Generate Heatmap and Save ============
# Create a heatmap for the lesional and non-lesional samples based on the top 300 most variable genes
# Adjust clustering and visual parameters such as label angles and font sizes
pheatmap(counts_matrix_log[top_var_genes_log, ], 
         cluster_rows = TRUE, cluster_cols = TRUE, 
         show_rownames = FALSE, show_colnames = TRUE,  
         angle_col = 90,  # Rotate x-axis labels for better readability
         fontsize_col = 4,  # Reduce font size for the column labels
         filename = "Figures/3.2.5.Exploratory_Heatmaps/3.2.5.2.Unsupervised_Heatmap_Baseline_Samples.png")  # Save the heatmap as an image
    

Code Block 3.2.5.3. Grouped Heatmap for All Groups

# =========== Load Required Libraries ===========
library(matrixStats)
library(pheatmap)

# =========== Step 1: Apply Log2 Transformation ===========
counts_matrix_log <- log2(as.matrix(counts_matrix) + 1)

# =========== Step 2: Extract Sample Grouping Information ===========
sample_groups <- metadata$source_name_ch1

# =========== Step 3: Define Descriptive Labels ===========
new_labels <- c(
  "m0_AN" = "Baseline Non-lesional",
  "cyclosporine_m3_AN" = "Cyclosporine Non-lesional",
  "dupilumab_m3_AN" = "Dupilumab Non-lesional",
  "m0_AL" = "Baseline Lesional",
  "cyclosporine_m3_AL" = "Cyclosporine Lesional",
  "dupilumab_m3_AL" = "Dupilumab Lesional"
)
	
# =================================================================================
# =========== Step 3.1: Select Groups for Analysis ===========
selected_groups <- c("m0_AN", "cyclosporine_m3_AN", "dupilumab_m3_AN", "m0_AL", "cyclosporine_m3_AL", "dupilumab_m3_AL")

# ===== we can replace the above selected_groups with one of the bellow  ===== 
# 1. Heatmap for all groups:
# c("m0_AN", "cyclosporine_m3_AN", "dupilumab_m3_AN", "m0_AL", "cyclosporine_m3_AL", "dupilumab_m3_AL")

# 2. (Lesional vs. Non-lesional) Compare "Baseline Lesional" with "Baseline Non-lesional": 
# c("m0_AN", "m0_AL")

# 3. (Lesional vs. Non-lesional) Compare Cyclosporine Lesional with Cyclosporine Non-lesional:
# c("cyclosporine_m3_AN", "cyclosporine_m3_AL") 

# 4. (Lesional vs. Non-lesional) Compare Dupilumab Lesional with Dupilumab Non-lesional:
# c("dupilumab_m3_AN", "dupilumab_m3_AL") 

# 5. (Therapy Comparisons) Compare Baseline Lesional, Cyclosporine Lesional, and Dupilumab Lesional (compared to Baseline Non-Lesional).
# c("m0_AN", "m0_AL", "cyclosporine_m3_AL", "dupilumab_m3_AL")  

# 6. (Therapy Comparisons) Compare Dupilumab Lesional and Cyclosporine Lesional.
# c("dupilumab_m3_AL", "cyclosporine_m3_AL")

# 7. (Pre-treatment vs. Post-treatment) Compare Baseline Lesional & Cyclosporine Lesional.
# c("m0_AL", "cyclosporine_m3_AL")

# 8. (Pre-treatment vs. Post-treatment) Compare Baseline Lesional & Dupilumab Lesional.
# c("m0_AL", "dupilumab_m3_AL")
# =================================================================================

	
# =========== Step 4: Reorder Columns Based on Group Selection ===========
desired_order <- selected_groups
reordered_columns <- match(desired_order, sample_groups)
counts_matrix_log <- counts_matrix_log[, reordered_columns]
colnames(counts_matrix_log) <- new_labels[sample_groups[reordered_columns]]

# =========== Step 5: Calculate Top Variable Genes ===========
top_var_genes_log <- head(order(rowVars(counts_matrix_log), decreasing = TRUE), 300)

# =========== Step 6: Generate the Heatmap ===========
pheatmap(counts_matrix_log[top_var_genes_log, ], 
         cluster_rows = TRUE, cluster_cols = TRUE,
         show_rownames = FALSE, show_colnames = TRUE,
         angle = 45, vjust = 0.5, hjust = 0.6,
         fontsize_col = 8, 
         filename = "Figures/3.2.5.Exploratory_Heatmaps/3.2.5.3.Grouped_Heatmap_All_Groups.png")
    

Code Block 3.2.5.4. Grouped Heatmap: Compare Baseline Lesional with Baseline Non-lesional

# ============================
selected_groups <- c("m0_AN", "m0_AL")
# ============================
    

Code Block 3.2.5.5. Grouped Heatmap: Compare Cyclosporine Lesional vs. Non-lesional

# ============================
selected_groups <- c("cyclosporine_m3_AN", "cyclosporine_m3_AL")
# ============================
    

Code Block 3.2.5.6. Grouped Heatmap: Compare Dupilumab Lesional vs. Dupilumab Non-lesional

# ============================
selected_groups <- c("dupilumab_m3_AN", "dupilumab_m3_AL")
# ============================
    

Code Block 3.2.5.7. Grouped Heatmap: Compare Baseline Non-Lesional vs. Baseline Lesional, Cyclosporine Lesional, and Dupilumab Lesional

# ============================
selected_groups <- c("m0_AN", "m0_AL", "cyclosporine_m3_AL", "dupilumab_m3_AL") 
# ============================
    

Code Block 3.2.5.8. Grouped Heatmap: Compare Dupilumab Lesional vs. Cyclosporine Lesional

# ============================
selected_groups <- c("dupilumab_m3_AL", "cyclosporine_m3_AL")
# ============================
    

Code Block 3.2.5.9. Grouped Heatmap: Compare Pre-treatment vs. Post-treatment (Baseline Lesional & Cyclosporine Lesional)

# ============================
selected_groups <- c("m0_AL", "cyclosporine_m3_AL")
# ============================
    

Code Block 3.2.5.10. Grouped Heatmap: Compare Pre-treatment vs. Post-treatment (Baseline Lesional & Dupilumab Lesional)

	# ============================
	selected_groups <- c("m0_AL", "dupilumab_m3_AL")
	# ============================
	    

Code Block 4. Filtering and Subsetting Samples for DESeq2 Analysis

				
# Step 4: Subset baseline lesional/non‑lesional (m0) and align with counts

# 4.1 keep only m0_AL / m0_AN
metadata_subset <- metadata[metadata$source_name_ch1 %in% c("m0_AL", "m0_AN"), ]

# 4.2 subset counts by titles (your counts colnames are titles)
keep <- intersect(colnames(counts_matrix), metadata_subset$title)
counts_subset <- counts_matrix[, keep, drop = FALSE]

# 4.3 reorder metadata to match counts columns
metadata_subset <- metadata_subset[match(colnames(counts_subset), metadata_subset$title), ]

# 4.4 condition factor
metadata_subset$condition <- ifelse(metadata_subset$source_name_ch1 == "m0_AL", "lesional", "non_lesional")
metadata_subset$condition <- factor(metadata_subset$condition, levels = c("non_lesional", "lesional"))

	
	
		    

Code Block 5.1.1. Creating the DESeq2 Dataset Object (dds)

# Step 5 – Create DESeq2 Dataset Object

# Step 5.1 — Assign row names to metadata
rownames(metadata_subset) <- metadata_subset$title

# Step 5.2 — Convert count matrix to integer mode	
counts_subset <- as.matrix(counts_subset)
storage.mode(counts_subset) <- "integer"

# Step 5.3 — Create the DESeq2 dataset object
dds <- DESeqDataSetFromMatrix(
  countData = counts_subset,
  colData   = metadata_subset,
  design    = ~ condition
)
  

6. Quality Control & Data Exploration (Post-DESeq2 Object Creation)

• This section assesses the quality and structure of the dataset after creating the DESeqDataSet object (dds). Earlier, in Section 3 (Exploratory Data Analysis — Before dds), we explored all samples to understand metadata structure, sample distributions, and the raw count matrix prior to modelling. Here, we shift focus to the selected analysis subset (baseline lesional and non-lesional samples) and perform post-dds quality checks to ensure the data are biologically consistent and free from technical artefacts or outliers before running DESeq().
• The analyses in this section begin with the Variance Stabilizing Transformation (VST), which converts raw count data into a stabilized scale suitable for visualization and exploratory analysis. Subsequent steps include generating a PCA plot to reveal overall clustering patterns, and creating a sample-to-sample distance heatmap to confirm these relationships numerically.
• Later in Section 8 (after differential expression analysis), a separate PCA and heatmap (DEG-focused PCA/heatmap) will be generated using only the most significantly differentially expressed genes (for example, the top 500 DEGs based on adjusted p-value). These DEG-focused plots emphasize the strongest condition-specific transcriptional changes and complement the global PCA by zooming in on the genes most responsible for those observed differences.
• The PCA and heatmap presented in section 6 are referred to as “global exploratory visualizations” because they use all genes in the dataset, typically after variance-stabilizing transformation (VST). These plots provide an overview of the dataset’s overall structure, helping identify potential batch effects, technical artifacts, or broad separation by biological condition.
• The section 6, finally, concludes with checks on library sizes and count distributions to validate normalization consistency and ensure that expression values are evenly distributed after transformation.
• Together, these quality control steps confirm that the dataset is technically reliable, biologically coherent, and properly normalized—providing a solid foundation for downstream differential expression analysis and pathway interpretation.

6.1.Variance Stabilizing Transformation (VST):

• VST transformation prepares data for visualization (Code Block 6.1). DESeq2 provides two transformation methods for this purpose: the variance stabilizing transformation (vst()) and the regularized logarithm transformation (rlog()). Both are designed to make the variance roughly independent of the mean across genes, which is essential for reliable visualizations and distance-based analyses. The VST is typically recommended for medium-to-large datasets (n > 30) due to its speed and efficiency, whereas the rlog transformation can be more suitable for smaller datasets (n < 30), offering slightly improved variance stabilization at the cost of higher computational time. These transformations are used solely for visualization and exploratory quality control; the actual differential expression analysis (DESeq()) is always performed on the raw count data.

Code Block 6.1. Variance Stabilizing Transformation (VST) Applied to DESeq2 Dataset

			  
# Apply variance stabilizing transformation (VST) to the DESeq2 dataset for downstream visualization
vsd <- vst(dds, blind=TRUE)

		  

6.2. Run PCA:

• Principal Component Analysis (PCA) is an essential exploratory step that helps reveal clustering and dominant sources of variation. It also helps visualize overall similarities and differences among samples based on their global gene expression profiles. By reducing the dimensionality of the variance-stabilized dataset, PCA projects each sample into a two-dimensional space defined by first two principal components the (PC1 and PC2), which capture the majority of variation in the data.
• In the resulting plot, as shown in Figure 6.2, samples are coloured according to their condition (lesional or non-lesional). A moderate separation can be observed along PC1, where non-lesional samples cluster predominantly on one side while lesional samples shift towards the other. However, some overlap remains, indicating natural biological variability within each group or potential subtle technical effects. This pattern confirms that condition explains a meaningful portion of the variance (26% on PC1 and 15% on PC2), while additional sources of variation contribute to finer sample differences.
• In PCA, each axis (principal component) represents a direction of maximum variation in the data, and the percentage of variance associated with it indicates how much of the overall gene expression variability that axis captures. Generally, when PC1 explains more than 40–50% of the variance, it suggests a strong and clear distinction between groups, often seen in datasets comparing highly distinct biological states such as different tissues or disease vs. healthy controls. Values between 20–40% usually reflect a moderate separation, typical for complex transcriptomic datasets, where biological differences exist but are not absolute. Values below 20% are considered weak separation, indicating subtle effects or highly homogeneous samples. In this dataset, PC1 (26%) and PC2 (15%) together explain 41% of the total variance, representing a meaningful yet moderate level of distinction between groups.
• The visual separation observed here aligns with what is expected for human RNA-seq data collected from matched samples. The partial clustering of lesional and non-lesional samples suggests that the condition is a major source of variation, while the remaining variability likely arises from patient-specific factors, biological heterogeneity, and minor technical differences. This balance between separation and overlap demonstrates that the transformation and normalization have effectively preserved biological structure while minimizing technical bias, confirming that the dataset is well-prepared for downstream differential expression analysis.
• The PCA plot is generated using the DESeq2 built-in function plotPCA(), which automatically computes the principal components from the variance-stabilized data and labels the axes with the percentage of total variance explained by each component.



Code Block 6.2. PCA Plot of Variance-Stabilized Data Grouped by Condition

 # ==============================================================
# 6.2. PCA of Variance-Stabilized Data (Global Overview)
# ==============================================================

library(ggplot2)
library(DESeq2)

# Extract PCA data from the variance-stabilized matrix (vsd)
pca_df <- plotPCA(vsd, intgroup = "condition", returnData = TRUE)
percentVar <- round(100 * attr(pca_df, "percentVar"))

# Build PCA plot manually for full control over aesthetics
p_pca_global <- ggplot(pca_df, aes(x = PC1, y = PC2, colour = condition)) +
  geom_point(size = 2.5, alpha = 0.85) +
  scale_colour_manual(values = c("lesional" = "#0072FF",
                                 "non_lesional" = "#66FF00")) +
  xlab(paste0("PC1: ", percentVar[1], "% variance")) +
  ylab(paste0("PC2: ", percentVar[2], "% variance")) +
  labs(title = "PCA of Variance-Stabilized Data (Global Gene Set)",
       colour = "Condition") +
  theme_bw(base_size = 14) +
  theme(
    plot.title = element_text(size = 18, face = "bold", hjust = 0.5),
    axis.title = element_text(size = 14, face = "bold"),
    axis.text = element_text(size = 12),
    legend.title = element_text(size = 13),
    legend.text = element_text(size = 11)
  )

# Save as high-resolution PNG
ggsave("C:/Users/YOUR/PATH/Figure_6.2.png",
       plot = p_pca_global, width = 10, height = 8, dpi = 300)

			  
		  

6.3. Identify possible batch effects:

• In RNA-seq analysis, a batch effect refers to any systematic technical variation unrelated to the biological question of interest. Such effects typically arise from differences in sequencing runs, library preparation dates, reagents, technicians, or instruments used. They can make samples appear more similar or different for technical rather than biological reasons, potentially biasing the results if not detected and corrected. Batch effects can be investigated using PCA patterns by colouring samples by batch, patient, or technical factor.
• In this study, all samples were generated within a single sequencing project (platform GPL21290) using a consistent experimental protocol, and the metadata does not indicate the presence of multiple batches or sequencing runs. Therefore, no explicit technical batch effect is expected. However, since each patient contributed both lesional and non-lesional samples collected at the same time, the primary source of variation may stem from biological pairing rather than technical processing. This represents a within-patient effect rather than a batch effect.
• As a result, it is not necessary to perform additional batch correction in this analysis. Nevertheless, if exploratory analyses such as PCA reveal clustering primarily by patient rather than by condition, the design formula in DESeq2 can be adjusted to include patient as a blocking factor (for example, ~ patient + condition) to control for within-patient variability.

6.4. Sample-to-Sample Distance Heatmap:

• In the resulting heatmap, as shown in Figure 6.4, each cell represents the Euclidean distance between a pair of samples, calculated using VST-transformed counts across all genes. The color gradient from blue (0) to red (150) reflects how similar or dissimilar two samples are in their overall gene expression profiles: blue indicates high similarity (low distance), while red denotes greater dissimilarity (high distance).



• The accompanying dendrograms display hierarchical clustering of samples based on these distances, allowing visual identification of groups with similar expression patterns. The diagonal line of dark blue squares corresponds to each sample compared to itself, producing a distance of zero, which naturally appears along the main diagonal of the matrix. The diagonal line is always deep blue (the smallest possible value) because the distance between a sample and itself is 0.
• In this dataset, several blue blocks along the diagonal suggest clusters of samples that are more alike — for instance, many lesional or non-lesional samples tend to cluster together — whereas red regions between clusters indicate samples that differ substantially, reflecting the biological contrast between conditions. This structure mirrors the moderate separation observed in the PCA plot, confirming that condition is a key contributor to the variation while additional factors such as inter-patient variability explain finer differences.
• Overall, the structure of the heatmap supports clear separation of lesional and non-lesional skin which confirms the dataset’s technical soundness and biological coherence, ensuring that it is well-suited for downstream differential expression analysis.

Code Block 6.4. Generate and Visualize Sample-to-Sample Distance Heatmap Using VST-Transformed Counts

# ==============================================================
# 6.4. Sample-to-Sample Distance Heatmap
# ==============================================================

# ========= Load Required Package =========
library(pheatmap)

# ========= Compute Pairwise Sample Distances =========
# Calculate distances between samples using VST-transformed counts.
# Transpose the assay matrix so that distance is computed across samples.
dists <- dist(t(assay(vsd)))

# Convert distance object to a full matrix for visualization
mat <- as.matrix(dists)

# ========= Generate and Save Heatmap =========
# Create a sample-to-sample distance heatmap and save as a high-resolution PNG.
# This heatmap helps visualize overall similarity patterns between samples.
pheatmap(
  mat,
  clustering_distance_rows = dists,  # hierarchical clustering for rows
  clustering_distance_cols = dists,  # hierarchical clustering for columns
  main = "Sample-to-sample distances",  # heatmap title
  color = colorRampPalette(c("navy", "white", "firebrick3"))(50),  # blue-white-red color scale
  
  # ========= Font and Label Settings =========
  fontsize = 28,       # controls general font size (including title)
  fontsize_row = 6,    # smaller font for row names
  fontsize_col = 6,    # smaller font for column names
  angle_col = 90,      # rotate column labels for better visibility
  
  # ========= Dendrogram and Legend Settings =========
  treeheight_row = 80, # row dendrogram height
  treeheight_col = 80, # column dendrogram height
  legend = TRUE,       # display color legend
  
  # ========= Display and Output Settings =========
  show_rownames = TRUE,  # show sample names on rows
  show_colnames = TRUE,  # show sample names on columns
  filename = "C:/Users/omidm/OneDrive/Desktop/RNA_Seq_Task1/6.EDA/Figure_6.4.png",  # output path
  width = 18, height = 16  # figure size in inches for crisp rendering
) 
			  


			  
		  

6.5. Check library sizes:

• The Library Size Check plot confirms normalization and sequencing depth uniformity and visualizes the total number of reads (counts) obtained for each sample after alignment and counting. Each bar represents one sample, and its height reflects the sequencing depth — the total number of reads mapped to genes in that sample. The y-axis indicates the total number of reads, and in this dataset, ranging from approximately 13 million to 32 million per sample, while the x-axis lists the samples grouped by condition (baseline lesional and non-lesional). The bars are color-coded to distinguish between conditions (lesional = blue , non-lesional = yellow). This visualization provides an immediate overview of sequencing consistency across all samples.
• Ideally, RNA-seq libraries should have comparable sequencing depths across all samples, as major deviations can bias normalization and affect downstream analyses. In Figure 6.5, most samples fall between 20 and 32 million reads, indicating a balanced and uniform sequencing effort. A few samples toward the left exhibit smaller library sizes (around 13–16 million reads), but these are not drastically lower and therefore remain within an acceptable range. This suggests that sequencing was consistent overall and no major technical issues occurred during library preparation.

Figure 6.5. Library Size Check per Sample. Each bar represents the total number of mapped reads (library size) for one sample after alignment and counting. Bars are colored by condition (lesional = blue, non-lesional = yellow). The y-axis shows total read counts ranging from approximately 13 to 32 million reads per sample, illustrating consistent sequencing depth across all libraries. No extreme deviations are observed, confirming balanced sequencing coverage suitable for reliable normalization and downstream differential expression analysis.
• Uniform library sizes are essential because methods such as DESeq2 assume similar sequencing depths across samples when performing normalization. Large discrepancies (e.g., one sample with 5M reads and another with 80M) could bias normalization and lead to inaccurate estimation of differential expression results. Ensuring that all libraries fall within a comparable range helps maintain fair comparisons between conditions and supports robust statistical inference.
• From a quality control standpoint, it is important to monitor both the lower and upper extremes of library sizes. Libraries with very low read counts (e.g., below 10 million) often indicate technical problems such as poor RNA quality, failed library preparation, or insufficient sequencing depth. Conversely, extremely large libraries (above 60–70 million reads) can disproportionately influence normalization, masking subtle expression differences among the rest of the samples. In this dataset, however, all libraries fall comfortably within a healthy range (13–32 million reads), reflecting excellent sequencing quality.
• As illustrated earlier in Figure 2.2.1.2, six outlier samples were identified and removed during preprocessing — three with unusually low library sizes (approximately 10 million reads) and three with abnormally high library sizes. Although high-count samples can sometimes represent biologically interesting cases, they can also distort global normalization and obscure population-wide trends. Removing both the low- and high-end outliers improved the uniformity of the dataset, ensuring that subsequent analyses reflect the core biological signal rather than technical variation.
• This diagnostic visualization confirms that no major sequencing depth imbalance exists, no sample dominates or underperforms significantly, and the dataset is technically robust. The resulting distribution of library sizes demonstrates high data integrity and supports reliable normalization, making the dataset ready for downstream steps such as differential expression analysis and functional exploration.

Code Block 6.5. Visualizing Library Sizes Across Samples Using ggplot2

											  
# ==============================================================
# 6.5. Library Size Check (Enhanced Visualization)
# ==============================================================

library(ggplot2)

# Extract library sizes (total read counts per sample)
library_sizes <- colSums(counts(dds))

# Create a data frame for plotting
library_df <- data.frame(
  Sample = names(library_sizes),
  LibrarySize = library_sizes,
  Condition = colData(dds)$condition  # assumes condition column exists
)

# Order samples by library size for visual clarity
library_df$Sample <- factor(library_df$Sample, levels = library_df$Sample[order(library_df$LibrarySize)])

# Generate the bar plot
ggplot(library_df, aes(x = Sample, y = LibrarySize, fill = Condition)) +
  geom_bar(stat = "identity", color = "black", width = 0.9) +
   geom_text(aes(label = scales::comma(LibrarySize)), 
            vjust = 0.4, hjust= 0,size = 3.5, angle = 90, fontface = "bold") +  # add labels
  scale_fill_manual(values = c("lesional" = "#0072FF", "non_lesional" = "#FFFF00")) +
  theme_bw() +
  theme(
    axis.text.x = element_text(angle = 90, vjust = 0.5, hjust = 1, size = 14),
    axis.text.y = element_text(size = 14),
    axis.title = element_text(size = 24, face = "bold"),
    plot.title = element_text(size = 34, face = "bold", hjust = 0.5),
    legend.title = element_text(size = 24),
    legend.text = element_text(size = 14)
  ) +
  labs(
    title = "Library Sizes per Sample",
    x = "Samples",
    y = "Total Read Counts",
    fill = "Condition"
  )

# Optional: Save as a high-resolution PNG
ggsave("C:/Users/omidm/OneDrive/Desktop/RNA_Seq_Task1/6.EDA/Figure_6.5.png",
       width = 25, height = 8, dpi = 300)
											  
										  

6.6. Count Distribution — Pre- and Post-Normalization.

Inspecting count distributions with two complementary approaches:

• To diagnose data quality after creating the DESeq2 object, it is useful to visualise how per-sample expression values are distributed both before and after normalisation. The first approach applies a simple log₂ transform to the raw counts, log2(counts(dds) + 1), which compresses large values but still reflects differences in sequencing depth and composition. The second approach uses the DESeq2 variance stabilising transformation (VST), assay(vsd), which removes the mean–variance relationship and places samples on a common, log2-like scale that is more appropriate for distance-based exploratory analysis. Together, these plots provide a before/after view that helps confirm normalisation has behaved as expected and that no samples have atypical distributions.



How each view helps:

• The pre-normalisation boxes in Figure 6.6.1 typically show slightly different medians and spreads across samples because technical factors (library size and composition) are still present. This is a healthy sign at this stage and can reveal obviously under-sequenced or over-sequenced libraries if they exist. The post-normalisation boxes in Figure 6.6.2 should align closely across samples, with similar medians and interquartile ranges. That tighter alignment indicates that size-factor normalisation and VST have removed unwanted technical variability while preserving the underlying biological signal for PCA and downstream differential testing.

Interpreting the results in this dataset:

• The pre-normalisation plot shows modest variability in medians and spreads, which is expected given natural differences in sequencing depth. After VST, the boxplots cluster tightly with a shared central tendency, and the global median line falls sensibly through the middle of the boxes. This pattern confirms that the dataset is well normalised and suitable for subsequent steps such as PCA, sample-to-sample distances, and DE analysis. Note that VST does not aim to separate lesional from non-lesional groups at this stage; it aims to equalise technical variation so that group structure can be evaluated reliably by PCA and DEA rather than by raw distribution shape.
• 

  Figure 6.6.1. Pre-normalisation count distributions using log2(counts(dds) + 1). Boxes reflect per-sample distributions of log-transformed raw counts. Modest shifts in medians and box widths indicate residual technical variation prior to normalisation.
• 

  Figure 6.6.2. Post-normalisation count distributions using VST (assay(vsd)). Boxes align closely across samples, demonstrating successful removal of mean–variance dependence and effective size-factor normalisation. The dashed line marks the global median VST value.


Below is a step-by-step code walkthrough for generating both figures:
• Bring in the plotting and reshaping toolkits: We load ggplot2 for building publication-quality figures and reshape2 for converting matrices/data frames into a long format that ggplot understands. This keeps the plotting pipeline tidy and consistent across both approaches.

  
  library(ggplot2)
  library(reshape2)
  	

• Start from the raw material: We extract the integer count matrix from the DESeqDataSet. At this stage nothing is normalised; the values reflect the number of reads assigned to each gene per sample.

  
  raw_counts <- counts(dds)
  	

• Compress the dynamic range for visibility: Taking log base 2 of the counts makes large values manageable and separates low values enough to see structure. Adding 1 prevents problems where a gene has zero counts in a sample. This is a visual convenience, not a normalisation step.

  
  log_counts <- log2(raw_counts + 1)
  	

• Switch to a data-frame workflow and retain gene IDs: ggplot works most smoothly with data frames. We convert the matrix and preserve the gene names in a dedicated column so they are not lost when reshaping.

  
  log_counts_df <- as.data.frame(log_counts)
  log_counts_df$Gene <- rownames(log_counts_df)
  	

• Reshape to long format for sample-wise boxplots: Each row becomes one gene–sample measurement. The sample name moves into a variable column, and the transformed expression values live in a single numeric column. This is the grammar ggplot expects for faceted or grouped summaries.

  
  log_counts_long <- melt(
    log_counts_df,
    id.vars = "Gene",
    variable.name = "Sample",
    value.name = "Log2Counts"
  )
  	

• Attach biological labels for colouring: We join each sample to its condition (lesional or non_lesional) using a safe match against the column metadata. This enables clear colour grouping in the boxplots without changing the underlying data.

  
  log_counts_long$Condition <- colData(dds)$condition[
    match(log_counts_long$Sample, rownames(colData(dds)))
  ]
  	

• Build the pre-normalisation (Raw Count Distribution) boxplot using log₂(counts + 1): We plot the log₂ counts per sample with boxes coloured by condition. Outliers are shown lightly to reduce clutter on wide plots. Labels and theme choices aim for legibility with many samples.

  
  p1 <- ggplot(log_counts_long, aes(x = Sample, y = Log2Counts, fill = Condition)) +
    geom_boxplot(outlier.alpha = 0.1, lwd = 0.3) +
    scale_fill_manual(values = c("lesional" = "#0072FF", "non_lesional" = "#FFFF00")) +
    labs(
      title = "Pre-normalization Count Distributions (log2(counts + 1))",
      x = "Samples",
      y = "Log2 Expression (Raw Counts)",
      fill = "Condition"
    ) +
    theme_bw() +
    theme(
      plot.title  = element_text(size = 16, face = "bold", hjust = 0.5),
      axis.text.x = element_text(angle = 90, vjust = 0.5, hjust = 1, size = 6),
      axis.title  = element_text(size = 12, face = "bold"),
      legend.title= element_text(size = 10),
      legend.text = element_text(size = 9)
    )
  	

• Save a crisp, wide image: We write the pre-normalisation panel to disk at high DPI to keep labels sharp on a wide canvas.

  
  ggsave("C:/Users/YOUR/PATH/Figure_6.6.1.png",
         plot = p1, width = 18, height = 8, dpi = 300)
  	

• Move to the normalised scale (VST): while we have already created vsd with vst(dds, blind = TRUE), we pull the variance-stabilised matrix. VST removes library-size effects and stabilises variance across the mean, making samples comparable.

  
  vst_mat <- assay(vsd)
  	

• Prepare and reshape the VST matrix for plotting: As before, we convert to a data frame, preserve gene IDs, and melt to a long layout. The resulting table mirrors the pre-normalisation data shape, which simplifies plotting and comparison.

  
  vst_df <- as.data.frame(vst_mat)
  vst_df$Gene <- rownames(vst_df)
  vst_long <- melt(
    vst_df,
    id.vars = "Gene",
    variable.name = "Sample",
    value.name = "VST"
  )
  	

• Reapply condition labels and impose a meaningful sample order: We colour by condition, then compute each sample’s median VST and order the x-axis by this statistic. Ordering reduces visual noise and helps we scan for global shifts.

  
  vst_long$Condition <- colData(dds)$condition[
    match(vst_long$Sample, rownames(colData(dds)))
  ]
  sample_order <- aggregate(VST ~ Sample, data = vst_long, median)
  vst_long$Sample <- factor(vst_long$Sample,
                            levels = sample_order$Sample[order(sample_order$VST)])
  	

• Add a global reference line: We compute the overall median VST across all samples. A dashed horizontal line helps us visually confirm that medians sit on a common scale after normalisation.

  
  global_median <- median(vst_long$VST, na.rm = TRUE)
  	

• Build the post-normalisation boxplot using VST (Variance Stabilizing Transformation): With VST values, boxes should align closely if normalisation succeeded. The dashed line is the global median for quick reference.

  
  p2 <- ggplot(vst_long, aes(x = Sample, y = VST, fill = Condition)) +
    geom_boxplot(outlier.alpha = 0.1, lwd = 0.3) +
    geom_hline(yintercept = global_median, linetype = "dashed", linewidth = 0.4, alpha = 0.6) +
    scale_fill_manual(values = c("lesional" = "#0072FF", "non_lesional" = "#FFFF00")) +
    labs(
      title = "Post-normalization Count Distributions (VST)",
      x = "Samples (ordered by median VST)",
      y = "VST Expression (log2-like Units)",
      fill = "Condition"
    ) +
    theme_bw() +
    theme(
      plot.title  = element_text(size = 16, face = "bold", hjust = 0.5),
      axis.text.x = element_text(angle = 90, vjust = 0.5, hjust = 1, size = 6),
      axis.title  = element_text(size = 12, face = "bold"),
      legend.title= element_text(size = 10),
      legend.text = element_text(size = 9)
    )
  	

• Export the VST panel: Same geometry as the first figure so the two views are directly comparable on the page.

  
  ggsave("C:/Users/YOUR/PATH/Figure_6.6.2.png",
         plot = p2, width = 18, height = 8, dpi = 300)
  	

Code Block 6.6.Inspecting Count Distributions — Comparing Raw and VST-Transformed Data

											  
# ==============================================================
# 6.6. Count Distribution — Pre- and Post-Normalization Comparison
# ==============================================================

library(ggplot2)
library(reshape2)

# -------------------------
# Step 1: Raw log2(counts + 1)
# -------------------------

# Extract raw counts from DESeqDataSet
raw_counts <- counts(dds)

# Log2 transform with +1 to avoid log(0)
log_counts <- log2(raw_counts + 1)

# Convert to data frame (samples as columns)
log_counts_df <- as.data.frame(log_counts)
log_counts_df$Gene <- rownames(log_counts_df)

# Melt to long format for ggplot2
log_counts_long <- melt(log_counts_df, id.vars = "Gene",
                        variable.name = "Sample", value.name = "Log2Counts")

# Add condition information
log_counts_long$Condition <- colData(dds)$condition[match(log_counts_long$Sample, rownames(colData(dds)))]

# Create pre-normalization boxplot
p1 <- ggplot(log_counts_long, aes(x = Sample, y = Log2Counts, fill = Condition)) +
  geom_boxplot(outlier.alpha = 0.1, lwd = 0.3) +
  scale_fill_manual(values = c("lesional" = "#0072FF", "non_lesional" = "#FFFF00")) +
  labs(
    title = "Pre-normalization Count Distributions (log2(counts + 1))",
    x = "Samples",
    y = "Log2 Expression (Raw Counts)",
    fill = "Condition"
  ) +
  theme_bw() +
  theme(
    plot.title  = element_text(size = 16, face = "bold", hjust = 0.5),
    axis.text.x = element_text(angle = 90, vjust = 0.5, hjust = 1, size = 6),
    axis.title  = element_text(size = 12, face = "bold"),
    legend.title= element_text(size = 10),
    legend.text = element_text(size = 9)
  )

# Save pre-normalization plot
ggsave("C:/Users/YOUR/PATH/ / / /Figure_6.6.1_PreNormalization.png",
       plot = p1, width = 18, height = 8, dpi = 300)


# -------------------------
# Step 2: Post-normalization (VST)
# -------------------------

# Compute variance-stabilized data if not already done
# vsd <- vst(dds, blind = TRUE)

vst_mat <- assay(vsd)

# Convert to data frame for ggplot2
vst_df <- as.data.frame(vst_mat)
vst_df$Gene <- rownames(vst_df)

# Melt to long format
vst_long <- melt(vst_df, id.vars = "Gene",
                 variable.name = "Sample", value.name = "VST")

# Add condition info
vst_long$Condition <- colData(dds)$condition[match(vst_long$Sample, rownames(colData(dds)))]

# Compute median per sample for ordering
sample_order <- aggregate(VST ~ Sample, data = vst_long, median)
vst_long$Sample <- factor(vst_long$Sample, levels = sample_order$Sample[order(sample_order$VST)])

# Global median line
global_median <- median(vst_long$VST, na.rm = TRUE)

# Create post-normalization boxplot
p2 <- ggplot(vst_long, aes(x = Sample, y = VST, fill = Condition)) +
  geom_boxplot(outlier.alpha = 0.1, lwd = 0.3) +
  geom_hline(yintercept = global_median, linetype = "dashed", linewidth = 0.4, alpha = 0.6) +
  scale_fill_manual(values = c("lesional" = "#0072FF", "non_lesional" = "#FFFF00")) +
  labs(
    title = "Post-normalization Count Distributions (VST)",
    x = "Samples (ordered by median VST)",
    y = "VST Expression (log2-like Units)",
    fill = "Condition"
  ) +
  theme_bw() +
  theme(
    plot.title  = element_text(size = 16, face = "bold", hjust = 0.5),
    axis.text.x = element_text(angle = 90, vjust = 0.5, hjust = 1, size = 6),
    axis.title  = element_text(size = 12, face = "bold"),
    legend.title= element_text(size = 10),
    legend.text = element_text(size = 9)
  )

# Save post-normalization plot
ggsave("C:/Users/YOUR/PATHS/Figure_6.6.2_PostNormalization.png",
       plot = p2, width = 18, height = 8, dpi = 300)

7. Differential Expression Analysis:

7.1. Fit the DESeq2 model

• This step runs DESeq() on the prepared dds object. DESeq2 estimates size factors, dispersions, and fits a negative binomial GLM per gene based on the design (~ condition). The fitted model is then used for statistical testing in the next steps.

7.2. Extract results for lesional vs non_lesional and filter significant genes

• We extract the contrast lesional vs non_lesional, which yields log2 fold changes interpreted as lesional / non_lesional. We set common thresholds for significance, for example adjusted P value < 0.05 and absolute log2 fold change ≥ 1. The summary shows overall test diagnostics, and we count how many genes pass the thresholds.

7.3. Shrink log2 fold changes with apeglm for more stable effect sizes

• Log2 fold changes from the raw Wald test can be noisy for lowly expressed genes. apeglm shrinkage stabilizes effect sizes while keeping direction and relative ranking useful for interpretation and plotting. We identify the correct coefficient name from resultsNames(dds), apply lfcShrink(), then reapply the same significance thresholds to obtain a shrunk set of DEGs.

8. DEG Interpretation, Visualization, and Export

• After completing DESeq2 model fitting and log2 fold change shrinkage (Section 7), this section focuses on the interpretation and visualization of differential expression results. The purpose is to move beyond raw statistics and uncover biologically meaningful patterns among the identified genes. We begin by summarizing the number of significantly up- and down-regulated genes and visualizing their overall distribution using summary bar plots.
• Next, we explore the structure of the dataset using DEG-focused PCA and sample-to-sample distance heatmap. Unlike the global exploratory plots in Section 6 (which included all genes to assess general structure and batch effects), these DEG-focused visualizations emphasize only the top differentially expressed genes—those driving the most distinct transcriptional differences between lesional and non-lesional skin. This provides a clearer view of how samples cluster once only the most biologically responsive genes are considered.
• Subsequent subsections present DEG-level visualizations including the Volcano plot, MA plot, heatmap of top DEGs, and gene expression density plots. Together, these highlight the direction, magnitude, and consistency of differential expression across genes and conditions. Finally, the complete and filtered DEG tables (including human-readable gene symbols) are exported for downstream analyses such as pathway enrichment, network mapping, and biomarker discovery.

8.0. Setup & Common Helpers

• This short setup chunk is separated from the rest so everything the later plots and tables depend on is defined once, in one place. It loads the packages we actually call in Section 8; sets a single output folder and a single colour palette so all figures look consistent; declares the analysis thresholds used across the section (adjusted P value and |log2FC|); constructs a tidy results data frame from the shrunken DESeq2 results; brings the VST expression matrix into scope for plotting; and prepares helper functions to map Ensembl gene IDs to human gene symbols. Keeping these shared pieces at the top prevents duplication and makes it easy to change a setting (for example, the colour for “lesional” or the |log2FC| cut-off) without chasing it through multiple blocks.
• The libraries are loaded inside suppressPackageStartupMessages({ ... }) to keep the console clean. DESeq2 supplies assay(vsd) and the shrunken results object we created in Section 7. ggplot2 is used for the PCA, volcano, and density figures. pheatmap builds the heatmaps; RColorBrewer provides palettes for them. reshape2 contributes melt() when we reshape matrices to long format for ggplot. dplyr powers the data-wrangling verbs (mutate(), filter(), case_when(), relocate()). readr writes the CSVs. Finally, AnnotationDbi + org.Hs.eg.db perform the Ensembl→HGNC symbol mapping used in plot labels and exported tables.
• out_dir defines where Section-8 outputs are saved and is created if missing; cond_cols fixes a single, consistent mapping from condition values to colours. The shared thresholds alpha (adjusted P value) and lfc_cut (absolute log2 fold change) are declared once and reused everywhere so filtering logic is consistent across the section.
• res_df is the tidy data frame most plots read from. We start by coercing the shrunken DESeq2 results (res_shrunk) to a data frame, then add a gene column from the row names (mutate(gene = rownames(...))). We remove rows without an adjusted P value (filter(!is.na(padj))). Next, case_when() creates a simple categorical status: “Up” if padj < alpha and log2FoldChange >= lfc_cut; “Down” if padj < alpha and log2FoldChange <= -lfc_cut; otherwise “NS” (not significant). This three-way label drives the volcano plot colouring and helps with quick summaries.
• vsd_mat <- assay(vsd) pulls in the variance-stabilised expression matrix (rows = genes, columns = samples) prepared earlier; it’s the common input for the DEG-focused PCA, distance heatmaps, and expression density panels. For human-readable labels, the mapping helpers strip Ensembl version suffixes (strip_ver() turns ENSG000001234.5 into ENSG000001234) and then look up symbols with AnnotationDbi::mapIds() from org.Hs.eg.db. The wrapper ens_to_symbol() converts any vector of Ensembl IDs to symbols, falling back to the original Ensembl ID if a symbol is unavailable. Finally, we append a gene_symbol column to res_df so symbols are available to every downstream plot and export without repeating the mapping step.

8.1. Overview and Summary of DEG Results

• After performing differential expression analysis and log2 fold-change shrinkage, we identified significant genes using two main thresholds: an adjusted p-value (< 0.05) and an absolute log2 fold change (≥ 1). This filtering step highlights genes showing biologically meaningful and statistically reliable changes between lesional and non-lesional skin.
• This section provides a concise overview and summary of the differential expression results generated by DESeq2. Before diving into visual exploration, it is essential to quantify how many genes are significantly up-regulated or down-regulated after model fitting and shrinkage. These simple counts give an immediate sense of the overall transcriptional activity—whether the dataset reflects widespread activation or repression of gene expression under the studied condition.
• The code begins by calculating basic DEG statistics from res_df, the processed results data frame created earlier. Using sum() in combination with logical conditions, it computes:
• deg_total – the total number of significant DEGs (res_df$sig != "NS").
• deg_up – the number of up-regulated genes (res_df$sig == "Up").
• deg_down – the number of down-regulated genes (res_df$sig == "Down").
• The cat() command then prints this summary directly to the console, providing a quick textual check of the number of DEGs before moving to graphical representation.
• Next, the results are visualized in a simple yet informative bar plot summary (Figure 8.1). The code constructs a small data frame called deg_df that stores the directions (“Up-regulated” and “Down-regulated”) along with their respective counts. The ggplot2 function geom_col() is used to generate vertical bars whose heights correspond to the number of genes in each category. The scale_fill_manual() line assigns fixed colours to the two groups—pink (#EF476F) for up-regulated and blue (#118AB2) for down-regulated genes—ensuring visual consistency across all figures.



• The theme_bw() function applies a clean white background for clarity, and the theme() call centers the title while adjusting font sizes for readability. The resulting plot, titled “Significant DEGs (after shrinkage)”, provides a quick visual impression of the dataset’s directionality—whether up-regulation dominates, down-regulation dominates, or the balance is roughly equal. Finally, ggsave() exports the figure as a high-resolution PNG file (Figure_8.1_DEG_Summary_Bar.png) to the designated out_dir folder for documentation and reporting.

8.2. DEG-Focused Diagnostic Visualizations (QC & Data Patterns)

• In Section 6, we generated global PCA and sample-to-sample distance heatmaps using all VST-transformed genes to assess overall dataset structure, normalization, and possible batch effects. In this section, we shift to a DEG-focused view by restricting the analysis to the top differentially expressed genes (ranked by adjusted p-value after LFC shrinkage). These genes are the main contributors to the iological separationb between lesional and non-lesional samples. This approach highlights patterns driven specifically by disease-associated genes and helps confirm that the identified DEGs form a biologically coherent structure across samples.

8.2.1 PCA Plot on Top DEGs:

• The Principal Component Analysis (PCA) shown in Figure 8.2.1 was computed using only the top 500 most significant DEGs (ranked by adjusted p-value after shrinkage). Before performing PCA, the gene expression matrix underwent a variance stabilizing transformation (VST) to make variance approximately independent of the mean, followed by row-wise scaling (z-score normalization of each gene). This scaling step ensures that all genes contribute equally to the principal components, regardless of their absolute expression levels — otherwise, highly expressed genes could dominate the analysis. Row-wise scaling means that each row (gene) in the expression matrix is standardized independently — that is, its expression values across all samples are rescaled so that the gene’s mean becomes 0 and its standard deviation becomes 1.
• Each point in the PCA plot represents one sample, and its coordinates (PC1 and PC2) summarize the overall expression pattern of these top DEGs. PC1 captures the largest source of variation (≈ 61%) and PC2 the second largest (≈ 6%). The clear separation along PC1 between lesional (blue) and non_lesional (light green) samples indicates that the expression of these 500 DEGs alone is sufficient to distinguish the two biological states with minimal overlap. This strong separation demonstrates that the DESeq2 model successfully identified genes whose combined expression profiles consistently encode the disease signature of atopic dermatitis lesions.
• Biologically, this pattern means that most lesional samples share a coherent transcriptional program distinct from non-lesional skin, likely reflecting the activation of inflammatory, immune-response, and barrier-disruption pathways. The limited overlap between groups suggests that inter-patient variability exists but does not obscure the main disease signal. From a data-quality perspective, such a distinct clustering also implies that normalization and model fitting were successful: the primary source of variance in the dataset corresponds to true biological differences rather than technical artifacts. PCA at this stage therefore serves both as a biological confirmation of DEG relevance and a diagnostic quality-control step, showing that the identified DEGs are robust markers separating lesional and non-lesional transcriptomes.

8.2.2 Using PC1 Loadings to Identify Disease-Relevant Driver Genes

• PCA (Principal Component Analysis) identifies the principal components (PCs) that explain the greatest variance in the data. In this analysis, PCA on the top differentially expressed genes (DEGs) provided a strong separation (61%) between lesional and non-lesional samples, with PC1 explaining of the variance (see Figure 8.2.1). Because this PCA was fit only on the top 500 DEGs, the variance captured by PC1 is dominantly biological rather than technical. The natural next step is to exploit PC1 to find the specific genes that drive this separation.
• In PCA, applied to gene expression data, each principal component is a weighted combination of gene expression variables. Those weights are called loadings, and in R’s prcomp() they are stored in the rotation matrix (pca$rotation). For PC1, the loading of a gene tells us both how strongly and in which direction that gene influences the first principal component. High-magnitude loadings mark the genes that most define PC1 and, therefore, the main axis that separates lesional from non-lesional samples in our DEG-focused space.
• Concretely, the PC1 loadings are the numerical coefficients showing how much each gene “loads onto” PC1. A larger absolute loading means a stronger contribution to that component. The loadings are unit-normalized per component: the sum of squared loadings (∑loading²) across genes equals 1 for each PC. The raw sum of loadings does not carry a constraint or interpretation.
• The sign of a loading indicates the direction relative to the sample scores. In our analysis, lesional samples cluster on the negative side of PC1, and non-lesional on the positive side. That implies genes with negative PC1 loadings are lesional-enriched drivers, while genes with positive PC1 loadings are non-lesional-enriched. This sign convention is not universal: PCA axes can flip with no change to the model’s meaning. We verified the direction by comparing mean PC1 scores across groups and only then mapped negative loadings to lesional enrichment and positive loadings to non-lesional enrichment.
• With that orientation in hand, we extract three gene sets: 1️⃣ all PC1-contributing genes (See “All” CSV; ⬇️ Download CSV ) ranked by absolute loading (the comprehensive view), 2️⃣ the lesional-enriched set (negative loadings) (See “Lesional-enriched” CSV; ⬇️ Download CSV), and 3️⃣ the non-lesional-enriched set (positive loadings) (See “Non-lesional-enriched” CSV; ⬇️ Download CSV). Sorting by loading magnitude within each sign gives a biologically ordered list of candidate drivers: the most negative loadings likely represent disease-state drivers in lesions, while the most positive loadings reflect homeostatic or non-lesional programs.
• When molecular profiles from treated samples are available, they can be compared against these PCA-derived disease and healthy gene lists to determine which treatments most effectively shift expression patterns toward the non-diseased state. This can be done by measuring the degree of overlap, similarity, or enrichment between treatment-responsive genes and the reference signatures, allowing for an objective and biologically grounded assessment of therapeutic effectiveness.

• Code Block 8.2.2.1 selects the top 500 DEGs based on adjusted p-value, subsets the variance-stabilized expression matrix to those genes, and performs PCA on the transposed matrix using prcomp(t(vst_top), scale.=TRUE). The PC1 loadings are extracted from the rotation matrix (pca$rotation[, 1]), and a tidy data frame is created containing Ensembl IDs, corresponding gene symbols, and both the raw and absolute loading values. The genes are then sorted by the magnitude of their loadings, so that the most influential contributors to PC1 appear at the top. The resulting table provides a comprehensive, ranked overview of all genes contributing to PC1 and is saved as a CSV file for downstream analysis. This dataset serves as the foundation for distinguishing between lesional- and non-lesional-enriched subsets and for conducting subsequent functional or pathway analyses focused on the most biologically impactful genes.

• Code Block 8.2.2.2 builds upon the complete PC1 table and isolates genes with negative loading values, which correspond to lesional-enriched profiles based on the PCA orientation. The genes are sorted in ascending order of PC1 loading so that those with the most negative loadings—hence the strongest association with lesional samples—appear first. The filtered dataset is then saved as a CSV file containing all lesional-enriched genes. These genes represent the key drivers of the lesional phenotype, capturing transcriptional programs that characterize the disease state. They form an essential resource for downstream pathway analysis, identification of disease-related signatures, and prioritization of biomarkers potentially relevant to lesional tissue biology.

• Code Block 8.2.2.3 mirrors the previous one but focuses on genes with positive PC1 loadings, representing those associated with non-lesional or homeostatic expression patterns. By sorting the data in descending order of PC1 loading, the genes most positively influencing PC1 are placed at the top. The resulting list is saved as a CSV file containing all non-lesional-enriched genes. This gene set provides a counterpart to the lesional-enriched list, highlighting biological processes that define the healthy or non-diseased state. Comparing these two sets of genes—those driving the lesional and non-lesional ends of PC1—offers a balanced and biologically interpretable framework for understanding how transcriptional programs diverge between disease and normal skin.

• Together, these three analyses provide a comprehensive view of PC1-driven biology, where the full loading table captures the global structure of variation, and the sign-specific subsets reveal directionally distinct gene groups underpinning the separation between lesional and non-lesional samples observed in Figure 8.2.1.

8.2.3 Sample-to-Sample Distance Heatmap on Top DEGs:

• 8.2.3 Sample-to-Sample Distance Heatmap on Top DEGs — Figure 8.2.3 displays a Euclidean distance heatmap based on the VST expression of the top DEGs after row-wise z-score normalization. In this visualization, darker navy areas indicate higher similarity between samples, while red and white regions represent greater distances (lower similarity). Two well-defined low-distance blocks align with the lesional and non-lesional groups, confirming that DEGs capture a clear and consistent condition-driven separation. This distinct grouping, also illustrated by the dendrograms along the axes, confirms that the identified top DEGs are biologically meaningful and effectively discriminate between the two different condition groups (lesional and non-lesional skin), thereby validating the strength of the differential expression analysis and the coherence of the disease-associated gene signature across samples.
• Data-Driven Sample Clustering for Stratified Differential Analysis: In future analyses, a useful approach is to separate these identified clusters—where, for example, non-lesional and lesional samples form their own natural groups—and treat them as distinct stratified categories. By doing this, each cluster can be analyzed separately for differential expression and functional insights. This stratified method helps uncover unique molecular signatures for each subgroup, allowing for meaningful comparisons across different sample categories.

HEEEREEE: 8.3. DEG-Focused Visualizations (Result Exploration)

• After confirming data structure & quality, we focused on the DEGs themselves to better understand their biological patterns. Several visualization approaches were used to interpret the direction, magnitude, and distribution of gene expression changes across samples.

8.3.1 Volcano Plot of DEGs:

• The volcano plot displays each gene’s log2 fold change (x-axis) versus its statistical significance (–log10 adjusted p-value; y-axis). Upregulated genes appear on the right, downregulated genes on the left, and highly significant DEGs toward the top. This plot allows quick identification of genes with strong and consistent expression differences between lesional and non-lesional skin.

8.3.2 MA Plot:

• The MA plot displays each gene’s log2 fold change (M; Minus) (difference in expression between conditions) on the y-axis and its average normalized expression (A; Average) on the x-axis. It provides an overview of how gene expression changes relate to expression level. It confirms whether differential expression occurs across a wide expression range, meaning that both highly expressed and lowly expressed genes can show significant changes between conditions. In this dataset, DEGs were distributed symmetrically around zero, indicating balanced up- and downregulation and reliable normalization.
• In this dataset, most genes cluster closely around zero log2 fold change, indicating that expression levels are largely consistent between lesional and non-lesional samples after normalization. Only a smaller subset of genes deviates notably from zero, reflecting the differentially expressed genes (DEGs). Applying shrinkage (via apeglm) has stabilized these fold-change estimates, especially for low-count genes, resulting in a smoother and more reliable representation of overall gene expression trends.
• An ideal MA plot for well-normalized data would show the majority of genes near the log fold change of zero, with outliers representing differentially expressed genes. The plot shows a good distribution of points with no systematic bias, suggesting the data has been properly normalized. In an MA plot, the data has been properly normalized if the majority of data points are clustered symmetrically around the horizontal line at y=0 across the entire range of x-axis values. Conversely, the data has not been properly normalized if the points form a diagonal slope or a curve, indicating a systematic bias where differences between samples depend on the magnitude of the values.

8.3.3 Heatmap of Top Differentially Expressed Genes:

• A clustered heatmap of the top DEGs was generated to visualize overall expression trends. Rows represent genes, columns represent samples, and colors indicate normalized expression. The clear separation between lesional and non-lesional samples highlights the discriminatory power of these top DEGs.

8.3.4 Expression Density Plots (Top DEGs):

• Density plots of expression values for selected top DEGs illustrate how normalized expression distributions differ between conditions. These plots confirm that lesional samples show consistent directional shifts in expression, reinforcing the statistical findings.

8.3.5 Enhanced Heatmap with Annotations:

• An enhanced version of the DEG heatmap was created, combining expression data with sample annotations (such as condition or patient ID). This layout enhances interpretability by linking molecular expression trends to sample metadata, making the visualization suitable for publication and interactive web display.

8.4. Exporting DEG Tables

• To ensure reproducibility and facilitate downstream analyses, both full and filtered DEG tables were exported. The full table includes all genes tested along with their log2 fold changes, standard errors, p-values, and adjusted p-values. This comprehensive dataset is particularly suitable for analyses such as GSEA that require ranked gene lists.
• The filtered table includes only the significant DEGs that meet the specified thresholds, offering a concise dataset ideal for reporting, visualization, or functional analysis. Both outputs are formatted for reuse, ensuring consistency and transparency in downstream steps.

10. Functional Enrichment Analysis

• !!! This section is still a work in progress !!!
• This section links differential expression to biological meaning by testing whether specific pathways, processes, and networks are over-represented among DEGs.
• It includes Protein–Protein Interaction (PPI), Over-Representation Analysis (ORA) for GO terms & KEGG pathways, and Gene Set Enrichment Analysis (GSEA) for GO terms & KEGG pathways

Status:

• Live now: 10.1 PPI subnetwork of top DEGs using STRING, community detection, GO checks within modules, and centrality-based ranking.
• In progress: Over-Representation Analysis (ORA) for GO terms & KEGG pathways, and Gene Set Enrichment Analysis (GSEA) for GO terms & KEGG pathways. Figures and write-ups are being prepared and will be added here.

10.1. Protein-Protein Interaction (PPI) subnetwork of top DEGs (STRING)

• Differential Expression Analysis (DEA) identifies genes (DEGs) that show significant changes in expression between different conditions, but disease biology emerges from how these gene products work together as genes rarely act in isolation. A PPI subnetwork contextualises the DEGs within a curated map of protein interactions, so we can see which parts of the cellular machinery are co-affected, which groups of proteins move together, and which specific proteins sit at influential positions in the network. This makes PPI analysis a natural bridge from statistical signal to biological mechanism.
• This analysis serves the project objectives in three ways. First, it clusters DEGs into communities (modules) that often correspond to pathways or processes, which we test with GO enrichment to expose themes relevant to AD pathogenesis. Second, it ranks proteins by network centrality to highlight candidates that could modulate broader parts of the system when perturbed. Third, it links expression change with network position, which is useful for prioritising biomarkers and potential therapeutic targets rather than treating all DEGs equally.
• We use several centrality measures because each captures a different notion of influence. 1️⃣High degree (hubs) identifies classic hub genes, proteins that interact with many others. 2️⃣High betweenness (bottlenecks/bridges) highlights genes that sit on many shortest paths and can control information flow between modules. 3️⃣High closeness points to genes that can quickly reach the rest of the network, useful for broad signal propagation. 4️⃣High eigenvector centrality marks influential genes that are connected to other influential genes. In text, it is best to reserve the term hub for high-degree nodes, and refer to others as bottleneck genes (betweenness), high-closeness genes, or high-eigenvector-centrality genes.
• To keep interpretation coherent, we build a high-confidence STRING subnetwork from the top DEGs, detect Louvain communities, and run GO over-representation within each community. We then compute centralities on the weighted graph, using STRING interaction scores and an inverse-weight distance for shortest-path metrics. Candidate biomarkers are prioritised where ① expression change is strong, ② network centrality is high in at least one metric, and ③ the gene sits in a module enriched for processes relevant to the disease context. Candidate drug targets are those central genes that also ④ have druggable properties in the literature or databases, or that act as bridges between key modules.
• Two caveats guide interpretation. Central genes can be essential housekeeping components, so network prominence alone is not sufficient for nomination. PPI knowledgebases are incomplete and biased, so results are cross-checked with enrichment, direction of change, and the known biology of the system. With these controls, the PPI subnetwork converts a flat DEG list into structured hypotheses: modules that hint at perturbed pathways, and a short list of influential genes to carry into biomarker discovery and target assessment.
• The following subsections implement this plan step by step, from common setup and DEG mapping to network construction, module enrichment, and centrality-based ranking.

10.1.1. Common setup and configuration

• This block (Code Block 10.1.1) prepares the PPI workflow’s shared environment so everything downstream uses the same thresholds, folders, and dependencies. It loads packages for STRING queries (STRINGdb), (GO enrichment clusterProfiler, org.Hs.eg.db(human genome annotation package)), graph construction and layouts (igraph, tidygraph, ggraph), and plotting utilities (ggplot2, ggrepel, grid, scales). It then defines the analysis “knobs” such as the STRING confidence cutoff, DEG filters, and how many genes to carry forward, creates a single output directory, registers filenames for tables generated later, and sets GO enrichment options. Finally, it checks that the DESeq2 objects exist. This section doesn’t write results; it standardizes configuration so all later sections read and write consistently under fig_dir.
• Regarding the above-mentioned STRING confidence cutoff, In the raw STRING database, the combined score ranges from 0 to 1000, and common confidence thresholds are >150 (low), >400 (medium), >700 (high), and >900 (very high) (as used by STRING wrappers like rbioapi). In DEG-based PPI network analyses, many studies use ≥400 (medium confidence) to retain a broader network, whereas ≥700 (high confidence) is preferred when focusing on more reliable interactions, and ≥900 (very high confidence) is sometimes used as a strict filter for validation or visualization (rbioapi docs; STRING help).

Code Block 10.1.1. Common Setup & Configuration — Packages, parameters, paths, and input checks

				 
############################################################
## STRING PPI + Modules (Louvain) + Enrichment-driven figures
##
## Context:
## Differential expression between lesional (AL) vs non-lesional (AN) skin,
## followed here by PPI construction (STRING), module detection (Louvain),
## GO-based functional panels, and centrality exports.
##
## Requirements in your R session:
## - 'res'        : DESeq2 results (rownames = ENSEMBL IDs, possibly with version)
## - 'res_shrunk' : shrinked results with 'log2FoldChange' aligned to 'res'
############################################################


# ==========================================================
# [Section 1/9] Common setup & configuration
#    Thresholds, output folder, filenames, enrichment options, input checks
# ==========================================================


# 1.1 Load packages
suppressPackageStartupMessages({          # keep console clean when loading libraries
  library(STRINGdb)                       # STRING API for mapping + interaction retrieval
  library(clusterProfiler)                # functional enrichment (e.g., GO)
  library(org.Hs.eg.db)                   # gene ID mapping for Homo sapiens
  library(igraph)                         # network construction + centralities
  library(ggraph)                         # grammar-of-graphics plotting for graphs
  library(tidygraph)                      # tidy interface around igraph
  library(ggplot2)                        # general plotting
  library(scales)                         # palettes, transformations
  library(ggrepel)                        # non-overlapping labels
  library(grid)                           # units for padding/spacing (unit())
})


# 1.2 User knobs (analysis thresholds and plotting controls)
min_score <- 700        # STRING combined_score threshold
top_k     <- 400        # top N DEGs to start with (after sorting by |LFC|)
padj_thr  <- 0.05       # adjusted p-value threshold
lfc_thr   <- 0.5        # absolute log2 fold-change threshold


# 1.3 Output location
fig_dir <- "C:/Users/Desktop/RNA_Seq_Task1/PPI"
dir.create(fig_dir, showWarnings = FALSE, recursive = TRUE)


# 1.4 Output filenames (node table, module lists, and module summaries)
out_nodes_csv  <- file.path(fig_dir, "ppi_modules_louvain_nodes.csv")
out_lists_csv  <- file.path(fig_dir, "ppi_modules_louvain_gene_lists.csv")
out_summary_csv <- file.path(fig_dir, "ppi_module_summary_full.csv")
out_summary_preview_csv <- file.path(fig_dir, "ppi_module_summary_preview.csv")


# 1.5 Enrichment figure settings
ontol             <- "BP"  # BP, MF, CC: Biological Process (BP), Molecular Function (MF), and Cellular Component (CC) 
top_terms_show    <- 10
min_genes_module  <- 5
top_n_enriched    <- 3
n_labels_per_net  <- 20


# 1.6 Input checks
# Ensure DESeq2 result objects exist in the environment
if (!exists("res") || !exists("res_shrunk")) {
  stop("Objects 'res' and 'res_shrunk' must exist in the environment.")
}

		    

10.1.2. Select DEGs and map ENSEMBL to gene symbols

• This section (Code Block 10.1.2.) selects the gene set that seeds the PPI network (STRING). Using base R operations, it removes ENSEMBL version suffixes, collects the adjusted P values and shrunken log2 fold changes, filters by significance and effect size, ranks by absolute fold change, and trims to the top k. It then converts ENSEMBL IDs to HGNC symbols with clusterProfiler::bitr and org.Hs.eg.db (the Human Genome Annotation), enforcing a one-to-one mapping to avoid ambiguous identifiers. The output is a clean symbol-level table (df2) that preserves the DE statistics (Figure 10.1.2.1). This table is the hand-off to STRING in the next section and defines the biological scope of the network.

Code Block 10.1.2. DEG Selection & ENSEMBL → SYMBOL Mapping

# ==========================================================
# [Section 2/9] DEG selection & ENSEMBL → SYMBOL mapping
#    Filter DEGs, rank by |LFC|, keep top_k; 1:1 map to HGNC symbols
# ==========================================================
 

# 2.1 Pick DEGs and prepare IDs
# -  - - - - - - - - 
#    Pipeline:
#      a) strip ENSEMBL version
#      b) collect padj and LFC from DESeq2 results
#      c) filter by padj and |LFC|
#      d) sort by |LFC| desc and keep top_k
# -  - - - - - - - - 

# Collect all gene rownames (ENSEMBL IDs, possibly with version suffix).
ens_all  <- rownames(res)
	
# Strip Ensembl version suffix if present (e.g., ENSG000001234.12 → ENSG000001234) to match mapping keys.
ens_base <- sub("\\.\\d+$", "", ens_all)
# - - - - - - 
# The above regex  removes a trailing "." (e.g., ".12"); otherwise no change:

#   \\.   -> literal dot
#   \\d+  -> one or more digits
#   $     -> end of string
# - - - - - - 

# Build a compact table of IDs and DE statistics that will feed the PPI steps.
stats_df <- data.frame(
  ENSEMBL = ens_base,
  padj    = res[ens_all, "padj"],
  LFC     = res_shrunk[ens_all, "log2FoldChange"],
  stringsAsFactors = FALSE
)

# Filter DEGs by combining significance (padj) and effect size (|LFC|)
stats_df <- subset(stats_df, !is.na(padj) & !is.na(LFC) & padj <= padj_thr & abs(LFC) >= lfc_thr)
if (nrow(stats_df) < 10) stop("Too few DEGs after filtering; adjust thresholds.")

# Rank by absolute effect size (|LFC|) and keep top_k to control network size/complexity.
stats_df <- stats_df[order(-abs(stats_df$LFC)), ]
stats_df <- head(stats_df, top_k)


# 2.2 Map ENSEMBL → SYMBOL using org.Hs.eg.db (via clusterProfiler::bitr).
# Expect occasional one-to-many mappings; we collapse to a simple 1:1 in the next line.
map <- suppressWarnings(
  bitr(stats_df$ENSEMBL, fromType = "ENSEMBL", toType = "SYMBOL", OrgDb = org.Hs.eg.db)
)
# Enforce one SYMBOL per ENSEMBL to keep node identifiers clean and unique.
map <- map[!duplicated(map$ENSEMBL), c("ENSEMBL","SYMBOL")]

# Merge DE stats with symbols; this table (df2) is the input to STRING mapping.
df2 <- merge(stats_df, map, by = "ENSEMBL")
if (nrow(df2) < 10) stop("Too few SYMBOLs after mapping; increase top_k or check IDs.")


		    

10.1.3. Section 3: Map to STRING and fetch interactions

• This step (Code Block 10.1.3) first creates a STRING object (string_db) for the human database (version 12), setting a confidence cutoff (min_score <- 700), that instructs STRING to return only interactions at or above our chosen confidence, so that only reliable, high confidence interactions are considered. Then, using this string_db object, we map our HGNC symbols from df2 to STRING protein IDs with the function string_db$map(). Any genes (rows) that cannot be mapped are dropped thanks to removeUnmappedRows = TRUE. Because some biological identifiers can collapse onto the same protein record, the script imposes a deterministic choice among duplicates by ordering rows by adjusted p-value with mapped <- mapped[order(mapped$padj), ] and then keeping one row per protein via mapped <- mapped[!duplicated(mapped$STRING_id), ]. With this ordering, the kept row is the one with the smallest padj for that protein, which favors the most statistically compelling representative when multiple symbols would otherwise create duplicate nodes in the graph (Figure 10.1.3.1).



• After the nodes are defined, we fetch all high-confidence interactions among the mapped proteins using edges <- string_db$get_interactions(mapped$STRING_id). This produces an edge list containing protein–protein pairs including three columns from, to, and combined_score (Figure 10.1.3.2). The combined_score column is not computed in the script—it comes precomputed from the STRING database. It summarizes confidence by combining multiple evidence channels (experiments, curated databases, co-expression, text mining, genomic neighborhood, gene fusion, and co-occurrence) and, in STRING v12, is typically scaled from 0 to 1000 (higher means stronger support). Then, using subset() function, we filter this list to keep only interactions above the chosen score threshold (min_score > 700), and only those where both proteins are part of our mapped set. At this stage, the table may still contain duplicated interactions written in opposite directions (e.g., A→B and B→A), which are later removed, reducing the total number of edges from 1,552 (Figure 10.1.3.2) to 776 (Figure 10.1.3.3).



• Finally, because we are treating the network as undirected, the same interaction can appear twice in the raw table (A→B and B→A). To avoid duplicates, the code builds an undirected key by lexicographically ordering the protein IDs: key <- ifelse(edges$from < edges$to, paste(edges$from, edges$to), paste(edges$to, edges$from)), and then removes repeated pairs with edges <- edges[!duplicated(key), ]. The cleaned edges list (Figure 10.1.3.3)—now limited to high-confidence, non-self, unique undirected interactions with an attached combined_score—is ready for the next section, where it’s passed to igraph to build the weighted PPI network.

Code Block 10.1.3. STRING Mapping & Interaction Retrieval — Map symbols to STRING IDs and fetch high-confidence edges

			


# ==========================================================
# [Section 3/9] STRING mapping & interaction retrieval
#    Map symbols to STRING IDs; fetch high-confidence edges; deduplicate pairs
# ==========================================================


# 3.1 Map to STRING identifiers
# -  - - - - - - - - 
# SYMBOL -> STRING IDs and INTERACTIONS
# -  - - - - - - - - 
# Initialize STRING for human (9606), requesting only edges at or above min_score.
# This setting pre-filters returned interactions by confidence.
string_db <- STRINGdb$new(
  version = "12", species = 9606,
  score_threshold = min_score, input_directory = ""
)

# map symbols to STRING IDs; remove proteins with no mapping
mapped <- string_db$map(df2, "SYMBOL", removeUnmappedRows = TRUE)


# Set a deterministic priority for many-to-one mappings (multiple SYMBOLs → same STRING_id).
# Here we prioritise the smallest padj (most significant). This only affects which duplicate is kept.
mapped <- mapped[order(mapped$padj), ]
# Keep a single row per STRING protein to avoid duplicate nodes in the graph.
# Because of the sort above, the retained row is the one with the smallest padj for that STRING_id.
mapped <- mapped[!duplicated(mapped$STRING_id), ]

	
if (nrow(mapped) < 10) stop("Too few mapped STRING IDs; adjust thresholds.")


# 3.2 Fetch high-confidence interactions among mapped proteins and deduplicate
edges <- string_db$get_interactions(mapped$STRING_id)
# Retain only the minimal columns needed downstream.
edges <- subset(
  edges,
  combined_score >= min_score &
    from %in% mapped$STRING_id &
    to   %in% mapped$STRING_id &
    from != to,
  select = c("from","to","combined_score")
)

# Build an undirected key and drop duplicates (keep one edge per pair)
key <- ifelse(edges$from < edges$to,
              paste(edges$from, edges$to),
              paste(edges$to, edges$from))

edges <- edges[!duplicated(key), ]



		    

10.1.4. Build the PPI graph and keep meaningful components

• This step (Code Block 10.1.4) begins by creating a vertex table (verts) using the mapped results from the previous section. We select the key columns — "STRING_id", "SYMBOL", "LFC", and "padj" — to retain both the protein identifiers and their differential expression statistics. The "STRING_id" is then renamed to "name" (colnames(verts)[1] <- "name") because the igraph package expects a vertex attribute called "name" to uniquely identify each node. This ensures compatibility when constructing the network graph. The resulting table (Figure 10.1.4.1) forms the foundation for node attributes, linking STRING protein IDs with their corresponding gene symbols, log fold change (LFC), and adjusted p-values.



• Next, we build an undirected graph (g) from the edge list and vertex table (verts) using the igraph::graph_from_data_frame() function. This function converts tabular edge and vertex data into an igraph network object in R. The d argument (edges, shown in Figure 10.1.3.3) in the igraph graph_from_data_frame() function is a data frame that serves as the edge list for the graph. It must contain at least two columns, which represent the "from" and "to" nodes for each edge. Any additional columns in this data frame are automatically used as edge attributes, such as combined_score. Setting directed = FALSE specifies that the graph is undirected, reflecting that STRING protein–protein interactions (PPIs) are mutual relationships without directionality. The vertices argument is an optional data frame that defines the graph’s nodes and their attributes. Its first column (name = verts$name) specifies the vertex identifiers, while the remaining columns—such as label, lfc, and padj—are automatically added as vertex attributes.
• The graph is then simplified using the simplify() function to remove loops (edges connecting a vertex to itself) and multiple edges (two or more edges connecting the same pair of vertices). Loops and redundant edges are biologically meaningless in this context, as a protein cannot interact with itself, and duplicate edges inflate connectivity artificially. Simplifying ensures a clean, non-redundant network ready for component analysis.
• The next step (Code Block 10.1.4; 4.3), by removing components with fewer than three nodes, identifies and retains only meaningful connected components. The list object comp (Figure 10.1.4.2) is created using the function components(g), which detects all connected components in the initial protein–protein interaction (PPI) network before any filtering is applied. Each connected component represents a group of proteins that are directly or indirectly connected through interaction edges. The resulting object comp contains three key elements: membership, csize, and no. The membership vector (length = 362) assigns each vertex to one of 155 components, while csize lists the size of each component (number of vertices per group), and no indicates the total number of detected components. As shown in (Figure 10.1.4.2), most components are very small, often containing only one or two nodes, whereas the largest connected component (LCC) contains 185 proteins. This distribution reveals that the initial network, built directly from the STRING database, consists of one large cluster of highly connected proteins alongside numerous small or isolated nodes. These small disconnected nodes have limited biological value in downstream network analyses and are therefore filtered out in the next step. In this step, the total number of nodes decreases from 362 to 211 after applying keep <- which(comp$csize[comp$membership] >= 3) and g <- induced_subgraph(g, vids = keep), meaning that singletons and pairs (components with fewer than three nodes) are excluded, keeping only larger, biologically meaningful clusters for downstream analysis.

Code Block 10.1.4. PPI Graph Construction & Component Filtering — Build igraph and retain components ≥ 3 nodes

			


# ==========================================================
# [Section 4/9] Graph build & component filtering
#    Build igraph with attributes; simplify; keep components ≥ 3 nodes
# ==========================================================
 

# ==========================================================
# 4.1 BUILD IGRAPH OBJECT with attributes (nodes: proteins; edges: STRING interactions)

# Assemble a vertex (node or protein) table with the attributes we’ll carry through (labels and DE stats).
verts <- unique(mapped[, c("STRING_id","SYMBOL","LFC","padj")])
# Rename STRING_id → name to meet igraph expectations.
colnames(verts)[1] <- "name"

# Build an undirected igraph from edges and vertices; 
# Attach gene label, LFC, and padj as node attributes
g <- graph_from_data_frame(
  d = edges, directed = FALSE,
  vertices = data.frame(
    name  = verts$name,
    label = verts$SYMBOL,
    lfc   = verts$LFC,
    padj  = verts$padj
  )
)


# 4.2  Simplify graph to remove multiedges and loops (average score if needed)
g <- simplify(
  g,
  remove.multiple = TRUE,
  remove.loops = TRUE,
  edge.attr.comb = list(combined_score = "mean", "ignore")
)


# 4.3  Identify connected components and keep meaningful components (≥ 3 nodes)
comp <- components(g)
keep <- which(comp$csize[comp$membership] >= 3)
g <- induced_subgraph(g, vids = keep)
if (vcount(g) == 0) stop("No components >= 3 nodes; lower score or increase top_k.")





		    

10.1.5. Section 5: Detect basic node metrics and Louvain modules and assemble tidy tables (CSV exports)

• This section (Code Block 10.1.5) computes simple per-node metrics (degree and up/down sign from LFC) and assigns community labels. It calculates degree and encodes regulation direction (up/down) from log2FC, then runs Louvain community detection (cluster_louvain) using the STRING scores as weights to define data-driven modules. From these assignments it creates two exportable tables: a long module–gene list CSV (out_lists_csv) and a node attribute CSV (DEGs with module assignments and statistics) (out_nodes_csv) that includes STRING ID, gene symbol, module, degree, log2FC, and padj.
• In this section, the network analysis proceeds with the calculation of the unweighted degree of each node in the graph g, which at this stage contains 211 proteins after filtering components smaller than three nodes. The function degree(g) is used to count the number of edges connected to each node, giving a measure of how many direct interactions each protein has in the PPI network. The term unweighted indicates that the edge weights are not considered in this calculation; every edge contributes equally to the degree, regardless of its STRING combined score. Viewing the resulting values with View(data.frame(deg)) Figure 10.1.5.1 shows that the degrees range from 45 to 1, meaning that while some proteins act as highly connected hubs, others interact with only one or a few partners.



• Examining the structure of the igraph object g (Figure 10.1.5.2) in RStudio confirms this relationship. The object appears as a list of 211 elements, each representing a vertex (node) and the indices of other vertices to which it is directly connected. For example, as shown in Figure 10.1.5.2, the node 9606.ENSP00000357726, visible as g[[6]], lists five integers — 17, 40, 41, 153, 154 — which correspond to the vertex indices of its direct neighbours, showing that this protein is connected to five others. In contrast, g[[5]], representing 9606.ENSP00000310170, contains only one index (189), indicating a single connection. These internal lists directly correspond to the degree values produced by degree(g): the number of integers within each list equals the number of connections per node.



• The computed degree values (the variable deg on the right side of the assignment) are then added to the graph as a vertex attribute using V(g)$deg <- deg. This stores the degree of each node within the graph object, making it available for later use in plots, tables, or filtering steps. The list of vertex attributes can be viewed with vertex_attr_names(g), which now includes “deg” alongside previously assigned attributes such as “name”, “label”, “lfc”, and “padj”. A categorical attribute named sign is then created from the log2 fold-change values (V(g)$lfc) using ifelse(V(g)$lfc >= 0, "Upregulated", "Downregulated"), assigning each protein a label based on its regulation direction in lesional versus non-lesional samples. At this point, each vertex (node) in the network carries both its degree (connectivity) and regulation status (sign), linking structural and biological information within the same object.
• Next, basic statistics are collected to describe the network’s scale and connectivity. The command components(g) recomputes connected components (subgraphs) in the current graph, returning the component membership of each node ($membership) (Figure 10.1.5.3), the size of every component ($csize) (Figure 10.1.5.4), and the total component count ($no).





• The object comp2, which is generated with the function components(g), summarizes the connected component structure of the network after filtering out components with fewer than three nodes. This updated graph (comp2) reflects only the meaningful clusters retained for downstream analyses. The membership vector now has 211 entries, corresponding to the remaining vertices, and the number of components (no) has been reduced from 155 to 8, indicating that most smaller or weakly connected clusters were removed during earlier filtering, leaving only the main, biologically relevant subnetworks for analysis. The csize values indicate that the largest connected component (LCC) still contains 185 nodes, while the smaller retained subnetworks each have between 3 and 6 nodes. As illustrated in Figure 10.1.5.5, compared with comp from Section 4, this reduction in component count shows that isolated or weakly connected nodes have been removed, leaving a more coherent and biologically interpretable network. This refined structure serves as the foundation for computing node metrics and identifying functional modules in the subsequent analysis steps.



• The largest connected component size is determined with max(comp2$csize), while vcount(g) and ecount(g) report the total numbers of vertices and edges. These values are useful for figure subtitles or descriptive text accompanying visualisations of the network.
• The list object comm is created using the cluster_louvain(g, weights = E(g)$combined_score) function, which performs community detection on the filtered protein–protein interaction network. This object represents the Louvain clustering output and contains six key elements: membership, memberships, modularity, names, vcount, and algorithm. The membership vector (length = 211) assigns each protein (vertex) to one of the detected modules, while memberships is a 3×211 matrix showing the hierarchical progression of clustering across three optimization levels, where the final row corresponds to the final module assignment. The modularity vector records the modularity score at each level, with the final value (≈ 0.542) indicating a strong modular organization of the network. The names field contains the STRING protein IDs associated with each vertex, and vcount confirms the number of vertices clustered (211). Finally, algorithm specifies the method used (“multi level”), which is the Louvain algorithm implemented in igraph. This structure reveals that the Louvain method successfully partitioned the weighted network into multiple biologically coherent modules based on connection strength (STRING scores), which are later used for node-level and module-level analyses (Figure 10.1.5.6).



• Each node’s community ID is stored as a factor attribute using V(g)$module <- factor(membership(comm)), and sizes(comm) provides a summary of module sizes before any filtering.
• To prepare results for export, gene symbols are grouped by their assigned module with split(V(g)$label, V(g)$module), producing a list where each element corresponds to one module and contains its constituent genes (Figure 10.1.5.7). By using split(V(g)$label, V(g)$module), we’re grouping the gene symbols by their module ID. So, the output modules_list becomes a list of modules, where each list element corresponds to one module. Only modules with at least three genes are retained (modules_list[sapply(modules_list, length) >= 3]), and the list is sorted by decreasing size (modules_list[order(sapply(modules_list, length), decreasing = TRUE)]) for easier interpretation.



• When gene symbols are grouped by their module ID, filtered, and ordered in decreasing size, the next step prepares the results for export by creating a structured node table. A vector of retained module identifiers is first extracted using keep_mods <- names(modules_list), ensuring that only the modules meeting the size threshold are carried forward. A new data frame named vdf is then constructed, containing one row per protein and summarising key node-level information, including its STRING identifier (V(g)$name), gene symbol (V(g)$label), module assignment (V(g)$module), degree (V(g)$deg), log2 fold change (V(g)$lfc), and adjusted p-value (V(g)$padj) (Figure 10.1.5.8). To ensure that the table includes only relevant proteins, it is filtered to retain rows corresponding to the selected modules using vdf <- subset(vdf, module %in% keep_mods) (Figure 10.1.5.9), and finally sorted by module and decreasing degree so that the most connected hub proteins appear first. The resulting node table is exported as a CSV file node attribute CSV for documentation, visualisation, and downstream analysis.





• Finally, a long-format (long module–gene list CSV) mapping of modules to genes is created using do.call(rbind, lapply(...)), producing a two-column data frame, mod_map (Figure 10.1.5.10), with module IDs and gene names. This file is also exported as a CSV and is particularly useful for pathway-enrichment analysis or supplementary reporting.



• Together, these steps quantify the connectivity of each protein, associate it with its regulation direction, identify community structures within the network, and generate organised outputs that capture both the topological and biological organisation of the protein-protein interaction graph.

Code Block 10.1.5. Node Metrics & Louvain Modules — Degree/sign, community detection, and CSV exports

			
# ==========================================================
# [Section 5/9] Node metrics & Louvain modules (CSV exports)
#    Degree/sign; Louvain modules; node table and module list
# ==========================================================


# ==========================================================
# 5.1 NODE METRICS (degree, sign) + network stats for figure subtitle
# ==========================================================
# Compute unweighted degree per node (vertex); store as attribute for plotting and tables.
deg <- degree(g)

# ----------------------------------------------------------
# Creating two vertex attributes on g:
# 1.deg (numeric): degree of each node
# 2. sign (character): "Upregulated" or "Downregulated" from V(g)$lfc
# ----------------------------------------------------------


# Store the degree vector as a vertex attribute named "deg".
# After this, we can access it as V(g)$deg and use it for plotting/tables.
V(g)$deg  <- deg 


# Create a categorical vertex (node) attribute "sign" from log2 fold-change (V(g)$lfc).
# Nodes with lfc >= 0 are labeled "Upregulated", others "Downregulated".
# (Assumes V(g)$lfc was set earlier when the graph was built.)
V(g)$sign <- ifelse(V(g)$lfc >= 0, "Upregulated", "Downregulated")


# ----------------------------------------------------------
# Summary counts for subtitle/context: size of largest component, total nodes, total edges.
# ----------------------------------------------------------
comp2   <- components(g) # connected-components summary: $membership, $csize (sizes), $no (count)
LCC     <- max(comp2$csize) # size of the Largest Connected Component (number of nodes in biggest component)
nodes   <- vcount(g) # total number of nodes in the current graph
edges_n <- ecount(g) # total number of edges in the current graph


# - - - - - - - 
# comp2$membership → component ID for each node

# comp2$csize → sizes of each component

# comp2$no → total count of components
# - - - - - - -

 
# 5.2 COMMUNITY DETECTION (modules via Louvain)
# - - - - - -  - - -
#    - Assign module labels per node
#    - Build a "module -> genes" list and keep modules with >= 3 genes
#    - Prepare tidy per-node table (vdf) and save CSVs
# - - - - - -  - - -

# detect Louvain communities using STRING scores as edge weights
comm <- cluster_louvain(g, weights = E(g)$combined_score) 

# store each vertex's community ID as factor attribute "module"
V(g)$module <- factor(membership(comm))

# vector of community sizes (number of vertices in each)
sizes_comm  <- sizes(comm)


# 5.3 Exports: node table and module→gene mapping
# ----------------------------------------------------------

# group gene symbols by their module ID
# split gene symbols (V(g)$label) into a list grouped by each vertex’s module ID (V(g)$module); one element per module
modules_list <- split(V(g)$label, V(g)$module)


# keep only modules with ≥ 3 genes
modules_list <- modules_list[sapply(modules_list, length) >= 3]

# sort modules by size (descending)
modules_list <- modules_list[order(sapply(modules_list, length), decreasing = TRUE)]

# vector of retained module IDs
keep_mods <- names(modules_list)


# ----------------------------------------------------------
# Build a tidy node table (one row per gene) for export/inspection and downstream linking.
# ----------------------------------------------------------

vdf <- data.frame(
  STRING_id = V(g)$name,
  gene      = V(g)$label,
  module    = V(g)$module,
  deg       = V(g)$deg,
  lfc       = V(g)$lfc,
  padj      = V(g)$padj,
  stringsAsFactors = FALSE
)


# Keep only nodes in retained modules; 
vdf <- subset(vdf, module %in% keep_mods)


# sort by module then by degree to surface hubs first. 
## Save node table and module gene lists
write.csv(vdf[order(as.numeric(vdf$module), -vdf$deg), ], out_nodes_csv, row.names = FALSE)


# Long-format module→gene map (useful for enrichment, reporting, and reproducibility).
# Save module -> gene mapping (long format)
mod_map <- do.call(rbind, lapply(names(modules_list), function(m) {
  data.frame(module = m, gene = modules_list[[m]], stringsAsFactors = FALSE)
}))

write.csv(mod_map, out_lists_csv, row.names = FALSE)


		    

10.1.6. Section 6: Module summaries and size overview

• The section begins by normalizing data types so later tabulations behave predictably. Converting module IDs and gene symbols to plain character with mod_map$module <- as.character(mod_map$module) and mod_map$gene <- as.character(mod_map$gene) prevents factor-level issues when splitting and pasting strings. We then count how many genes fall into each module using mod_counts <- sort(table(mod_map$module), decreasing = TRUE), which yields a size-ordered list of module counts so the largest modules appear first throughout the summaries (Figure 10.1.6.1).



• Next, a comprehensive, module-level summary is created. The call split(mod_map$gene, mod_map$module) groups genes by module; Within each group, we use lapply(..., function(x) { ux <- unique(x); paste(ux, collapse = ", ") }) to create a single comma-separated string of unique gene names. The lapply() function applies the small function function(x) to every element of the list created by split(). Inside it, unique(x) removes any repeated gene names, and paste(..., collapse = ", ") joins them together into one line separated by commas (Figure 10.1.6.2). The list examples_full is reordered to follow the same module order as mod_counts, which was already sorted from the largest to smallest module (Figure 10.1.6.3).





• This ensures that when we later combine these results into module_summary_full (Figure 10.1.6.4), the rows appear in descending order of module size, with the largest modules listed first. The result is written to disk with write.csv(module_summary_full, out_summary_csv, row.names = FALSE) and exposed on our site as a browsable “full” CSV; ⬇️ Download CSV (PPI module full summary).



• For fast scanning, a compact preview table is also generated. Using the same split-then-collapse pattern also generated. Using the same split-then-collapse pattern, examples_preview (Figure 10.1.6.5) keeps only the first 10 unique symbols per module: ux[seq_len(min(10, length(ux)))]. These snippets become the example_genes column in module_summary_preview, paired with the ordered module and n_genes. It is saved via write.csv(module_summary_preview, out_summary_preview_csv, row.names = FALSE) and published as a browsable “preview” CSV; ⬇️ Download CSV (PPI module preview summary).



• Finally, the section renders a quick visual overview of module sizes. A bar chart is built from module_summary_preview using ggplot2: ggplot(..., aes(x = reorder(module, -n_genes), y = n_genes, fill = factor(module))) draws one column per module, labels the bars with counts via geom_text, and applies a readable theme and viridis palette. The figure is saved as a high-resolution PNG with ggsave(file.path(fig_dir, "ppi_module_sizes_barplot_colored.png"), p_sizes, width = 8, height = 4.5, dpi = 320, bg = "white"). This produces a clear size overview suitable for reports (Figure 10.1.6.6).

Code Block 10.1.6. Module Summaries & Overview Figure — Full/preview CSVs and module size bar plot

# ==========================================================
# [Section 6/9] Module summaries & overview figure
#    - Full/preview CSVs; module size barplot (PNG)
# ==========================================================


# 6.1 MODULE SUMMARIES
# - - - -  - - -
#     - Full list (all genes in a single cell)
#     - Preview (first 10 genes) + size bar plot
# - - - -  - - -
# Make sure types are plain character before tabulating to avoid factor pitfalls.
mod_map$module <- as.character(mod_map$module) 
mod_map$gene   <- as.character(mod_map$gene)

# Count genes per module and order modules by size (largest first).
mod_counts <- sort(table(mod_map$module), decreasing = TRUE)

# Build a “full” summary: one row per module with all its genes collapsed to a single string.
examples_full <- lapply(split(mod_map$gene, mod_map$module), function(x) {
  ux <- unique(x); paste(ux, collapse = ", ")
})

# Reorder the list so its module names follow the same order as mod_counts 
# (largest modules first, ensuring consistency in later tables and plots)	  
examples_full <- examples_full[names(mod_counts)]


module_summary_full <- data.frame(
  module        = names(mod_counts),
  n_genes       = as.integer(mod_counts),
  genes_all     = unname(unlist(examples_full)),
  stringsAsFactors = FALSE
)

write.csv(module_summary_full, out_summary_csv, row.names = FALSE)

# Build a short “preview” summary: first 8 genes per module for quick scanning.
examples_preview <- lapply(split(mod_map$gene, mod_map$module), function(x) {
  ux <- unique(x); paste(ux[seq_len(min(10, length(ux)))], collapse = ", ")
})
examples_preview <- examples_preview[names(mod_counts)]
module_summary_preview <- data.frame(
  module        = names(mod_counts),
  n_genes       = as.integer(mod_counts),
  example_genes = unname(unlist(examples_preview)),
  stringsAsFactors = FALSE
)

write.csv(module_summary_preview, out_summary_preview_csv, row.names = FALSE)



# 6.2 Quick bar plot of module sizes (for overview)
library(ggplot2)

p_sizes <- ggplot(module_summary_preview,
                  aes(x = reorder(module, -n_genes), y = n_genes, fill = factor(module))) +
  geom_col(width = 0.8, color = "grey60") +
  geom_text(aes(label = n_genes), vjust = -0.35, size = 3) +
  scale_fill_viridis_d(option = "C", end = 0.95) +  # rich, non-white colors
  labs(x = "Module ID", y = "Number of genes", title = "PPI module sizes") +
  theme_minimal(base_size = 12) +
  theme(legend.position = "none",
        panel.grid.major.x = element_blank(),
        plot.title = element_text(face = "bold", size = 16),
        axis.title = element_text(face = "bold")) +
  expand_limits(y = max(module_summary_preview$n_genes) * 1.12)

# saved as PNG for documentation.
ggsave(file.path(fig_dir, "ppi_module_sizes_barplot_colored.png"),
       p_sizes, width = 8, height = 4.5, dpi = 320, bg = "white")





		    

10.1.7. Section 7: Full-network visualizations

• To contextualize the network at a glance, this section renders two complementary overviews using tidygraph::as_tbl_graph and ggraph with a Fruchterman–Reingold layout. One plot colors nodes by regulation direction to highlight expression patterns across the network; the other colors by module to display community structure. In both, node size reflects degree, a fixed number of high-degree hubs are labeled with ggraph::geom_node_text, and output images are saved as high-resolution PNGs. These figures help readers see overall architecture—dense regions, cross-module links, and prominent hubs—before drilling into functional details.
• Regarding the Outputs, the first visualization (Figure 10.1.7.1) shows the protein–protein interaction (PPI) subnetwork of the top differentially expressed genes (DEGs), with nodes sized according to their connectivity (degree) and colored by regulation status. Red nodes indicate upregulated genes, while blue nodes indicate downregulated genes. The labels highlight the top 25 hub genes with the highest connectivity, providing insight into potential key drivers within the network. This plot emphasizes how highly connected upregulated genes dominate the central network structure, while downregulated genes appear less frequent and more peripheral.
• The second visualization (Figure 10.1.7.2) displays the same network, but this time nodes are colored by Louvain-detected modules, which represent functionally or topologically distinct gene communities. Each color corresponds to a different module, with the legend identifying all 16 modules that meet the minimum size threshold (≥ 3 genes). This view highlights modular structure within the PPI network, making it easier to identify groups of genes that may be co-regulated or share related biological functions.
• Regarding the code, the code begins by converting the network object (`g`) into a tidygraph format (`as_tbl_graph`) suitable for visualization with ggraph. It defines the number of hub labels (`n_labels_target = 25`) and selects the top-degree genes to reduce clutter. For the first plot, the layout is fixed with a random seed for reproducibility, and edges are drawn in faint black lines. Nodes are sized by degree and colored according to their regulation sign (`Upregulated` or `Downregulated`), with hub labels placed using repulsion (repel = TRUE in geom_node_text() to space them apart) to avoid overlap. Titles and subtitles provide key statistics such as the number of nodes, edges, and largest connected component (LCC).
• For the second plot, the focus shifts to community structure. The Louvain module assignments are ensured to be factors for consistent categorical coloring. A custom color palette is generated using RColorBrewer’s `"Set2"` scheme, extended with `colorRampPalette` to produce distinct colors for all modules. Each module is then assigned a unique color, and nodes are drawn with size proportional to degree and colored by module. Hub labels are again selectively displayed to maintain readability. The result is a pair of complementary visualizations: one showing regulatory status across the full PPI network, and the other revealing its modular community structure.

Code Block 10.1.7. Full-Network Visualizations — FR layout colored by LFC sign and by module

			
# ==========================================================
# [Section 7/9] Full-network visualizations (PNGs)
#    - FR layout, (a) colored by LFC sign, (b) colored by module
# ==========================================================

# 7.1 Plot: full network colored by LFC sign
# Convert to a tidygraph object for ggraph plotting.
tg <- as_tbl_graph(g) 
n_labels_target <- 25
# Choose which nodes to label in the plot: top-degree hubs to reduce text clutter.
label_nodes <- names(sort(deg, decreasing = TRUE))[seq_len(min(n_labels_target, length(deg)))]

edge_col <- "black"
# Main title and subtitle showing key network stats for context.
title_main <- sprintf("PPI subnetwork of top DEGs (STRING score >= %d)", min_score)
subtitle_main <- sprintf(
  "Full network (components >= 3). LCC = %d | nodes = %d | edges = %d. labels = top %d hubs.",
  LCC, nodes, edges_n, n_labels_target
)

# Fix the force-directed layout seed for reproducible node placement across runs.
set.seed(42)
p_lfc <- ggraph(tg, layout = "fr") +
  geom_edge_link0(edge_width = 0.10, edge_alpha = 0.18, lineend = "round", colour = edge_col) +
  geom_node_point(aes(size = deg, color = sign)) +
  ggraph::geom_node_text(
    aes(label = ifelse(name %in% label_nodes, label, "")),
    repel = TRUE, max.overlaps = Inf, force = 2,
    box.padding = unit(0.35, "lines"), point.padding = unit(0.30, "lines"),
    segment.size = 0.12, size = 3.0
  ) +
  coord_cartesian(clip = "off") +
  scale_size(range = c(2.2, 6.0), name = "Degree") +
  scale_color_manual(values = c("Upregulated" = "firebrick", "Downregulated" = "steelblue"), name = "LFC sign") +
  ggtitle(title_main, subtitle = subtitle_main) +
  theme_void() +
  theme(plot.title = element_text(hjust = 0.5, face = "bold", size = 16),
        plot.subtitle = element_text(hjust = 0.5, size = 10),
        legend.position = "right",
        plot.margin = margin(20, 60, 20, 20))

ggsave(file.path(fig_dir, "7.1.ppi_full_network_lfc.png"), p_lfc, width = 8, height = 6, dpi = 300, bg = "white")


# 7.2 Plot: full network colored by module (categorical palette per Louvain community)
# Ensure module is a factor (so colors are discrete and stable).
V(g)$module <- as.factor(V(g)$module)
module_levels <- levels(V(g)$module)


# Distinct categorical palette for modules; one color per Louvain community.
# Here we use RColorBrewer 'Set2' and extend it with colorRampPalette 
# to provide clearer separation when many modules are present.
n_mod <- length(module_levels)
pal_mod <- colorRampPalette(RColorBrewer::brewer.pal(8, "Set2"))(n_mod)
names(pal_mod) <- module_levels

set.seed(42)
p_mod <- ggraph(tg, layout = "fr") +
  geom_edge_link0(edge_width = 0.10, edge_alpha = 0.18, lineend = "round", colour = "black") +
  geom_node_point(aes(size = deg, color = module)) +
  ggraph::geom_node_text(
    aes(label = ifelse(name %in% label_nodes, label, "")),
    repel = TRUE, max.overlaps = Inf, force = 2,
    box.padding = unit(0.35, "lines"), point.padding = unit(0.30, "lines"),
    segment.size = 0.12, size = 3.0
  ) +
  scale_size(range = c(2.2, 6.0), name = "Degree") +
  scale_color_manual(values = pal_mod, breaks = module_levels, name = "Module") +
  ggtitle(sprintf("PPI modules (Louvain). %d modules (>= 3 genes)", length(keep_mods))) +
  theme_void() +
  theme(plot.title = element_text(hjust = 0.5, face = "bold", size = 16),
        legend.position = "right",
        plot.margin = margin(20, 60, 20, 20))

ggsave(file.path(fig_dir, "7.2.ppi_full_network_modules.png"), p_mod, width = 8, height = 6, dpi = 300, bg = "white")

		    

10.1.8. Section 8: GO enrichment panels per module

• This part examines biological coherence within modules with GO over-representation analysis via clusterProfiler::enrichGO for the chosen ontology. It ranks modules by their best adjusted P value, writes a CSV overview of enrichment scores, and, for the top modules, saves two visuals per module: a focused network subgraph colored by log2FC and a dotplot summarizing the most significant GO terms. The code also uses the patchwork package to combine the network plot and the GO enrichment dot plot into a single side-by-side panel for each top module. Full enrichment tables per module are also exported to CSV. These outputs connect structure to function by showing what processes dominate each module and which genes drive those terms—useful for interpreting pathways relevant to atopic dermatitis.
• What this section does & what it produces (in the order the code runs). In 8.1, the code performs GO enrichment analysis for each network module. At the start of 8.1, a helper function, make_enrich_plot(), is defined. This function takes an enrichment result from enrichGO and turns it into a compact GO dot plot. It computes −log10(adj.P), parses the GeneRatio into numeric form, orders terms by adjusted P-value (with GeneRatio as a tie-breaker), and keeps only the top N (top_n = 10). It then draws a dot plot where the x-axis is −log10(adj.P), the y-axis is the GO terms, and dot size corresponds to the hit count. This helper will later be used in 8.2 to visualize selected modules. In 8.1.1, the code performs GO enrichment analysis for each network module, but only keeps modules with at least min_genes_module ≥ 5 genes, since very small modules are not reliable for ORA. With the current settings (min_genes_module = 5, ontology ontol = "BP"), it applies clusterProfiler::enrichGO() to each remaining module. From each enrichment result, it extracts the best term—the row with the smallest adjusted p-value (using GeneRatio as a tie-breaker)—and calculates a module enrichment score as −log10(best adjusted p-value). In 8.1.2, these one-line summaries are merged into a single module enrichment overview (CSV). This file includes the columns module, n_genes, best_term, best_padj, n_sig_terms (number of significant GO terms with adjusted p ≤ 0.05), and score. Each of these columns provides complementary insights: best_term is the GO term with the smallest adjusted p-value in that module, making it the most statistically significant biological process identified, best_padj is the smallest adjusted p-value (multiple-testing corrected p-value) among all enriched GO terms in a module, and it serves as the module’s “strength of enrichment” marker. Additionally, score is defined as −log10 of the adjusted p-value, which transforms the adjusted p-value into an intuitive scale where higher values indicate stronger enrichment and easier comparison across modules. n_sig_terms is the number of GO terms in the module with adjusted p ≤ 0.05. It reflects how many different biological processes are significantly enriched in that module, giving a sense of whether enrichment is narrow (few terms) or broad (many terms). Finally, n_genes is the number of genes contained in each module, reflecting the module’s size and the potential statistical power of its enrichment analysis. In 8.1.3, the code generates two overview figures from this table: (1) a bar plot of module enrichment scores (Figure 10.1.8.1) where higher bars indicate stronger enrichment; (2) a plot of the best GO term per module (Figure 10.1.8.2), which positions each module by its score and labels it with its top BP term. These visualizations provide an intuitive overview of which modules are statistically strongest and what biological processes they represent. In 8.1.4, the code selects the top N enriched modules for detailed visualization using mods_to_plot <- head(enrich_summary$module, top_n_enriched). With the setting top_n_enriched = 3, the 3 highest-scoring modules are selected for detailed outputs. Finally, in 8.2, the code produces per-module outputs for the top three enriched modules (11, 1, and 7). For every selected module, it writes an enrichment CSV file (Module 11, Module 1, and Module 7) containing all enriched GO terms, with standard clusterProfiler fields: ID (GO identifier), Description (term name), GeneRatio (proportion of module genes involved), BgRatio (proportion in the background set), p.adjust (adjusted p-value, the Benjamini-Hochberg (BH) correction), geneID (list of associated genes), and Count (number of genes contributing). In addition, a combined panel image is generated for each module that shows side-by-side: (1) the sub-network of the module (FR layout; node size = within-module degree; node color = log2FC), and (2) the GO dot plot of the top enriched terms (x = −log10 adj.P, dot size = hit count). These combined figures integrate structural network context with functional biological interpretation. For the top three enriched modules, see the figures for Module 11 (Figure 10.1.8.3), Module 1 (Figure 10.1.8.4), and Module 7 (Figure 10.1.8.5) .
• How it works, technically — The implementation uses igraph, tidygraph, and ggraph for network handling and visualization, and clusterProfiler plus org.Hs.eg.db for GO enrichment analysis. In 8.1.1, the code builds modules_list_f by filtering modules_list to retain only modules that meet the size threshold. It then loops over these modules, applies enrichGO(), and coerces each result into a data frame. The rows are ordered by adjusted p-value (with GeneRatioNum as a tie-breaker if needed), and the “best term” is extracted. A module’s enrichment score is then calculated as −log10(best adjusted p-value). In 8.1.2, the summaries are combined into enrich_summary, sorted by score (with module size as a secondary factor), and written as module enrichment overview.csv. In 8.1.3, the code builds two ggplot2 outputs: a geom_col() bar chart of enrichment scores and a geom_point() + geom_text() visualization of best terms. In 8.1.4, the top N modules are selected via head(). Then in 8.2, the code induces the subgraph for each selected module (induced_subgraph()), computes within-module degree, and labels hub nodes to avoid clutter. The module network is drawn with a force-directed FR layout, using faint edges, degree-scaled nodes, and a diverging log2FC color gradient. The companion GO plot is generated by make_enrich_plot(), which ranks enriched terms by adjusted p-value, keeps the top terms, and plots them against −log10(adj.P) with dot size = Count. If the patchwork package is available, it combines both plots into a single combo panel PNG. Overall, this section integrates a global overview of module enrichment with detailed structural and functional panels for the top modules, giving both high-level prioritization and deep biological insight.

Code Block 10.1.8. GO Enrichment per Module — Rank modules and plot network + GO dotplot panels

			
# enrich_summary 
# ==========================================================
# [Section 8/9] Enrichment-driven module panels (PNGs + CSVs)
#    - Rank modules by best adj.P (GO); per-module network + GO dotplot
# ==========================================================
 
 # 8.1 Rank modules by enrichment and write overview CSV
# -----------------------------------------------------
# Helper function: Create GO dotplot for a module
make_enrich_plot <- function(eg, mod_id, top_n = 10) {
  # Convert enrichGO result to a data frame (returns NULL if error)
  df <- tryCatch(as.data.frame(eg), error = function(e) NULL)
  
  # If no enrichment results, return a placeholder empty plot
  if (is.null(df) || nrow(df) == 0) {
    return(ggplot() + theme_void() +
             ggtitle(paste("No significant GO terms for module", mod_id)))
  }
  
  # Add –log10 adjusted p-values for ranking/plotting
  df$neglog10padj <- -log10(df$p.adjust)
  
  # Convert GeneRatio from "5/20" to numeric (0.25)
  df$GeneRatioNum <- vapply(strsplit(df$GeneRatio, "/"),
                            function(x) as.numeric(x[1]) / as.numeric(x[2]),
                            numeric(1))
  
  # Rank terms by adjusted p-value (smallest first), then by gene ratio
  df <- df[order(df$p.adjust, -df$GeneRatioNum), ]
  
  # Keep only the top N terms (default 10)
  df <- head(df, top_n)

  # Create dot plot:
  #   x = –log10(adj P), y = GO terms, dot size = gene count
  ggplot(df, aes(x = neglog10padj,
                 y = reorder(Description, neglog10padj),
                 size = Count)) +
    geom_point() +
    labs(x = "-log10(adj P)", y = NULL,
         title = paste0("GO ", ontol, " enrichment • module ", mod_id)) +
    theme_minimal(base_size = 11)
}
	    

# 8.1.1 Run enrichment per module and extract one-line summaries:
#   - Start with modules_list, which contains all detected modules (module ID → vector of genes).
#   - Very small modules (fewer than min_genes_module genes) usually cannot produce reliable GO enrichment.
#   - Therefore, we filter modules_list to keep only those with at least min_genes_module genes (e.g., ≥5 genes).
#   - The result, modules_list_f, is the filtered list of modules used for all downstream enrichment analysis.
modules_list_f <- modules_list[sapply(modules_list, length) >= min_genes_module]

# For each module, run enrichGO and extract a one-line summary (best term and its adjusted p-value). 
# tryCatch ensures a single module failure doesn't stop the whole panel generation.
enrich_per_module <- lapply(names(modules_list_f), function(m) {
  gsym <- unique(modules_list_f[[m]])
  eg <- tryCatch(
    enrichGO(gene = gsym, OrgDb = org.Hs.eg.db, keyType = "SYMBOL",
             ont = ontol, pAdjustMethod = "BH", pvalueCutoff = 0.05, readable = TRUE),
    error = function(e) NULL
  )
  df <- tryCatch(as.data.frame(eg), error = function(e) NULL)
  if (!is.null(df) && nrow(df) > 0) {
    best <- df[order(df$p.adjust, df$pvalue), ][1, ]
    data.frame(module = m, n_genes = length(gsym),
               best_term = best$Description, best_padj = best$p.adjust,
               n_sig_terms = sum(df$p.adjust <= 0.05),
               score = -log10(best$p.adjust), stringsAsFactors = FALSE)
  } else {
    data.frame(module = m, n_genes = length(gsym),
               best_term = NA_character_, best_padj = NA_real_,
               n_sig_terms = 0, score = 0, stringsAsFactors = FALSE)
  }
})

# 8.1.2 Generate enrichment summary table (sortable CSV):
#  - Overview of enrichment across modules (sortable CSV used to pick top modules for plotting).
enrich_summary <- do.call(rbind, enrich_per_module)
enrich_summary <- enrich_summary[order(-enrich_summary$score, -enrich_summary$n_genes), ]
write.csv(enrich_summary, file.path(fig_dir, "module_enrichment_overview.csv"), row.names = FALSE)

 
# 8.1.3 Create overview plots from enrich_summary
library(ggplot2)
library(stringr)

# keep modules that were tested and have a score (score = -log10(best_padj))
enrich_plot_df <- subset(enrich_summary, is.finite(score) & score > 0)

# order modules by score (desc), then by size to keep ties stable
enrich_plot_df$module <- factor(
  enrich_plot_df$module,
  levels = enrich_plot_df$module[order(-enrich_plot_df$score, -enrich_plot_df$n_genes)]
)

# 8.1.3 (A) Bar chart of enrichment scores per module
p_scores <- ggplot(enrich_plot_df,
                   aes(x = module, y = score, fill = module)) +
  geom_col(width = 0.8, color = "grey60") +
  geom_text(aes(label = round(score, 1)), vjust = -0.35, size = 3) +
  scale_fill_viridis_d(option = "C", end = 0.95, guide = "none") +
  labs(x = "Module ID", y = "-log10(best adj P)",
       title = "Module enrichment scores") +
  theme_minimal(base_size = 12) +
  theme(legend.position = "none",
        panel.grid.major.x = element_blank(),
        plot.title = element_text(face = "bold", size = 16),
        axis.title = element_text(face = "bold")) +
  expand_limits(y = max(enrich_plot_df$score, na.rm = TRUE) * 1.12)

ggsave(file.path(fig_dir, "module_enrichment_scores_barplot.png"),
       p_scores, width = 8, height = 4.5, dpi = 320, bg = "white")


# 8.1.3 (B) Best term per module (dot + label), positioned by score
# Wrap long term names so labels don’t run off the page
enrich_plot_df$best_term_wrapped <- stringr::str_wrap(enrich_plot_df$best_term, width = 40)

p_best <- ggplot(enrich_plot_df,
                 aes(x = score, y = reorder(module, score))) +
  geom_point(size = 2.8) +
  geom_text(aes(label = best_term_wrapped),
            hjust = 0, nudge_x = 0.25, size = 3) +
  labs(x = "-log10(best adj P)", y = "Module ID",
       title = "Top GO BP term per module") +
  theme_minimal(base_size = 12) +
  theme(plot.title = element_text(face = "bold", size = 16)) +
  # leave space on the right for labels
  coord_cartesian(xlim = c(0, max(enrich_plot_df$score, na.rm = TRUE) * 1.35))

ggsave(file.path(fig_dir, "module_enrichment_best_terms.png"),
       p_best, width = 10, height = 5.5, dpi = 320, bg = "white")
# ============================================================================




# 8.1.4 Pick the top N modules by enrichment score, -log10(best adjusted p-value),  for detailed panels
mods_to_plot <- head(enrich_summary$module, top_n_enriched)



# 8.2 Per-module network + GO dotplot (+ optional combined)
# Loop over selected modules: save network + GO dotplot (+ combined)
for (m in mods_to_plot) {
  # 8.2.1 Write full enrichment table for this module
  gsym <- unique(modules_list_f[[m]])
  eg <- tryCatch(
    enrichGO(gene = gsym, OrgDb = org.Hs.eg.db, keyType = "SYMBOL",
             ont = ontol, pAdjustMethod = "BH", pvalueCutoff = 0.05, readable = TRUE),
    error = function(e) NULL
  )
  
  # Inside the loop: write full enrichment table for transparency and later interpretation.
  df_full <- tryCatch(as.data.frame(eg), error = function(e) NULL)
  if (!is.null(df_full) && nrow(df_full) > 0) {
    write.csv(df_full, file.path(fig_dir, paste0("module_", m, "_GO_", ontol, ".csv")), row.names = FALSE)
  }

  # 8.2.2 Build and plot network for this module, sized by within-module degree.
  g_sub <- induced_subgraph(g, vids = V(g)[module == m])
  deg_sub <- degree(g_sub)
  V(g_sub)$deg_sub <- deg_sub

  # Label only the highest-degree nodes within this module to keep the network legible.
  label_ids <- names(sort(deg_sub, decreasing = TRUE))[seq_len(min(n_labels_per_net, length(deg_sub)))]
  label_flag <- ifelse(V(g_sub)$name %in% label_ids, V(g_sub)$label, "")

  # Network plot for this module (FR layout, sized by degree, colored by LFC)
  tg_sub <- as_tbl_graph(g_sub)
  set.seed(42)
  p_net <- ggraph(tg_sub, layout = "fr") +
    geom_edge_link0(edge_width = 0.25, edge_alpha = 0.25, colour = "gray40") +
    geom_node_point(aes(size = deg_sub, color = lfc)) +
    ggraph::geom_node_text(aes(label = label_flag), repel = TRUE, size = 3) +
    scale_size(range = c(2.5, 7), name = "Degree") +
    scale_color_gradient2(low = "steelblue", mid = "grey85", high = "firebrick",
                          midpoint = 0, name = "log2FC") +
    ggtitle(paste0("Module ", m, " network (n=", vcount(g_sub), ")")) +
    theme_void() +
  theme(
    plot.title       = element_text(hjust = 0.5, face = "bold", size = 14),
    plot.background  = element_rect(fill = "white", colour = NA),
    panel.background = element_rect(fill = "white", colour = NA)
  )
  

  # 8.2.3 GO dotplot for this module (top terms)
  p_enr <- make_enrich_plot(eg, m, top_n = top_terms_show) +
    theme(
    plot.background  = element_rect(fill = "white", colour = NA),
    panel.background = element_rect(fill = "white", colour = NA)
  )

  # 8.2.4 Save outputs:
  # Save per-module network and GO dotplot figures.
  ggsave(file.path(fig_dir, paste0("module_", m, "_network.png")), p_net, width = 7, height = 5, dpi = 300, bg = "white")
  ggsave(file.path(fig_dir, paste0("module_", m, "_GO_", ontol, "_dotplot.png")), p_enr, width = 6, height = 5, dpi = 300, bg = "white")

    # 8.2.5 combined panel (network + enrichment)
  # If patchwork is available, also save a side-by-side composite (network + enrichment) for quick review.
  if (requireNamespace("patchwork", quietly = TRUE)) {
    combined <- p_net + p_enr + patchwork::plot_layout(widths = c(1.6, 1))
    ggsave(file.path(fig_dir, paste0("module_", m, "_combo.png")), combined, width = 12, height = 5, dpi = 300, bg = "white")
  }
}


 


		    

10.1.9. Section 9: Centrality metrics and module-level summaries

• Finally, the workflow quantifies node importance using standard centralities from igraph. It defines a shortest-path distance as the inverse of STRING score so stronger interactions are “closer,” then computes a panel of centrality metrics on the protein–protein interaction (PPI) network: degree (and normalised degree), weighted strength, betweenness and closeness on the distance metric, and eigenvector centrality on the raw scores. These measures quantify how important individual genes (nodes) are in the network and how modules differ in their connectivity. A high degree means “many interaction partners”, betweenness means “acts as a bridge”, eigenvector means “influential through influential neighbours.” These map intuitively to biological concepts like hubs, bottlenecks, and master regulators. However, closeness just means “on average near others”, which is harder to tie directly to a functional role like control or influence, so less direct biological interpretation.
• The results are written to two output tables: A node-level CSV records all metrics per gene, and a module-level summary aggregates means and medians and identifies the top-degree hub in each module. These tables support ranking candidate biomarkers or targets, comparing modules by their internal connectivity, and cross-referencing highly central genes with enriched processes discovered earlier.
• The first file, ppi_node_centrality.csv, provides a per-gene table that links network properties with differential expression results. This table allows direct ranking of genes within or across modules to prioritize candidates for follow-up analysis (e.g., identifying hub genes that are also strongly upregulated). Each entry begins with the STRING identifier (in the format 9606.ENSP…, denoting a human Ensembl protein ID) alongside the corresponding gene symbol and its assigned module. The subsequent columns capture a panel of centrality measures: Degree centrality (the count of direct connections), Normalised degree (degree scaled to network size for comparability), Strength (a weighted degree summing STRING confidence scores), Betweenness centrality (how often a gene sits on shortest paths, reflecting its bridging role), Closeness centrality (how close a gene is, on average, to all others), and Eigenvector centrality (influence through connections to other influential nodes). Finally, differential expression values, the log fold change (lfc) and the adjusted p-value (padj), are appended. Together, these measures allow each gene to be ranked not only by how central it is within the protein–protein interaction network but also by how strongly and significantly its expression is altered, offering a combined view of structural importance and transcriptional relevance.
• Degree values span from 48 down to 1, normalised degree from about 0.21 to 0.004, and weighted strength from roughly 36,752 to 931. Because degree, normalised degree, and strength are mathematically linked, they follow the same ordering: nodes with more connections also accumulate higher total STRING confidence scores, and normalised degree is simply degree scaled to the size of the network. This pattern shows that a handful of genes function as highly connected hubs, while the majority maintain only a few interactions, reflecting a hub-and-spoke organisation (a “spoke” is a protein that is connected to a highly-connected “hub” protein but has few other connections of its own). Betweenness ranges from about 0.24 down to 0, with only a subset of genes reaching notable values. This indicates that while many nodes connect locally, relatively few serve as key brokers that lie on shortest paths between otherwise separate parts of the network. The presence of such high-betweenness nodes suggests points of control or vulnerability, where removing or perturbing them could disrupt communication across modules. Genes with the highest betweenness scores (up to 0.24 in this case) are the hubs or bottlenecks. They control the flow of information because they lie on a high proportion of the shortest paths between other pairs of genes. Closeness spans from roughly 996 at the high end to 110 at the low end. Higher values mark genes that are centrally embedded and can quickly reach most others, while lower values reflect more peripheral positions. This distribution shows a gradient in accessibility: some genes are located at the heart of the interaction map, while others sit at the edges, requiring longer chains of intermediaries to connect. Eigenvector values run from 1 down to 0, highlighting a sharp hierarchy. A small group of genes emerge as maximally influential, connected not only to many others but specifically to other well-connected hubs. Most genes, by contrast, cluster near zero, indicating limited influence. This separation underlines the dominance of a tightly interconnected core that exerts disproportionate structural control within the network.
• To highlight genes that are not only central in the PPI network but also strongly altered in expression, a combined filter was applied: all genes were first sorted in descending order by Degree, then genes in the top 10% for Degree centrality (highly connected hubs), Betweenness centrality (key bridges between modules), and Eigenvector centrality (influential nodes linked to other influential nodes) were selected. From this overlap, we further required a log fold change greater than 1 & adjusted p-value less than 0.0001 (Figure 10.1.9.1). Closeness centrality, however, was not included as a filter because it can give unstable values in disconnected networks and is harder to link directly to functional influence than degree, betweenness, or eigenvector centrality.
• This stringent selection highlights a focused set of genes such as IL1B, CCR7, PTPRC, FCGR3A, CD28, CCL2, CTLA4, CD2, CXCL10, GZMB, CCL5, CD69, and CCR2, which consistently rank as hubs, bridges, and core influencers while also showing strong upregulation with statistically robust adjusted p-values. Their biological roles further underline their relevance to the disease biology and collectively, these candidates map tightly onto atopic dermatitis biology:
• IL1B encodes interleukin-1β, a potent pro-inflammatory cytokine central to the initiation and maintenance of skin inflammation; CCR7 and CCR2, together with their ligands CCL2 and CCL5, coordinate leukocyte trafficking—specifically guiding the migration of T cells, dendritic cells, and monocytes into inflamed skin—thereby sustaining the immune infiltrate in lesional tissue; CXCL10 attracts activated T cells; PTPRC (CD45), CD2, CD28, and CD69 mark and modulate T-cell activation states, while CTLA4 provides inhibitory counter-signalling; GZMB denotes cytotoxic effector activity; and FCGR3A (CD16) links antibody-dependent pathways to innate immune responses. Within the context of this analysis—contrasting lesional (AL) and non-lesional (AN) skin before treatment in atopic dermatitis, these genes stand out as network-integrated biomarkers—they sit at hub/bridge/core positions in the PPI network and are robustly up-regulated—supporting their relevance for AD diagnosis and monitoring, and flagging them as promising entry points for therapeutic exploration.

• The second file, ppi_module_centrality_summary.csv, shifts the perspective from individual genes to entire network modules (Figure 10.1.9.2). Each row represents a module and reports the average and median values of all centrality measures across its member genes, together with the number of genes it contains and the single top hub gene (the individual gene within the module that has the highest unweighted degree). This summary captures not only how strongly individual genes behave as hubs, but also which clusters of genes form structurally dominant “neighborhoods” in the protein–protein interaction network, as illustrated in Figure 10.1.9.2.
• When modules are sorted in descending order by mean degree, the upper quartile (top 25%) serves as a benchmark to visualise where the highest values cluster. In Excel, cutoffs were calculated using a formula, =BYCOL($C$2:$O$18, LAMBDA(c, QUARTILE.INC(FILTER(c, c<>0), 3))), that identifies the 75th percentile (Q3) for each centrality measure while excluding zeros, and conditional formatting was applied to highlight modules exceeding these thresholds. This revealed that the same top modules repeatedly concentrate high values across multiple centrality measures. Notably, modules 11, 7, 1, 8, and 2 emerge together across degree, strength, betweenness, closeness, and eigenvector statistics. This indicates that they are not only hubs (high degree) but also bridges (betweenness) and globally central clusters (eigenvector), making them structurally dominant groups likely to coordinate functionally important pathways.
• Our network analysis therefore highlights these five modules and their hub genes—CCR7, PTPRC, SPRR1B, IL1B, and MMP9—as biologically coherent drivers of atopic dermatitis. CCR7 and PTPRC point to immune cell recruitment and activation, IL1B and MMP9 represent classic mediators of tissue inflammation and remodelling, while SPRR1B reflects barrier dysfunction through its role in epidermal differentiation. Together, these hub modules capture the dual pathology of AD: immune dysregulation and barrier compromise, reinforcing their relevance as candidate biomarkers and potential therapeutic targets.
• Strikingly, the same key players that were prioritised at the per-gene level (Figure 10.1.9.1) reappear here as module hubs (Figure 10.1.9.2): CCR7 directs T-cell and dendritic-cell trafficking, PTPRC (CD45) regulates T-cell activation signalling, and IL1B encodes a pro-inflammatory cytokine central to the initiation and persistence of skin inflammation. Their recurrence across both analyses underscores their functional centrality in atopic dermatitis.
• Other highly ranked modules add complementary biology. Module 1 is led by SPRR1B, a skin barrier protein often dysregulated in barrier-defective states, while Module 2 centres on MMP9, a matrix metalloproteinase that drives extracellular matrix breakdown and tissue remodelling. These hubs point directly to the barrier dysfunction and tissue-remodelling aspects of atopic dermatitis, complementing the immune-driven mechanisms highlighted by CCR7, PTPRC, and IL1B.
• Additionally, among these, Module 11 stands out as the most structurally dominant. With hub gene CCR7, it shows, by far, the highest mean and median values for multiple metrics: degree, normalised degree, strength, and eigenvector centrality. This consistency across both mean- and median-based statistics highlights that Module 11 is not only driven by a few extreme hubs but is broadly composed of highly connected and influential genes. Its top hub, CCR7, therefore emerges as a particularly important node—representing a gene that anchors the most central and cohesive module of the network. In the disease context, this reinforces CCR7’s pivotal role in directing immune-cell trafficking into lesional skin and positions Module 11 as the key immune-anchored cluster in atopic dermatitis pathogenesis.
• Altogether, the module-level analysis reveals that the most central communities of the network align closely with the dual pathology of atopic dermatitis: immune hyperactivation anchored by cytokines and T-cell regulators, alongside barrier compromise and tissue remodelling driven by structural and enzymatic proteins. By examining not just individual nodes but also the modules that organise them, this file highlights coordinated groups of genes that may act as collective drivers of disease biology, strengthening their candidacy as biomarkers and potential therapeutic entry points.

• The following code (Code Block 10.1.9) implements these calculations step by step, beginning with the definition of edge distances for path-based centralities. The STRING database gives each protein–protein interaction a confidence score stored as E(g)$combined_score, where larger values mean stronger evidence. For centrality measures like betweenness and closeness, we need distances where stronger edges are shorter. The code achieves this with E(g)$dist <- 1 / pmax(E(g)$combined_score, .Machine$double.eps) .
• Here, 1 / … inverts the score, so that higher confidence produces a shorter distance, while pmax(…, .Machine$double.eps) prevents division by zero by replacing any zero or near-zero score with .Machine$double.eps, the smallest positive number R can store.
• With distances set, the code calculates a panel of node centralities. First, degree(g) measures how many direct edges each node has. This is also normalized (deg_norm) by dividing by the maximum possible degree (vcount(g) - 1), so values fall between 0 and 1. Next, strength(g, weights = E(g)$combined_score) computes the weighted degree, where edges contribute in proportion to their STRING confidence, so nodes with many strong links stand out.
• Path-based centralities follow. betweenness(g, weights = E(g)$dist, normalized = TRUE) counts how often a node lies on shortest paths, where “shortest” is defined by the high-confidence edges. closeness(g, weights = E(g)$dist, normalized = TRUE) scores how close a node is on average to all others; if the graph is disconnected, infinite values are replaced with NA. Finally, eigen_centrality(g, weights = E(g)$combined_score) evaluates a node’s influence based not only on its own connections but also on whether it connects to other highly connected nodes.
• The results are combined into a tidy node-level data frame, storing gene identifiers (STRING_id, gene), module membership, centrality scores (degree, normalized degree, weighted strength, betweenness, closeness, eigenvector), and also the gene’s log-fold change (lfc) and adjusted p-value (padj). Nodes are sorted within each module by degree so hubs appear first, and the file is exported as ppi_node_centrality.csv.
• Finally, the code summarizes these values at the module level. It records the module size (n_genes), calculates mean centralities and median centralities across genes in each module, and identifies the top hub gene by degree. These summaries are merged into one table, ordered by descending module size, and saved as ppi_module_centrality_summary.csv.

Code Block 10.1.9. Centrality Metrics & Module-Level Summaries — Degree/strength, betweenness/closeness/eigenvector, CSV exports

			
		 

# ==========================================================
# [Section 9/9] Centrality metrics & module-level summaries (CSVs)
#    - Degree/strength; betweenness/closeness on 1/score; eigenvector
# ==========================================================
# 
 
# 9.1 Define edge distance for shortest paths
# Higher STRING score = stronger/closer connection, so use 1/score as distance.
E(g)$dist <- 1 / pmax(E(g)$combined_score, .Machine$double.eps)

# 9.2 Compute a small panel of standard centralities:
# - Degree (unweighted) and normalized degree
# - degree (local connectivity)
deg_unw   <- degree(g)                       # already in V(g)$deg, recompute for clarity
# degree_norm (scaled 0–1)
deg_norm  <- if (vcount(g) > 1) deg_unw / (vcount(g) - 1) else deg_unw

# - Strength (weighted degree by STRING score)
strength_w <- strength(g, weights = E(g)$combined_score)

# - Betweenness (weighted by distance)
btw_w <- betweenness(g, weights = E(g)$dist, directed = FALSE, normalized = TRUE)

# - Closeness (weighted by distance)
clo_w <- closeness(g, weights = E(g)$dist, normalized = TRUE)
clo_w[!is.finite(clo_w)] <- NA  # guard against Inf/NaN in disconnected cases

# - Eigenvector centrality (use score as weight, not distance/ influence via well-connected neighbors)
eig_w <- eigen_centrality(g, weights = E(g)$combined_score,
                          directed = FALSE, scale = TRUE)$vector


# 9.3 Export node-level centrality table
# Assemble a tidy per-node centrality table for export and ranking/plotting in downstream tools.
cent_df <- data.frame(
  STRING_id     = V(g)$name,
  gene          = V(g)$label,
  module        = as.character(V(g)$module),
  degree        = deg_unw,
  degree_norm   = deg_norm,
  strength_w    = strength_w,
  betweenness_w = btw_w,
  closeness_w   = clo_w,
  eigenvector   = eig_w,
  lfc           = V(g)$lfc,
  padj          = V(g)$padj,
  stringsAsFactors = FALSE
)

# Sort within each module by degree so the top hubs float to the top of each group in the CSV.
cent_df <- cent_df[order(as.numeric(cent_df$module), -cent_df$degree), ]

# Save node-level centralities
write.csv(cent_df, file.path(fig_dir, "ppi_node_centrality.csv"), row.names = FALSE)


# 9.4 Export module-level summary (size, means, medians, top hub (by degree)
# Module size
mod_sizes <- as.data.frame(table(cent_df$module), stringsAsFactors = FALSE)
names(mod_sizes) <- c("module", "n_genes")


# Means across nodes in each module
mod_means <- aggregate(
  cent_df[, c("degree","degree_norm","strength_w","betweenness_w","closeness_w","eigenvector")],
  by = list(module = cent_df$module),
  FUN = function(x) mean(x, na.rm = TRUE)
)
names(mod_means)[-1] <- paste0(names(mod_means)[-1], "_mean")

# Medians across nodes in each module
mod_meds <- aggregate(
  cent_df[, c("degree","degree_norm","strength_w","betweenness_w","closeness_w","eigenvector")],
  by = list(module = cent_df$module),
  FUN = function(x) median(x, na.rm = TRUE)
)
names(mod_meds)[-1] <- paste0(names(mod_meds)[-1], "_median")

# Identify a simple top hub per module (by highest degree). Other choices (e.g., strength) are possible.
top_hub <- tapply(seq_len(nrow(cent_df)), cent_df$module, function(idx) {
  i <- idx[which.max(cent_df$degree[idx])]
  cent_df$gene[i]
})
top_hub_df <- data.frame(module = names(top_hub), top_hub_gene = unname(unlist(top_hub)),
                         stringsAsFactors = FALSE)

# Merge all summaries
module_centrality_summary <- Reduce(function(x, y) merge(x, y, by = "module", all = TRUE),
                                    list(mod_sizes, mod_means, mod_meds, top_hub_df))

# Order modules by descending size (then numeric id) in the final summary for quick, consistent browsing.
module_centrality_summary$module_num <- as.numeric(module_centrality_summary$module)
module_centrality_summary <- module_centrality_summary[order(-module_centrality_summary$n_genes,
                                                             module_centrality_summary$module_num), ]
module_centrality_summary$module_num <- NULL

# Save per-module summaryLondon Luton Airport
write.csv(module_centrality_summary,
          file.path(fig_dir, "ppi_module_centrality_summary.csv"), row.names = FALSE)


		    

10.2. GSEA for Gene Ontology (GO) Enrichment (BP, MF, CC)

• TEXTTTTTTTTTTT FOR 10.2.

o 10.3. ORA for KEGG Pathway Enrichment

• TEXTTTTTTTTTTT FOR 10.3.

6. Create DESeqDataSet and Run Differential Expression Analysis

We create a DESeqDataSet object, normalize the data, and run the differential expression analysis using DESeq2. This model tests for differences in gene expression between the Lesional and Non-Lesional groups.

7. Extract and Process Results

We extract the results from DESeq2, separating the differentially expressed genes (DEGs) into upregulated and downregulated categories, and prepare the top genes for further analysis.

8. General Data Visualization

We generate various plots to visualize the data and the results of the differential expression analysis.

8.1 Box Plot

8.2 Dispersion Estimates

8.3 UMAP Plot

8.4 Histogram of P-values

8.5 Volcano Plot

8.6 MA Plot

9. Functional Annotation and Pathway Analysis

After identifying DEGs, we perform GO and KEGG pathway enrichment analyses to understand the biological significance of the results.

9.1 Gene Ontology (GO) Enrichment Analysis

9.2 KEGG Pathway Enrichment Analysis

10. Venn Diagram for Gene Sets

We visualize the overlap between different gene sets (total, upregulated, downregulated) using a Venn diagram.

11. Extract Top 20 Genes and Print Results

Finally, we extract the top 20 differentially expressed genes from both categories and print the results.

Conclusion

This complete pipeline guides you through performing differential gene expression analysis on the GSE157194 dataset using DESeq2. The analysis identifies significant gene expression changes between lesional and non-lesional skin samples from Atopic Dermatitis patients, and includes essential data visualizations and enrichment analyses to interpret the results biologically.