The airway (human smooth muscle) dataset contains raw RNA sequencing (RNA-seq) read counts measured for genes across eight (8) samples. Four (4) samples were untreated while four (4) were treated with dexamethasone (dex) in vitro. Dexamethasone is a glucocorticoid steroid hormone that binds the glucocorticoid receptor, and the receptor-dex complex moves into the nucleus. This compelx then binds short DNA sequences known as glucocorticoid response elements (GREs) near specific genes, altering transcription rates. Because these are raw counts rather than statistical results, differential expression analysis must first be performed.

Dexamethasone is the treatment variable. The DESeq() function estimates size factors for normalization and fits a negative binomial generalized linear model for each gene. It also performs hypothesis testing to determine whether gene expression differs significantly between treated and untreated samples.

Method: Differential expression analysis simultaneously tests more than 64,000 genes. If a significance level of \(\alpha = 0.05\) were applied independently to every test, approximately 5% of genes could appear significant purely by random chance, resulting in ~3,205 (64,102*0.05) false positive discoveries. To address this, DESeq2 applies the Benjamini–Hochberg procedure, which controls the False Discovery Rate (FDR). The adjusted p-values (padj) provide a more reliable measure of statistical significance than the raw p-values when many hypothesis tests are performed simultaneously.

Variance Stabilizing transformation applied, followed by principal component analysis(PCA). This is because the variance of RNA-seq counts increases with the mean (heteroscedasticity). This transforms the counts so that variance is approximately constant across expression levels, making PCA more meaningful.

PCA shows samples cluster by treatment (dex).

The results() function returns a table containing one row for each gene. Important variables include the estimated log-base2 fold change (log2FoldChange), the hypothesis test p-value (pvalue), and the Benjamini–Hochberg adjusted p-value (padj).

## log2 fold change (MLE): dex untrt vs trt 
## Wald test p-value: dex untrt vs trt 
## DataFrame with 6 rows and 6 columns
##                   baseMean log2FoldChange     lfcSE      stat    pvalue
##                  <numeric>      <numeric> <numeric> <numeric> <numeric>
## ENSG00000000003 708.602170      0.3788470  0.173141  2.188082 0.0286636
## ENSG00000000005   0.000000             NA        NA        NA        NA
## ENSG00000000419 520.297901     -0.2037604  0.100599 -2.025478 0.0428183
## ENSG00000000457 237.163037     -0.0340428  0.126279 -0.269584 0.7874802
## ENSG00000000460  57.932633      0.1171786  0.301237  0.388992 0.6972820
## ENSG00000000938   0.318098      1.7245505  3.493633  0.493627 0.6215698
##                      padj
##                 <numeric>
## ENSG00000000003  0.139308
## ENSG00000000005        NA
## ENSG00000000419  0.183359
## ENSG00000000457  0.930572
## ENSG00000000460  0.895441
## ENSG00000000938        NA

Results were converted to a data frame and cleaned by removing rows with missing values prior to plotting.

result_df <- as.data.frame(result)
result_df <- na.omit(result_df)

The Significant variable is created to classify genes as statistically significant or not. Two conditions are to be satisfied: an absolute log-base2 fold change >1 (indicating at least a two-fold increase or decrease in expression) and an adjusted p-value <0.05. This is used to determine only the color of each point in the volcano plot.

Figure 1. The volcano plot visualizes both the magnitude of gene expression changes and their statistical significance. The horizontal x-axis displays the log-base2 fold change in dexamethasone-treated samples, where positive values indicate genes that are upregulated and negative values indicate genes that are downregulated. The vertical y-axis displays -log-base10 (Benjamini–Hochberg adjusted p-value), so genes with stronger statistical evidence appear higher on the graph.

X-axis (log-base2 fold change): Effect size (how much expression changed)

Gene counts were normalized and then selected the top 50 differentially expressed genes by the adjusted p-values (padj).

Extracted the expression values for the selected genes and standardized each gene for heatmap visualization.

Created sample annotations for dexamethasone treatment and cell line.

Generated a clustered heatmap with sample annotations. Interpretation: Each column is an RNA-Seq sample, and each row represents one of the top differentially expressed genes. Red-blue scale shows relative expression, where the color is relative the gene’s average (due to standardization performed earlier.)

The dendrogram groups these samples into two main branches, which indicates treatment is the dominant source of variation. It appears the genes separate into two major clusters: The upper gene cluster is highly expressed in untreated samples (upper left) and low in treated samples (upper right). Conversely, the lower gene cluster shows the opposite pattern: low expression in untreated (bottom left) and high expression in treated (bottom right). Cell line on the right side appears to have less influence, because it is mixed within each cluster. This suggests dexamethasone treatment explains more variation than cell line for these genes.