The famous human airway smooth muscle dataset is called in. The 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). Because these are raw counts rather than statistical results, differential expression analysis must first be performed before creating a volcano plot.

## class: RangedSummarizedExperiment 
## dim: 63677 8 
## metadata(1): ''
## assays(1): counts
## rownames(63677): ENSG00000000003 ENSG00000000005 ... ENSG00000273492
##   ENSG00000273493
## rowData names(10): gene_id gene_name ... seq_coord_system symbol
## colnames(8): SRR1039508 SRR1039509 ... SRR1039520 SRR1039521
## colData names(9): SampleName cell ... Sample BioSample

Dexamethasone is the treatment variable. The DESeq() function estimates size factors for normalization, gene-wise dispersion, 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 α = 0.05 were applied independently to every test, approximately 5% of genes could appear significant purely by random chance, resulting in thousands of false positive discoveries. To reduce this problem, DESeq2 applies the Benjamini–Hochberg procedure, which controls the False Discovery Rate (FDR). The adjusted p-values (padj) therefore provide a more reliable measure of statistical significance than the raw p-values when many hypothesis tests are performed simultaneously.

Since the airway dataset contains 64,102 genes, using unadjusted p-values would be expected to produce approximately 64,102 × 0.05 ≈ 3,205 false positives under the null hypothesis.

The results() function returns a table containing one row for each gene. Important variables include the estimated log₂ fold change (log2FoldChange), the hypothesis test p-value (pvalue), and the Benjamini–Hochberg adjusted p-value (padj). The latter accounts for the large number of simultaneous hypothesis tests performed across all genes.

## 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 table was formatted as a data frame and cleaned by removing rows with missing values prior to plotting.

res_df <- as.data.frame(results)
res_df <- na.omit(res_df)

The Significant variable is created to classify genes as statistically significant or not. Two conditions are to be satisfied: an absolute log₂ fold change >1 (indicating at least a two-fold increase or decrease in expression) and an adjusted p-value <0.05 (indicating statistical significance after controlling the false discovery rate). This is later used only 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₂ fold change, where positive values indicate genes that are upregulated in the dexamethasone-treated samples and negative values indicate genes that are downregulated. The vertical y-axis displays −log₁₀(Benjamini–Hochberg adjusted p-value), so genes with stronger statistical evidence appear higher on the graph.

X-axis (log₂ fold change): Effect size (how much expression changed) Left = downregulated genes Right = upregulated genes