This report analyzes the E-MTAB-9838 Arabidopsis thaliana RNA-seq data: untreated controls (no), flg22 treatment (flg22), and AtPep1 treatment (pep1), with three samples per treatment. The plots are shown inside this HTML report, so they will also be visible when the HTML report is published to RPubs.

Before knitting: Extract all files from the accompanying ZIP archive into one folder. Open this .Rmd file from that extracted folder in RStudio and click Knit. Both TSV input files must be in the same folder as this .Rmd file. Do not run setwd() inside the report.

1 Package installation (run in RStudio Console only)

The chunk below is not run during knitting. Remove the leading # and run the commands in the RStudio Console if packages are missing.

# install.packages("BiocManager")
# BiocManager::install(c("DESeq2", "EnhancedVolcano", "org.At.tair.db",
#                       "clusterProfiler", "enrichplot"),
#                      update = FALSE, ask = FALSE)
# install.packages(c("ggplot2", "knitr", "rmarkdown"))

2 Load packages and configure the report

3 Step 1: Load the raw count data

The raw count table contains two gene-information columns followed by nine sample count columns. Gene.ID contains Arabidopsis TAIR locus identifiers.

counts <- read.delim(counts_file, check.names = TRUE)

# R changes the original column names "Gene ID" and "Gene Name"
# to "Gene.ID" and "Gene.Name" when check.names = TRUE.
rownames(counts) <- counts$Gene.ID

# Keep gene names in a separate table.
genes <- counts[, c("Gene.ID", "Gene.Name"), drop = FALSE]

# Keep all nine sample columns as a matrix.
counts2 <- as.matrix(counts[, !(names(counts) %in% c("Gene.ID", "Gene.Name"))])

# Check that counts are unique, numeric, nonmissing, nonnegative integers.
stopifnot(!anyDuplicated(rownames(counts2)))
stopifnot(!anyDuplicated(colnames(counts2)))
stopifnot(is.numeric(counts2), !anyNA(counts2))
stopifnot(all(counts2 >= 0))
stopifnot(all(counts2 == floor(counts2)))
storage.mode(counts2) <- "integer"

cat("Genes:", nrow(counts2), "\n")
## Genes: 32833
cat("Samples:", ncol(counts2), "\n")
## Samples: 9
knitr::kable(head(genes, 5), caption = "Examples of Arabidopsis gene identifiers")
Examples of Arabidopsis gene identifiers
Gene.ID Gene.Name
AT1G01010 AT1G01010 NAC001
AT1G01020 AT1G01020 ARV1
AT1G01030 AT1G01030 NGA3
AT1G01040 AT1G01040 DCL1
AT1G01046 AT1G01046 MIR838A

4 Step 2: Load the metadata

metadata <- read.delim(metadata_file, check.names = TRUE)
rownames(metadata) <- metadata$Run

# Use only the stimulus variable.
metadata <- metadata[, "Sample.Characteristic.stimulus.", drop = FALSE]
colnames(metadata) <- "treatment"

# Rename the treatments consistently; no unknown values are silently renamed.
metadata$treatment[metadata$treatment == "none"] <- "no"
metadata$treatment[metadata$treatment == "AtPep1"] <- "pep1"
stopifnot(all(metadata$treatment %in% c("no", "flg22", "pep1")))

knitr::kable(metadata, caption = "Sample treatment metadata")
Sample treatment metadata
treatment
ERR4901839 no
ERR4901840 no
ERR4901841 no
ERR4901842 flg22
ERR4901843 flg22
ERR4901844 flg22
ERR4901845 pep1
ERR4901846 pep1
ERR4901847 pep1

5 Step 3: Match the metadata to the count matrix

# Every count-matrix sample must exist in the metadata.
stopifnot(all(colnames(counts2) %in% rownames(metadata)))

# Arrange the metadata rows in exactly the same order as the count columns.
metadata <- metadata[colnames(counts2), , drop = FALSE]
stopifnot(identical(colnames(counts2), rownames(metadata)))

cat("Sample alignment correct:",
    identical(colnames(counts2), rownames(metadata)), "\n")
## Sample alignment correct: TRUE

6 Step 4: Set the control as the reference group

metadata$treatment <- factor(metadata$treatment,
                             levels = c("no", "flg22", "pep1"))

# There should be three samples in each treatment group.
print(table(metadata$treatment))
## 
##    no flg22  pep1 
##     3     3     3

7 Step 5: Create and filter the DESeq2 dataset

dds <- DESeqDataSetFromMatrix(
  countData = counts2,
  colData = metadata,
  design = ~ treatment
)

# Exclude genes with 10 or fewer total counts across the nine samples.
dds <- dds[rowSums(counts(dds)) > 10, ]
cat("Genes retained after filtering:", nrow(dds), "\n")
## Genes retained after filtering: 21786

8 Step 6: Run differential expression analysis

dds <- DESeq(dds)
print(resultsNames(dds))
## [1] "Intercept"             "treatment_flg22_vs_no" "treatment_pep1_vs_no"

9 Step 7: Obtain the three treatment comparisons

Positive log2 fold change means higher expression in the first treatment named in the contrast.

res_flg22 <- results(dds, contrast = c("treatment", "flg22", "no"), alpha = 0.05)
res_pep1 <- results(dds, contrast = c("treatment", "pep1", "no"), alpha = 0.05)
res_pep1_vs_flg22 <- results(dds, contrast = c("treatment", "pep1", "flg22"),
                             alpha = 0.05)

cat("flg22 versus control:\n")
## flg22 versus control:
print(summary(res_flg22, alpha = 0.05))
## 
## out of 21786 with nonzero total read count
## adjusted p-value < 0.05
## LFC > 0 (up)       : 2785, 13%
## LFC < 0 (down)     : 1205, 5.5%
## outliers [1]       : 0, 0%
## low counts [2]     : 423, 1.9%
## (mean count < 2)
## [1] see 'cooksCutoff' argument of ?results
## [2] see 'independentFiltering' argument of ?results
## 
## NULL
cat("\npep1 versus control:\n")
## 
## pep1 versus control:
print(summary(res_pep1, alpha = 0.05))
## 
## out of 21786 with nonzero total read count
## adjusted p-value < 0.05
## LFC > 0 (up)       : 3075, 14%
## LFC < 0 (down)     : 2131, 9.8%
## outliers [1]       : 0, 0%
## low counts [2]     : 423, 1.9%
## (mean count < 2)
## [1] see 'cooksCutoff' argument of ?results
## [2] see 'independentFiltering' argument of ?results
## 
## NULL
cat("\npep1 versus flg22:\n")
## 
## pep1 versus flg22:
print(summary(res_pep1_vs_flg22, alpha = 0.05))
## 
## out of 21786 with nonzero total read count
## adjusted p-value < 0.05
## LFC > 0 (up)       : 1818, 8.3%
## LFC < 0 (down)     : 2170, 10%
## outliers [1]       : 0, 0%
## low counts [2]     : 2957, 14%
## (mean count < 10)
## [1] see 'cooksCutoff' argument of ?results
## [2] see 'independentFiltering' argument of ?results
## 
## NULL

10 Step 8: Export the three raw DESeq2 result tables (optional)

CSV files are saved to Arabidopsis_DESeq2_results when that folder is writable. Saving CSV files is separate from showing graphs in HTML.

save_csv(data.frame(Gene.ID = rownames(res_flg22), as.data.frame(res_flg22)),
         "flg22_vs_none.csv")
save_csv(data.frame(Gene.ID = rownames(res_pep1), as.data.frame(res_pep1)),
         "pep1_vs_none.csv")
save_csv(data.frame(Gene.ID = rownames(res_pep1_vs_flg22),
                    as.data.frame(res_pep1_vs_flg22)),
         "pep1_vs_flg22.csv")

11 Step 9: Choose the comparison for the remaining analyses

Change only active_comparison to one of the three names shown below, then knit the report again.

active_comparison <- "flg22_vs_none"
# Other choices: "pep1_vs_none" or "pep1_vs_flg22"

results_list <- list(
  flg22_vs_none = res_flg22,
  pep1_vs_none = res_pep1,
  pep1_vs_flg22 = res_pep1_vs_flg22
)

stopifnot(active_comparison %in% names(results_list))
result <- results_list[[active_comparison]]
cat("Selected comparison:", active_comparison, "\n")
## Selected comparison: flg22_vs_none

12 Step 10: Add Arabidopsis gene names to the selected results

No extra TAIR column is needed. The Gene.ID column already contains TAIR locus identifiers.

result_df <- as.data.frame(result)
result_df$Gene.ID <- rownames(result_df)
result_df$Gene.Name <- genes$Gene.Name[match(result_df$Gene.ID, genes$Gene.ID)]

result_df <- result_df[, c("Gene.ID", "Gene.Name", "baseMean",
                           "log2FoldChange", "lfcSE", "stat", "pvalue", "padj")]

save_csv(result_df, paste0("annotated_", active_comparison, ".csv"))
knitr::kable(head(result_df, 10), digits = 3,
             caption = "Preview of the selected DESeq2 results")
Preview of the selected DESeq2 results
Gene.ID Gene.Name baseMean log2FoldChange lfcSE stat pvalue padj
AT1G01010 AT1G01010 NAC001 103.723 -0.230 0.186 -1.231 0.218 0.470
AT1G01020 AT1G01020 ARV1 125.826 -0.295 0.170 -1.734 0.083 0.258
AT1G01030 AT1G01030 NGA3 20.010 -0.184 0.395 -0.465 0.642 0.831
AT1G01040 AT1G01040 DCL1 351.590 -0.205 0.128 -1.606 0.108 0.306
AT1G01046 AT1G01046 MIR838A 2.946 -0.474 0.984 -0.482 0.630 0.825
AT1G01050 AT1G01050 PPA1 828.007 0.037 0.109 0.335 0.737 0.882
AT1G01060 AT1G01060 LHY 881.313 0.284 0.181 1.568 0.117 0.322
AT1G01070 AT1G01070 40.907 -0.075 0.264 -0.286 0.775 0.903
AT1G01080 AT1G01080 949.753 0.174 0.089 1.960 0.050 0.181
AT1G01090 AT1G01090 PDH-E1 ALPHA 1891.917 -0.026 0.090 -0.285 0.776 0.904

13 Step 11: Examine selected plant immunity genes

genes_to_check <- c("AT3G45640",  # MPK3
                    "AT2G43790",  # MPK6
                    "AT5G46330",  # FLS2
                    "AT1G73080",  # PEPR1
                    "AT2G38470")  # WRKY33

immune_genes <- result_df[result_df$Gene.ID %in% genes_to_check, ]
knitr::kable(immune_genes, digits = 3,
             caption = "Differential expression of selected immune-response genes")
Differential expression of selected immune-response genes
Gene.ID Gene.Name baseMean log2FoldChange lfcSE stat pvalue padj
AT1G73080 AT1G73080 PEPR1 1617.871 2.301 0.100 23.001 0.000 0.000
AT2G38470 AT2G38470 WRKY33 10129.680 3.779 0.164 23.085 0.000 0.000
AT2G43790 AT2G43790 MPK6 476.651 0.074 0.101 0.731 0.465 0.709
AT3G45640 AT3G45640 MPK3 7652.691 3.376 0.148 22.837 0.000 0.000
AT5G46330 AT5G46330 FLS2 583.786 1.195 0.101 11.778 0.000 0.000

13.1 Normalized MPK3 counts across treatments

Important: This function produces a plot directly in the knitted HTML.

plotCounts(dds, gene = "AT3G45640", intgroup = "treatment")

14 Step 12: MA plot — shown in HTML and RPubs

The key fix is to call plotMA() directly. Do not wrap it in pdf() / dev.off(), which would send the figure to a file rather than to the HTML.

plotMA(result, main = active_comparison, ylim = c(-5, 5))

15 Step 13: Volcano plot — shown in HTML and RPubs

EnhancedVolcano() creates a ggplot object. Calling print(volcano_plot) displays it in the knitted HTML. Genes without symbols use Gene.ID as labels.

volcano_df <- result_df[is.finite(result_df$log2FoldChange), ]

plot_labels <- volcano_df$Gene.Name
missing_labels <- is.na(plot_labels) | trimws(plot_labels) == ""
plot_labels[missing_labels] <- volcano_df$Gene.ID[missing_labels]

# Missing adjusted p-values are shown as nonsignificant in the visualization.
volcano_df$plot_padj <- volcano_df$padj
volcano_df$plot_padj[is.na(volcano_df$plot_padj)] <- 1
volcano_df$plot_padj <- pmax(volcano_df$plot_padj, .Machine$double.xmin)

volcano_plot <- EnhancedVolcano(
  volcano_df,
  lab = plot_labels,
  x = "log2FoldChange",
  y = "plot_padj",
  pCutoff = 0.05,
  FCcutoff = 1,
  title = active_comparison
)
## Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
## ℹ Please use `linewidth` instead.
## ℹ The deprecated feature was likely used in the EnhancedVolcano package.
##   Please report the issue to the authors.
## This warning is displayed once per session.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.
## Warning: The `size` argument of `element_line()` is deprecated as of ggplot2 3.4.0.
## ℹ Please use the `linewidth` argument instead.
## ℹ The deprecated feature was likely used in the EnhancedVolcano package.
##   Please report the issue to the authors.
## This warning is displayed once per session.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.
print(volcano_plot)

16 Step 14: PCA plot — shown in HTML and RPubs

PCA includes all nine samples, irrespective of the comparison selected in Step 9.

vsd <- vst(dds, blind = FALSE)
pca_plot <- plotPCA(vsd, intgroup = "treatment")
print(pca_plot)

17 Step 15: Select significant genes and display the top 20

Use a data frame (sig_df) for tables and a vector of gene IDs (sig_genes) for enrichment analysis.

sig_rows <- !is.na(result_df$padj) & result_df$padj < 0.05
sig_df <- result_df[sig_rows, , drop = FALSE]
sig_df <- sig_df[order(sig_df$padj), , drop = FALSE]
sig_df <- sig_df[!duplicated(sig_df$Gene.ID), , drop = FALSE]

top20 <- head(sig_df[, c("Gene.ID", "Gene.Name", "log2FoldChange", "padj")], 20)
knitr::kable(top20, digits = 4,
             caption = "Top 20 significant genes by adjusted p-value")
Top 20 significant genes by adjusted p-value
Gene.ID Gene.Name log2FoldChange padj
AT1G61340 AT1G61340 5.1209 0
AT1G65390 AT1G65390 PP2A5 5.0324 0
AT1G80840 AT1G80840 WRKY40 4.6888 0
AT3G21070 AT3G21070 NADK1 3.4850 0
AT3G57530 AT3G57530 CPK32 3.5035 0
AT4G29780 AT4G29780 5.4790 0
AT2G40000 AT2G40000 HSPRO2 3.9985 0
AT1G56250 AT1G56250 VBF 6.9428 0
AT5G15870 AT5G15870 3.7547 0
AT2G26190 AT2G26190 IQM4 3.4541 0
AT4G01250 AT4G01250 WRKY22 4.9642 0
AT4G24570 AT4G24570 PUMP4 5.3029 0
AT3G02840 AT3G02840 7.3515 0
AT1G09940 AT1G09940 HEMA2 3.3827 0
AT1G72520 AT1G72520 LOX4 5.1681 0
AT1G21110 AT1G21110 IGMT3 4.2942 0
AT3G18710 AT3G18710 PUB29 4.4350 0
AT4G31550 AT4G31550 WRKY11 3.9659 0
AT1G23710 AT1G23710 4.9187 0
AT1G21120 AT1G21120 5.8244 0 
sig_genes <- unique(sig_df$Gene.ID)
background <- unique(result_df$Gene.ID[!is.na(result_df$padj)])
cat("Significant genes in", active_comparison, ":", length(sig_genes), "\n")
## Significant genes in flg22_vs_none : 3990

18 Step 16: GO Biological Process enrichment — shown in HTML and RPubs

A few TAIR IDs may fail annotation mapping; that is typically a warning, not a reason to stop the report. If no enriched terms are returned, the report displays a message rather than an empty plot.

ego <- NULL

if (length(sig_genes) > 0) {
  ego <- tryCatch(
    enrichGO(
      gene = sig_genes,
      universe = background,
      OrgDb = org.At.tair.db,
      keyType = "TAIR",
      ont = "BP",
      pAdjustMethod = "BH",
      pvalueCutoff = 0.05,
      qvalueCutoff = 0.05,
      readable = TRUE
    ),
    error = function(e) {
      cat("GO enrichment unavailable:", conditionMessage(e), "\n")
      NULL
    }
  )
}
## Warning in bitr(gene, fromType = fromType, toType = "ENTREZID", OrgDb = OrgDb):
## 1.8% of input gene IDs are fail to map...
## Warning in bitr(gene, fromType = fromType, toType = "ENTREZID", OrgDb = OrgDb):
## 3.74% of input gene IDs are fail to map...
if (!is.null(ego) && nrow(as.data.frame(ego)) > 0) {
  save_csv(as.data.frame(ego), paste0("GO_BP_", active_comparison, ".csv"))
  go_plot <- enrichplot::dotplot(ego, showCategory = 20)
  print(go_plot)  # Required to show the graph inside the knitted HTML.
} else {
  cat("No significant GO Biological Process terms were found.\n")
}

19 Step 17: KEGG enrichment — shown in HTML and RPubs

KEGG access may require an internet connection while knitting. If KEGG is unavailable, this section prints an explanation and the rest of the report still completes.

ekegg <- NULL

if (length(sig_genes) > 0) {
  ekegg <- tryCatch(
    enrichKEGG(
      gene = sig_genes,
      universe = background,
      organism = "ath",
      keyType = "kegg",
      pAdjustMethod = "BH",
      pvalueCutoff = 0.05
    ),
    error = function(e) {
      cat("KEGG enrichment unavailable:", conditionMessage(e), "\n")
      NULL
    }
  )
}

if (!is.null(ekegg) && nrow(as.data.frame(ekegg)) > 0) {
  save_csv(as.data.frame(ekegg), paste0("KEGG_", active_comparison, ".csv"))
  kegg_plot <- barplot(ekegg, showCategory = 20)
  print(kegg_plot)  # Required to show the graph inside the knitted HTML.
} else {
  cat("No significant KEGG pathways were found, or the service was unavailable.\n")
}

20 Step 18: Create a ranked gene vector for possible GSEA

geneList <- result$stat
names(geneList) <- rownames(result)
geneList <- geneList[is.finite(geneList)]
geneList <- sort(geneList, decreasing = TRUE)
cat("Genes in the ranked list:", length(geneList), "\n")
## Genes in the ranked list: 21786

21 Step 19: Export a final annotated differential expression table

The Direction column names the treatment in the numerator of the selected comparison. Not_Significant includes genes with missing adjusted p-values. No fold-change cutoff is applied in this classification.

final_report <- result_df[, c("Gene.ID", "Gene.Name", "baseMean",
                              "log2FoldChange", "lfcSE", "stat",
                              "pvalue", "padj")]

sig <- !is.na(final_report$padj) & final_report$padj < 0.05
treatment <- strsplit(active_comparison, "_vs_")[[1]][1]

final_report$Direction <- ifelse(
  sig & !is.na(final_report$log2FoldChange) & final_report$log2FoldChange > 0,
  paste0("Up_in_", treatment),
  ifelse(
    sig & !is.na(final_report$log2FoldChange) & final_report$log2FoldChange < 0,
    paste0("Down_in_", treatment),
    "Not_Significant"
  )
)

file_name <- paste0("DESeq2_", active_comparison, "_results.csv")
save_csv(final_report, file_name)

cat("Direction counts:\n")
## Direction counts:
print(table(final_report$Direction))
## 
##   Down_in_flg22 Not_Significant     Up_in_flg22 
##            1205           17796            2785
knitr::kable(head(final_report, 10), digits = 5,
             caption = "Preview of the final differential expression table")
Preview of the final differential expression table
Gene.ID Gene.Name baseMean log2FoldChange lfcSE stat pvalue padj Direction
AT1G01010 AT1G01010 NAC001 103.72261 -0.22961 0.18647 -1.23137 0.21819 0.46988 Not_Significant
AT1G01020 AT1G01020 ARV1 125.82629 -0.29510 0.17023 -1.73355 0.08300 0.25784 Not_Significant
AT1G01030 AT1G01030 NGA3 20.00956 -0.18391 0.39533 -0.46522 0.64177 0.83141 Not_Significant
AT1G01040 AT1G01040 DCL1 351.58964 -0.20521 0.12775 -1.60642 0.10818 0.30594 Not_Significant
AT1G01046 AT1G01046 MIR838A 2.94594 -0.47444 0.98384 -0.48223 0.62964 0.82506 Not_Significant
AT1G01050 AT1G01050 PPA1 828.00725 0.03659 0.10909 0.33546 0.73727 0.88235 Not_Significant
AT1G01060 AT1G01060 LHY 881.31253 0.28402 0.18114 1.56794 0.11689 0.32156 Not_Significant
AT1G01070 AT1G01070 40.90688 -0.07545 0.26369 -0.28611 0.77479 0.90339 Not_Significant
AT1G01080 AT1G01080 949.75274 0.17419 0.08887 1.96004 0.04999 0.18101 Not_Significant
AT1G01090 AT1G01090 PDH-E1 ALPHA 1891.91692 -0.02564 0.09009 -0.28457 0.77597 0.90399 Not_Significant

22 Finished

All figures above are rendered directly in the knitted HTML and should be visible in RPubs after you use RStudio’s Publish option on the HTML report. RPubs displays the report and embedded plots; it does not host the separately exported CSV files from your computer.

cat("Selected comparison:", active_comparison, "\n")
## Selected comparison: flg22_vs_none
cat("HTML report includes: MPK3 counts, MA plot, volcano plot, PCA plot,\n")
## HTML report includes: MPK3 counts, MA plot, volcano plot, PCA plot,
cat("and GO/KEGG plots when enriched terms are available.\n")
## and GO/KEGG plots when enriched terms are available.
cat("Local CSV output folder:", output_dir, "\n")
## Local CSV output folder: Arabidopsis_DESeq2_results
cat("Output folder writable:", can_save, "\n")
## Output folder writable: TRUE