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.
• 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/omidm/OneDrive/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")
	# ============================
	    

# Load Required Libraries
library(DESeq2)
library(data.table)
library(umap)
library(ggplot2)
library(clusterProfiler)
library(org.Hs.eg.db)
library(VennDiagram)
library(knitr)
library(tinytex)
            

# Load counts table from GEO
urld <- "https://www.ncbi.nlm.nih.gov/geo/download/?format=file&type=rnaseq_counts"
path <- paste(urld, "acc=GSE157194", "file=GSE157194_raw_counts_GRCh38.p13_NCBI.tsv.gz", sep="&")
tbl <- as.matrix(data.table::fread(path, header=T, colClasses="integer"), rownames="GeneID")

# Load gene annotations 
apath <- paste(urld, "type=rnaseq_counts", "file=Human.GRCh38.p13.annot.tsv.gz", sep="&")
annot <- data.table::fread(apath, header=T, quote="", stringsAsFactors=F, data.table=F)
rownames(annot) <- annot$GeneID
            

# Sample selection
gsms <- paste0("01101001010101011010011001011001011001101001011010",
               "01100101011010101001011001010101010010101010101001",
               "01100100110XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX",
               "XXXXXXXXXXXXXXXX")
sml <- strsplit(gsms, split="")[[1]]

# Filter out excluded samples (marked as "X")
sel <- which(sml != "X")
sml <- sml[sel]
tbl <- tbl[ ,sel]

# Group membership for samples
gs <- factor(sml)
groups <- make.names(c("Lesional","Non-Lesional"))
levels(gs) <- groups
sample_info <- data.frame(Group = gs, row.names = colnames(tbl))
            

# Pre-filter low count genes
keep <- rowSums(tbl >= 10) >= min(table(gs))
tbl <- tbl[keep, ]
            

# PCA Plot
vsd <- vst(ds, blind=FALSE)
pcaData <- plotPCA(vsd, intgroup="Group", returnData=TRUE)
percentVar <- round(100 * attr(pcaData, "percentVar"))

ggplot(pcaData, aes(PC1, PC2, color=Group)) +
  geom_point(size=3) +
  xlab(paste0("PC1: ", percentVar[1], "% variance")) +
  ylab(paste0("PC2: ", percentVar[2], "% variance")) +
  coord_fixed() +
  theme_minimal()
            

# Create DESeqDataSet and run DESeq2
ds <- DESeqDataSetFromMatrix(countData=tbl, colData=sample_info, design= ~Group)
ds <- DESeq(ds, test="Wald", sfType="poscount")
            

# Extract results for top genes
r <- results(ds, contrast = c("Group", groups[1], groups[2]), alpha = 0.01, pAdjustMethod = "fdr")

# Separate downregulated and upregulated DEGs
downregulated <- r[r$log2FoldChange < 0 & r$padj < 0.01, ]
upregulated <- r[r$log2FoldChange > 0 & r$padj < 0.01, ]

# Take the top 10 genes for downregulated and upregulated
top_downregulated <- downregulated[order(downregulated$log2FoldChange)[1:10], ]
top_upregulated <- upregulated[order(upregulated$log2FoldChange, decreasing = TRUE)[1:10], ]

# Merge with gene annotations
top_downregulated <- merge(as.data.frame(top_downregulated), annot, by = 0, sort = FALSE)
top_upregulated <- merge(as.data.frame(top_upregulated), annot, by = 0, sort = FALSE)

# Round "padj" and "log2FoldChange" columns for readability
top_downregulated <- within(top_downregulated, {
  padj <- sprintf("%.4f", round(padj, 2))
  log2FoldChange <- sprintf("%.4f", round(log2FoldChange, 2))
})

top_upregulated <- within(top_upregulated, {
  padj <- sprintf("%.4f", round(padj, 2))
  log2FoldChange <- sprintf("%.4f", round(log2FoldChange, 2))
})

# Print the top_downregulated table using kable
kable(top_downregulated, format = "markdown")
kable(top_upregulated, format = "markdown")
            

# Box-and-whisker plot
dat <- log10(counts(ds, normalized = T) + 1)
ord <- order(gs)
palette(c("#1B9E77", "#7570B3"))
par(mar=c(7,4,2,1))
boxplot(dat[,ord], boxwex=0.6, notch=T, main="GSE157194", ylab="lg(norm.counts)", outline=F, las=2, col=gs[ord])
legend("topleft", groups, fill=palette(), bty="n")
            

# Dispersion estimates plot
plotDispEsts(ds, main="GSE157194 Dispersion Estimates")
            

# UMAP plot
dat <- dat[!duplicated(dat), ]
par(mar=c(3,3,2,6), xpd=TRUE, cex.main=1.5)
ump <- umap(t(dat), n_neighbors = 15, random_state = 123)
plot(ump$layout, main="UMAP plot, nbrs=15", xlab="", ylab="", col=gs, pch=20, cex=1.5)
legend("topright", inset=c(-0.15,0), legend=groups, pch=20, col=1:length(groups), title="Group", pt.cex=1.5)
            

# Histogram of adjusted p-values
hist(r$padj, breaks=seq(0, 1, length = 21), col = "grey", border = "white", 
     xlab = "", ylab = "", main = "GSE157194 Frequencies of padj-values")
            

# Volcano plot
old.pal <- palette(c("#00BFFF", "#FF3030"))
par(mar=c(4,4,2,1), cex.main=1.5)
plot(r$log2FoldChange, -log10(r$padj), main=paste(groups[1], "vs", groups[2]),
     xlab="log2FC", ylab="-log10(Padj)", pch=20, cex=0.5)
with(subset(r, padj<0.01 & abs(log2FoldChange) >= 1.5),
     points(log2FoldChange, -log10(padj), pch=20, col=(sign(log2FoldChange) + 3)/2, cex=1))
legend("bottomleft", title=paste("Padj<", 0.01, sep=""), legend=c("down", "up"), pch=20,col=1:2)
            

# MA plot
par(mar=c(4,4,2,1), cex.main=1.5)
plot(log10(r$baseMean), r$log2FoldChange, main=paste(groups[1], "vs", groups[2]),
     xlab="log10(mean of normalized counts)", ylab="log2FoldChange", pch=20, cex=0.5)
with(subset(r, padj<0.01 & abs(log2FoldChange) >= 1.5),
     points(log10(baseMean), log2FoldChange, pch=20, col=(sign(log2FoldChange) + 3)/2, cex=1))
legend("bottomleft", title=paste("Padj<", 0.01, sep=""), legend=c("down", "up"), pch=20,col=1:2)
abline(h=0)
palette(old.pal)
            

# GO enrichment analysis for upregulated and downregulated genes
up_genes <- top_upregulated$GeneID
down_genes <- top_downregulated$GeneID

# GO enrichment for upregulated genes
ego_up <- enrichGO(gene = up_genes, OrgDb = org.Hs.eg.db, keyType = "ENSEMBL",
                   ont = "BP", pAdjustMethod = "fdr", qvalueCutoff = 0.05)

# GO enrichment for downregulated genes
ego_down <- enrichGO(gene = down_genes, OrgDb = org.Hs.eg.db, keyType = "ENSEMBL",
                     ont = "BP", pAdjustMethod = "fdr", qvalueCutoff = 0.05)

# Plotting the top GO terms
barplot(ego_up, showCategory=10, title="Top GO terms for Upregulated Genes", drop=TRUE)
barplot(ego_down, showCategory=10, title="Top GO terms for Downregulated Genes", drop=TRUE)
            

# KEGG pathway enrichment analysis for upregulated and downregulated genes
kegg_up <- enrichKEGG(gene = up_genes, organism = "hsa", pAdjustMethod = "fdr", qvalueCutoff = 0.05)
kegg_down <- enrichKEGG(gene = down_genes, organism = "hsa", pAdjustMethod = "fdr", qvalueCutoff = 0.05)

# Plotting the top KEGG pathways
barplot(kegg_up, showCategory=10, title="Top KEGG Pathways for Upregulated Genes")
barplot(kegg_down, showCategory=10, title="Top KEGG Pathways for Downregulated Genes")
            

# Venn diagram for gene sets
total_genes <- rownames(tbl)
upregulated_genes <- rownames(top_upregulated)
downregulated_genes <- rownames(top_downregulated)

venn_data <- list(
  Total = total_genes,
  Upregulated = upregulated_genes,
  Downregulated = downregulated_genes
)

venn_plot <- venn.diagram(
  x = venn_data,
  category.names = c("Total", "Upregulated", "Downregulated"),
  filename = NULL,
  margin = 0,
  cex.prop = 0.8,
  offset.prop = 0.2
)

grid.draw(venn_plot)
            

# Extract results for top genes table (rounded)
r <- results(ds, contrast = c("Group", groups[1], groups[2]), alpha = 0.01, pAdjustMethod = "fdr")

# Separate downregulated and upregulated DEGs
downregulated <- r[r$log2FoldChange < 0 & r$padj < 0.01, ]
upregulated <- r[r$log2FoldChange > 0 & r$padj < 0.01, ]

# Take the top 20 genes for downregulated and upregulated
top_downregulated <- downregulated[order(downregulated$log2FoldChange)[1:20], ]
top_upregulated <- upregulated[order(upregulated$log2FoldChange, decreasing = TRUE)[1:20], ]

# Merge with gene annotations
top_downregulated <- merge(as.data.frame(top_downregulated), annot, by = 0, sort = FALSE)
top_upregulated <- merge(as.data.frame(top_upregulated), annot, by = 0, sort = FALSE)

# Round "padj" and "log2FoldChange" columns
top_downregulated <- within(top_downregulated, {
  padj <- sprintf("%.4f", round(padj, 4))
  log2FoldChange <- sprintf("%.4f", round(log2FoldChange, 4))
 })

top_upregulated <- within(top_upregulated, {
  padj <- sprintf("%.4f", round(padj, 4))
  log2FoldChange <- sprintf("%.4f", round(log2FoldChange, 4))
})

# Print the top_downregulated table using kable
kable(top_downregulated, format = "markdown") 
kable(top_upregulated, format = "markdown")