1. Introduction

This report analyzes TPM expression profiles across nine Arabidopsis thaliana seed developmental stages (E-GEOD-94457; Narsai et al., 2017), spanning dry seed, seed development, germination, and stratification phases. Using three biological replicates per stage, the workflow integrates PCA, limma-based differential expression testing, unsupervised clustering, and GO enrichment to characterize transcriptional reprogramming during seed development.

2. Methodology

2.1 Data Import and Library

Data were downloaded from the Expression Atlas repository (Papatheodorou et al., 2020; http://ftp.ebi.ac.uk/pub/databases/microarray/data/atlas/experiments/E-GEOD-94457/), originally generated by Narsai et al. (2017), and imported into R along with the necessary libraries. The selection of samples g1 through g9 ensured that the PCA comparison focused on the seed developmental trajectory, excluding outliers associated with post-germination growth conditions. This study utilized only three of the five replicates, as this was sufficient to represent the biological population.

library(dplyr)
library(tidyr)
library(ggplot2)
library(limma)
library(ComplexHeatmap)
library(circlize)
library(grid)
library(ggrepel)
library(clusterProfiler)
library(org.At.tair.db)
library(enrichplot)

file_path <- "D:\\OneDrive - Universitas Islam Indonesia\\College\\R for Transcriptomics\\Final Exam\\E-GEOD-94457-tpms.tsv"
expr_data <- read.table(file_path, header = TRUE, row.names = 1, sep = "\t")

# Ambil grup g1 sampai g9
selected_groups <- c("g1", "g2", "g3", "g4", "g5", "g6", "g7", "g8", "g9")
expr_data_subset <- expr_data[, selected_groups]

# Fungsi pemecah replikat (3 Replikat)
expand_column_safe <- function(col_name, vec) {
  split_list <- strsplit(as.character(vec), ",")
  max_rep <- 3

  padded_mat <- t(sapply(split_list, function(x) {
    length(x) <- max_rep
    as.numeric(x)
  }))

  df_res <- as.data.frame(padded_mat)
  colnames(df_res) <- paste0(col_name, "_rep", 1:max_rep)
  return(df_res)
}

expr_expanded_list <- lapply(selected_groups, function(col) {
  expand_column_safe(col, expr_data_subset[[col]])
})

expr_expanded <- do.call(cbind, expr_expanded_list)
rownames(expr_expanded) <- rownames(expr_data_subset)

2.2 Data Transformation and Filtering

Expression values were log2-transformed via \(\log_2(x+1)\), then genes with zero variance across samples were filtered out to remove uninformative and non-varying transcripts.

expr_log_rep <- log2(expr_expanded + 1)

gene_var <- apply(expr_log_rep, 1, var, na.rm = TRUE)
clean_log_rep <- expr_log_rep[gene_var > 0 & !is.na(gene_var), ]

cat("Dimensi sebelum filtering:", dim(expr_log_rep), "\n")
## Dimensi sebelum filtering: 21020 27
cat("Dimensi setelah filtering:", dim(clean_log_rep), "\n")
## Dimensi setelah filtering: 16100 27
head(clean_log_rep[, 1:6])
##             g1_rep1   g1_rep2   g1_rep3   g2_rep1   g2_rep2   g2_rep3
## AT1G01010 2.0000000 2.0000000 2.0000000 2.3219281 2.3219281 2.3219281
## AT1G01020 2.0000000 2.3219281 2.3219281 2.3219281 2.3219281 2.3219281
## AT1G01030 0.6780719 0.7655347 0.7655347 0.7655347 0.8479969 0.9259994
## AT1G01040 1.0000000 1.0000000 1.0000000 2.0000000 2.0000000 2.0000000
## AT1G01046 0.0000000 0.1375035 0.2630344 0.2630344 0.3785116 0.4854268
## AT1G01050 3.8073549 3.8073549 3.8073549 4.1699250 4.1699250 4.1699250

2.3 Principal Component Analysis

PCA (Jolliffe & Cadima, 2016) reduces dimensionality via eigendecomposition of the covariance matrix, \(\text{Cov}(X) = \frac{1}{n-1}X^TX\), projecting samples onto components \(Z = XW\) that capture maximal variance.

2.4 Differential Expression Analysis

Differential expression was tested with limma (Ritchie et al., 2015) using Bonferroni (\(p_{adj}=p \times m\); Dunn, 1961, family-wise error) and Benjamini-Hochberg (\(q_{(k)}=\min_{j \ge k}\frac{m}{j}p_{(j)}\), Benjamini & Hochberg, 1995), thresholded by fold change, \(\log_2FC=\log_2(g2/g1)\).

2.5 Clustering

Bonferroni-significant genes (adj. p<0.05, |log2FC|>1) were chosen for clustering for their conservative, high-confidence threshold. K-means (MacQueen, 1967) minimizes within-cluster variance (\(\sum_k\sum_{x_i \in C_k}||x_i-\mu_k||^2\)), while Ward’s method (Ward, 1963) merges clusters minimizing the increase in total within-cluster variance.

2.6 Functional Enrichment

GO Biological Process enrichment (clusterProfiler::enrichGO; Wu et al., 2021, based on the Gene Ontology; Ashburner et al., 2000) using org.At.tair.db was performed separately on K-means Cluster 8, Ward Cluster 9, and their combined gene set.

3. Results

3.1 Principal Component Analysis

pca_rep_res <- prcomp(t(clean_log_rep), center = TRUE, scale. = FALSE)
var_explained <- round(100 * (pca_rep_res$sdev^2 / sum(pca_rep_res$sdev^2)), 2)

pca_rep_df <- as.data.frame(pca_rep_res$x)
pca_rep_df$sample_id <- rownames(pca_rep_df)
pca_rep_df$group <- gsub("_rep.*", "", pca_rep_df$sample_id)

stage_map <- c(
  "g1" = "Dry Seed",
  "g2" = "Seed Development",
  "g3" = "Germination 1h",
  "g4" = "Germination 12h",
  "g5" = "Germination 24h",
  "g6" = "Germination 6h",
  "g7" = "Stratification 1h",
  "g8" = "Stratification 12h",
  "g9" = "Stratification 48h"
)
pca_rep_df$Stage <- stage_map[pca_rep_df$group]

ggplot(pca_rep_df, aes(x = PC1, y = PC2, color = Stage)) +
  geom_point(size = 3.5, alpha = 0.85) +
  geom_text_repel(aes(label = sample_id), size = 3, max.overlaps = 15) +
  labs(
    title = "PCA of Arabidopsis Seed Developmental Stages (3 Replicates)",
    x = paste0("PC1 (", var_explained[1], "% Variance)"),
    y = paste0("PC2 (", var_explained[2], "% Variance)")
  ) +
  theme_minimal() +
  theme(plot.title = element_text(face = "bold", size = 13))

The PCA plot shows two proximate clusters sharing the same quadrant: one grouping dry seed, seed development, and early stratification (1h); the other grouping late-stage stratification with early (1h) and late (24h) germination.

3.2 Differential Expression Analysis

Here, you will find the difference between the Benjamini-Hochberg method and the Bonferroni method that divides the p-value by the number of genes, making it more strict. This is what the volcano plot looks like.

selected_samples <- colnames(clean_log_rep)[grep("g1|g2", colnames(clean_log_rep))]
expr_deg_subset <- clean_log_rep[, selected_samples]

groups_deg <- factor(gsub("_rep.*", "", colnames(expr_deg_subset)))
design_deg <- model.matrix(~ 0 + groups_deg)
colnames(design_deg) <- levels(groups_deg)

fit <- lmFit(expr_deg_subset, design_deg)
contrast_matrix <- makeContrasts(g2_vs_g1 = g2 - g1, levels = design_deg)
fit_contrast <- contrasts.fit(fit, contrast_matrix)
fit_bayes <- eBayes(fit_contrast)

# Benjamini-Hochberg (FDR)
deg_bh <- topTable(fit_bayes, coef = "g2_vs_g1", number = Inf, adjust.method = "BH")
deg_bh$Gene <- rownames(deg_bh)
deg_bh <- deg_bh %>%
  mutate(Expression = case_when(
    logFC > 1 & adj.P.Val < 0.05 ~ "Up-regulated",
    logFC < -1 & adj.P.Val < 0.05 ~ "Down-regulated",
    TRUE ~ "Not Significant"
  ))

# Bonferroni
deg_bonf <- topTable(fit_bayes, coef = "g2_vs_g1", number = Inf, adjust.method = "bonferroni")
deg_bonf$Gene <- rownames(deg_bonf)
deg_bonf <- deg_bonf %>%
  mutate(Expression = case_when(
    logFC > 1 & adj.P.Val < 0.05 ~ "Up-regulated",
    logFC < -1 & adj.P.Val < 0.05 ~ "Down-regulated",
    TRUE ~ "Not Significant"
  ))

sig_genes_bonf <- deg_bonf$Gene[deg_bonf$Expression != "Not Significant"]

# Volcano plot: Benjamini-Hochberg
ggplot(deg_bh, aes(x = logFC, y = -log10(adj.P.Val), color = Expression)) +
  geom_point(alpha = 0.5, size = 1.5) +
  scale_color_manual(values = c("Up-regulated" = "#E41A1C",
                                "Down-regulated" = "#377EB8",
                                "Not Significant" = "grey70")) +
  geom_hline(yintercept = -log10(0.05), linetype = "dashed", color = "black") +
  geom_vline(xintercept = c(-1, 1), linetype = "dashed", color = "black") +
  labs(
    title = "Volcano Plot: Seed Development (g2) vs Dry Seed (g1)",
    subtitle = "Adjustment Method: Benjamini-Hochberg (FDR < 0.05, |Log2FC| > 1)",
    x = "Log2 Fold Change",
    y = "-Log10 Adjusted P-Value (BH)"
  ) +
  theme_minimal() +
  theme(plot.title = element_text(face = "bold", size = 12))

Based on the Benjamini-Hochberg plot, a greater number of seed development genes appear to be upregulated than downregulated to dry seeds. A “V” also appears, representing genes that passed the initial filter (filter 0) but did not qualify as differentially expressed genes (DEGs).

# Volcano plot: Bonferroni
ggplot(deg_bonf, aes(x = logFC, y = -log10(adj.P.Val), color = Expression)) +
  geom_point(alpha = 0.5, size = 1.5) +
  scale_color_manual(values = c("Up-regulated" = "#E41A1C",
                                "Down-regulated" = "#377EB8",
                                "Not Significant" = "grey70")) +
  geom_hline(yintercept = -log10(0.05), linetype = "dashed", color = "black") +
  geom_vline(xintercept = c(-1, 1), linetype = "dashed", color = "black") +
  labs(
    title = "Volcano Plot: Seed Development (g2) vs Dry Seed (g1)",
    subtitle = "Adjustment Method: Bonferroni (FWER < 0.05, |Log2FC| > 1)",
    x = "Log2 Fold Change",
    y = "-Log10 Adjusted P-Value (Bonferroni)"
  ) +
  theme_minimal() +
  theme(plot.title = element_text(face = "bold", size = 12))

With the Bonferroni correction, the number of DEGs appears lower than with the Benjamini method. Despite the stringency of the Bonferroni correction, DEGs were still identified; therefore, we will use these Bonferroni-derived DEGs for the next step.

3.3 Clustering

ordered_samples <- colnames(expr_deg_subset)

header_labels <- factor(
  ifelse(gsub("_rep.*", "", ordered_samples) == "g1", "Dry Seed", "Seed Development"),
  levels = c("Dry Seed", "Seed Development")
)

sample_info <- data.frame(
  Sample = ordered_samples,
  Group = header_labels,
  stringsAsFactors = FALSE
)

log2_tpm_mat <- as.matrix(expr_deg_subset[sig_genes_bonf, ordered_samples])

# --- K-Means ---
set.seed(123)
km_res_log2 <- kmeans(log2_tpm_mat, centers = 10, nstart = 25)

df_log2 <- as.data.frame(log2_tpm_mat)
df_log2$Gene_id <- rownames(df_log2)

km_log2_info <- data.frame(
  Gene_id = names(km_res_log2$cluster),
  K_Cluster = paste0("K", km_res_log2$cluster),
  stringsAsFactors = FALSE
)

summary_log2 <- df_log2 %>%
  inner_join(km_log2_info, by = "Gene_id") %>%
  group_by(K_Cluster) %>%
  summarise(across(all_of(ordered_samples), mean, na.rm = TRUE)) %>%
  as.data.frame()

rownames(summary_log2) <- summary_log2$K_Cluster
mat_summary_log2 <- as.matrix(summary_log2[, ordered_samples])
mat_summary_log2 <- mat_summary_log2[intersect(paste0("K", 1:10), rownames(mat_summary_log2)), ]

min_val_log2 <- min(mat_summary_log2, na.rm = TRUE)
max_val_log2 <- max(mat_summary_log2, na.rm = TRUE)
mid_val_log2 <- (min_val_log2 + max_val_log2) / 2

col_fun_log2 <- colorRamp2(
  c(min_val_log2, mid_val_log2, max_val_log2),
  c("#2166AC", "#F7F7F7", "#B2182B")
)

ht_km_log2_summary <- Heatmap(
  mat_summary_log2,
  name = "Log2 TPM",
  col = col_fun_log2,
  cluster_rows = FALSE, cluster_columns = FALSE,
  column_split = sample_info$Group,
  column_title_gp = gpar(fontface = "bold", fontsize = 12),
  cell_fun = function(j, i, x, y, width, height, fill) {
    grid.text(
      sprintf("%.2f", mat_summary_log2[i, j]), x, y,
      gp = gpar(fontsize = 9, fontface = "bold",
                col = ifelse(mat_summary_log2[i, j] > mid_val_log2, "white", "black"))
    )
  },
  show_row_names = TRUE, row_names_gp = gpar(fontface = "bold", fontsize = 10),
  show_column_names = TRUE, column_names_rot = 45, border = TRUE
)

draw(ht_km_log2_summary, column_title = "Expression Profile of 10 K-Means Clusters (3 Replicates)",
     column_title_gp = gpar(fontsize = 14, fontface = "bold"))

Using the k-means algorithm, two interesting gene clusters were identified. Gene Cluster 6 shows high expression during the dry seed phase and low expression during seed development, whereas Gene Cluster 8 shows the opposite pattern. Here, however, we will focus on Cluster 6 to investigate the biological processes associated with the dry seed stage.

# --- Ward's Method ---
dist_mat <- dist(log2_tpm_mat, method = "euclidean")
hc_ward  <- hclust(dist_mat, method = "ward.D2")
ward_clusters <- cutree(hc_ward, k = 10)

ward_log2_info <- data.frame(
  Gene_id = names(ward_clusters),
  W_Cluster = paste0("W", ward_clusters),
  stringsAsFactors = FALSE
)

summary_ward <- df_log2 %>%
  inner_join(ward_log2_info, by = "Gene_id") %>%
  group_by(W_Cluster) %>%
  summarise(across(all_of(ordered_samples), mean, na.rm = TRUE)) %>%
  as.data.frame()

rownames(summary_ward) <- summary_ward$W_Cluster
mat_summary_ward <- as.matrix(summary_ward[, ordered_samples])
hc_summary_ward <- hclust(dist(mat_summary_ward, method = "euclidean"), method = "ward.D2")

ht_ward_log2_summary <- Heatmap(
  mat_summary_ward,
  name = "Log2 TPM",
  col = col_fun_log2,
  cluster_rows = hc_summary_ward, show_row_dend = TRUE, row_dend_width = unit(1.5, "cm"),
  cluster_columns = FALSE,
  column_split = sample_info$Group,
  column_title_gp = gpar(fontface = "bold", fontsize = 12),
  cell_fun = function(j, i, x, y, width, height, fill) {
    grid.text(
      sprintf("%.2f", mat_summary_ward[i, j]), x, y,
      gp = gpar(fontsize = 9, fontface = "bold",
                col = ifelse(mat_summary_ward[i, j] > mid_val_log2, "white", "black"))
    )
  },
  show_row_names = TRUE, row_names_gp = gpar(fontface = "bold", fontsize = 10),
  show_column_names = TRUE, column_names_rot = 45, border = TRUE
)

draw(ht_ward_log2_summary, column_title = "Expression Profile of 10 Ward's Clusters (3 Replicates)",
     column_title_gp = gpar(fontsize = 14, fontface = "bold"))

The same applies to Ward’s algorithm, which considers the within-cluster sum of squares; clusters 6 and 9 are of particular interest. However, we will use cluster 9, which shows high values at the dry seed stage.

3.4 Functional Enrichment

genes_km_k6_clean    <- sub("\\.\\d+$", "", names(km_res_log2$cluster[km_res_log2$cluster == 6]))
genes_ward_w9_clean  <- sub("\\.\\d+$", "", names(ward_clusters[ward_clusters == 9]))
combined_genes_clean <- unique(c(genes_km_k6_clean, genes_ward_w9_clean))

cat("Number of gene in K-Means K6       :", length(genes_km_k6_clean), "\n")
## Number of gene in K-Means K6       : 12
cat("Number of gene in Ward W9          :", length(genes_ward_w9_clean), "\n")
## Number of gene in Ward W9          : 9
cat("Unique Total (K6 + W9)       :", length(combined_genes_clean), "\n\n")
## Unique Total (K6 + W9)       : 12

Here, we will create three biological pathway enrichment functions:

  1. 12 Genes in K-Means Cluster 6
  2. 9 Genes in Ward Cluster 9
  3. 12 Genes Unique in K-Means Cluster 6 and Ward Cluster 9
ego_kmeans_k6 <- enrichGO(
  gene          = genes_km_k6_clean,
  OrgDb         = org.At.tair.db,
  keyType       = "TAIR",
  ont           = "BP",               # Biological Process
  pAdjustMethod = "BH",
  pvalueCutoff  = 0.1,
  qvalueCutoff  = 0.2
)

if (!is.null(ego_kmeans_k6) && nrow(as.data.frame(ego_kmeans_k6)) > 0) {
  print(
    dotplot(ego_kmeans_k6, showCategory = 15, label_format = 30) +
      ggtitle("GO Enrichment (BP): K-Means Cluster 6") +
      theme_bw() +
      theme(
        axis.text.y = element_text(size = 9, color = "black"),
        plot.title  = element_text(face = "bold", size = 12)
      )
  )
  write.csv(as.data.frame(ego_kmeans_k6), "GO_Enrichment_KMeans_Cluster6.csv", row.names = FALSE)
} else {
  cat("-> Tidak ada istilah GO yang signifikan pada K-Means K6.\n")
}

K-means cluster 6 is enriched for chlorophyll, tetrapyrrole, and porphyrin biosynthesis, anthocyanin metabolism, and chloroplast protein import indicating coordinated regulation of plastid formation and photosynthetic apparatus assembly.

ego_ward_w9 <- enrichGO(
  gene          = genes_ward_w9_clean,
  OrgDb         = org.At.tair.db,
  keyType       = "TAIR",
  ont           = "BP",               # Biological Process
  pAdjustMethod = "BH",
  pvalueCutoff  = 0.1,
  qvalueCutoff  = 0.2
)

if (!is.null(ego_ward_w9) && nrow(as.data.frame(ego_ward_w9)) > 0) {
  print(
    dotplot(ego_ward_w9, showCategory = 15, label_format = 30) +
      ggtitle("GO Enrichment (BP): Ward Cluster 9") +
      theme_bw() +
      theme(
        axis.text.y = element_text(size = 9, color = "black"),
        plot.title  = element_text(face = "bold", size = 12)
      )
  )
  write.csv(as.data.frame(ego_ward_w9), "GO_Enrichment_Ward_Cluster9.csv", row.names = FALSE)
} else {
  cat("-> Tidak ada istilah GO yang signifikan pada Ward W9.\n")
}

For cluster 9, the GO terms are exactly the same as those for the K-Means cluster above, though the GeneRatio is slightly higher due to the smaller cluster size. The biological narrative remains the same.

ego_combined <- enrichGO(
  gene          = combined_genes_clean,
  OrgDb         = org.At.tair.db,
  keyType       = "TAIR",
  ont           = "BP",               # Biological Process
  pAdjustMethod = "BH",
  pvalueCutoff  = 0.1,
  qvalueCutoff  = 0.2
)

if (!is.null(ego_combined) && nrow(as.data.frame(ego_combined)) > 0) {
  print(
    dotplot(ego_combined, showCategory = 15, label_format = 30) +
      ggtitle("GO Enrichment (BP): Combined Genes (K6 + W9)") +
      theme_bw() +
      theme(
        axis.text.y = element_text(size = 9, color = "black"),
        plot.title  = element_text(face = "bold", size = 12)
      )
  )
  write.csv(as.data.frame(ego_combined), "GO_Enrichment_Combined_K6_W9.csv", row.names = FALSE)
} else {
  cat("-> Tidak ada istilah GO yang signifikan pada gen gabungan.\n")
}

The combined analysis yields identical GO terms to both individual clusters, confirming K-means and Ward capture the same gene module rather than a method-specific artifact. Note, however, that each term shows Count = 1 indicating the significant p.adjust reflects small gene-set size, not multiple converging genes.

4. Conclusion

All three plots converge on one theme: chloroplast biogenesis and pigment pathways aligning with the dry seed–to–seed development transition. This convergence across methods strengthens the finding, though the low gene count per term warrants caution.

5. Reference

Ashburner, M., Ball, C. A., Blake, J. A., et al. (2000). Gene ontology: tool for the unification of biology. Nature Genetics, 25(1), 25–29.

Benjamini, Y., & Hochberg, Y. (1995). Controlling the false discovery rate: a practical and powerful approach to multiple testing. Journal of the Royal Statistical Society: Series B, 57(1), 289–300.

Dunn, O. J. (1961). Multiple comparisons among means. Journal of the American Statistical Association, 56(293), 52–64.

Jolliffe, I. T., & Cadima, J. (2016). Principal component analysis: a review and recent developments. Philosophical Transactions of the Royal Society A, 374(2065), 20150202.

MacQueen, J. (1967). Some methods for classification and analysis of multivariate observations. Proceedings of the Fifth Berkeley Symposium on Mathematical Statistics and Probability, 1, 281–297.

Narsai, R., Law, S. R., Carrie, C., Xu, L., & Whelan, J. (2017). Extensive transcriptomic and epigenomic remodelling occurs during Arabidopsis thaliana germination. Genome Biology, 18, 172.

Papatheodorou, I., Moreno, P., Manning, J., et al. (2020). Expression Atlas update: from tissues to single cells. Nucleic Acids Research, 48(D1), D77–D83.

Ritchie, M. E., Phipson, B., Wu, D., Hu, Y., Law, C. W., Shi, W., & Smyth, G. K. (2015). limma powers differential expression analyses for RNA-sequencing and microarray studies. Nucleic Acids Research, 43(7), e47.

Ward, J. H. (1963). Hierarchical grouping to optimize an objective function. Journal of the American Statistical Association, 58(301), 236–244.

Wu, T., Hu, E., Xu, S., et al. (2021). clusterProfiler 4.0: A universal enrichment tool for interpreting omics data. The Innovation, 2(3), 100141.