0.1 Required packages

Install the following packages once in the RStudio Console if they are missing.

if (!requireNamespace("BiocManager", quietly = TRUE))
  install.packages("BiocManager")

BiocManager::install(c("DESeq2", "EnhancedVolcano"),
                     update = FALSE, ask = FALSE)

install.packages(c("ggplot2", "knitr"))

1 Problem 1: Load the count matrix correctly

Objective: Import the RNA-seq counts, retain gene identifiers and names in a separate annotation table, and create a sample-only count matrix.

counts <- read.csv("raw_counts.tsv", sep = "\t")

rownames(counts) <- counts$Gene.ID

genes = counts[, c("Gene.ID", "Gene.Name")]

counts = counts[, 3: ncol(counts)]

dim(counts)
## [1] 61860     6
head(rownames(counts), 3)
## [1] "ENSG00000000003" "ENSG00000000005" "ENSG00000000419"
head(colnames(counts), 3)
## [1] "ERR1736465" "ERR1736466" "ERR1736467"
head(genes, 3)
##                         Gene.ID Gene.Name
## ENSG00000000003 ENSG00000000003    TSPAN6
## ENSG00000000005 ENSG00000000005      TNMD
## ENSG00000000419 ENSG00000000419      DPM1

2 Problem 2: Load metadata and confirm sample matching

Objective: Match metadata rows to the exact sample order in the count matrix.

# Load metadata
metadata = read.csv("experimental_design.tsv", sep = "\t")

rownames(metadata) = metadata$Run

head(metadata)
##                   Run Sample.Characteristic.cell.line.
## ERR1736465 ERR1736465                           HS578T
## ERR1736466 ERR1736466                           HS578T
## ERR1736467 ERR1736467                           HS578T
## ERR1736468 ERR1736468                           HS578T
## ERR1736469 ERR1736469                           HS578T
## ERR1736470 ERR1736470                           HS578T
##            Sample.Characteristic.Ontology.Term.cell.line.
## ERR1736465           http://www.ebi.ac.uk/efo/EFO_0001192
## ERR1736466           http://www.ebi.ac.uk/efo/EFO_0001192
## ERR1736467           http://www.ebi.ac.uk/efo/EFO_0001192
## ERR1736468           http://www.ebi.ac.uk/efo/EFO_0001192
## ERR1736469           http://www.ebi.ac.uk/efo/EFO_0001192
## ERR1736470           http://www.ebi.ac.uk/efo/EFO_0001192
##            Sample.Characteristic.cell.type.
## ERR1736465                  epithelial cell
## ERR1736466                  epithelial cell
## ERR1736467                  epithelial cell
## ERR1736468                  epithelial cell
## ERR1736469                  epithelial cell
## ERR1736470                  epithelial cell
##            Sample.Characteristic.Ontology.Term.cell.type.
## ERR1736465      http://purl.obolibrary.org/obo/CL_0000066
## ERR1736466      http://purl.obolibrary.org/obo/CL_0000066
## ERR1736467      http://purl.obolibrary.org/obo/CL_0000066
## ERR1736468      http://purl.obolibrary.org/obo/CL_0000066
## ERR1736469      http://purl.obolibrary.org/obo/CL_0000066
## ERR1736470      http://purl.obolibrary.org/obo/CL_0000066
##            Sample.Characteristic.disease.
## ERR1736465               breast carcinoma
## ERR1736466               breast carcinoma
## ERR1736467               breast carcinoma
## ERR1736468               breast carcinoma
## ERR1736469               breast carcinoma
## ERR1736470               breast carcinoma
##            Sample.Characteristic.Ontology.Term.disease.
## ERR1736465         http://www.ebi.ac.uk/efo/EFO_0000305
## ERR1736466         http://www.ebi.ac.uk/efo/EFO_0000305
## ERR1736467         http://www.ebi.ac.uk/efo/EFO_0000305
## ERR1736468         http://www.ebi.ac.uk/efo/EFO_0000305
## ERR1736469         http://www.ebi.ac.uk/efo/EFO_0000305
## ERR1736470         http://www.ebi.ac.uk/efo/EFO_0000305
##            Sample.Characteristic.genetic.modification.
## ERR1736465                                        none
## ERR1736466                                        none
## ERR1736467                                        none
## ERR1736468                              gene knock out
## ERR1736469                              gene knock out
## ERR1736470                              gene knock out
##            Sample.Characteristic.Ontology.Term.genetic.modification.
## ERR1736465                      http://www.ebi.ac.uk/efo/EFO_0001461
## ERR1736466                      http://www.ebi.ac.uk/efo/EFO_0001461
## ERR1736467                      http://www.ebi.ac.uk/efo/EFO_0001461
## ERR1736468                      http://www.ebi.ac.uk/efo/EFO_0000506
## ERR1736469                      http://www.ebi.ac.uk/efo/EFO_0000506
## ERR1736470                      http://www.ebi.ac.uk/efo/EFO_0000506
##            Sample.Characteristic.genotype.
## ERR1736465              wild type genotype
## ERR1736466              wild type genotype
## ERR1736467              wild type genotype
## ERR1736468                  Snai1 knockout
## ERR1736469                  Snai1 knockout
## ERR1736470                  Snai1 knockout
##            Sample.Characteristic.Ontology.Term.genotype.
## ERR1736465                                            NA
## ERR1736466                                            NA
## ERR1736467                                            NA
## ERR1736468                                            NA
## ERR1736469                                            NA
## ERR1736470                                            NA
##            Sample.Characteristic.organism.
## ERR1736465                    Homo sapiens
## ERR1736466                    Homo sapiens
## ERR1736467                    Homo sapiens
## ERR1736468                    Homo sapiens
## ERR1736469                    Homo sapiens
## ERR1736470                    Homo sapiens
##            Sample.Characteristic.Ontology.Term.organism. Factor.Value.genotype.
## ERR1736465 http://purl.obolibrary.org/obo/NCBITaxon_9606     wild type genotype
## ERR1736466 http://purl.obolibrary.org/obo/NCBITaxon_9606     wild type genotype
## ERR1736467 http://purl.obolibrary.org/obo/NCBITaxon_9606     wild type genotype
## ERR1736468 http://purl.obolibrary.org/obo/NCBITaxon_9606         Snai1 knockout
## ERR1736469 http://purl.obolibrary.org/obo/NCBITaxon_9606         Snai1 knockout
## ERR1736470 http://purl.obolibrary.org/obo/NCBITaxon_9606         Snai1 knockout
##            Factor.Value.Ontology.Term.genotype. Analysed
## ERR1736465                                   NA      Yes
## ERR1736466                                   NA      Yes
## ERR1736467                                   NA      Yes
## ERR1736468                                   NA      Yes
## ERR1736469                                   NA      Yes
## ERR1736470                                   NA      Yes
# Only keep columns of interest
metadata = metadata[, c("Sample.Characteristic.genotype."), drop=FALSE]

# Rename column
colnames(metadata) = c("genotype")
metadata
##                      genotype
## ERR1736465 wild type genotype
## ERR1736466 wild type genotype
## ERR1736467 wild type genotype
## ERR1736468     Snai1 knockout
## ERR1736469     Snai1 knockout
## ERR1736470     Snai1 knockout
# Remove spaces in names to avoid DESeq warnings
metadata$genotype <- ifelse(metadata$genotype == "wild type genotype", "WT",
                            "KO")

# Make sure raw counts and metadata are in the same order / aligned
# Match metadata with the count matrix
x <- match(colnames(counts), rownames(metadata))

# Check that all samples were found
stopifnot(!anyNA(x))

# Reorder metadata
metadata <- metadata[x, , drop = FALSE]

# Confirm alignment
sample_alignment <- identical(colnames(counts), rownames(metadata))
cat("Sample alignment correct:", sample_alignment, "\n")
## Sample alignment correct: TRUE
stopifnot(sample_alignment)

3 Problem 3: Create Condition with WT as the reference

Objective: Create metadata$Condition and put WT first in the factor levels.

# Turn genotype into a factor
metadata$genotype = factor(metadata$genotype, levels=c("WT", "KO"))

metadata$genotype
## [1] WT WT WT KO KO KO
## Levels: WT KO

4 Problem 4: Build the DESeq2 dataset object

Objective: Combine the sample-aligned counts, metadata, and experimental design.

dds <- DESeqDataSetFromMatrix(
  countData = counts,
  colData = metadata,
  design = ~genotype
)
dds
## class: DESeqDataSet 
## dim: 61860 6 
## metadata(1): version
## assays(1): counts
## rownames(61860): ENSG00000000003 ENSG00000000005 ... ENSG00000290165
##   ENSG00000290166
## rowData names(0):
## colnames(6): ERR1736465 ERR1736466 ... ERR1736469 ERR1736470
## colData names(1): genotype

5 Problem 5: Filter lowly expressed genes

Objective: Remove genes with a total raw count of 10 or fewer across all samples.

# Ignore genes with low counts
dds <- dds[rowSums(counts(dds)) > 10, ]
dds
## class: DESeqDataSet 
## dim: 14548 6 
## metadata(1): version
## assays(1): counts
## rownames(14548): ENSG00000000003 ENSG00000000419 ... ENSG00000289554
##   ENSG00000289707
## rowData names(0):
## colnames(6): ERR1736465 ERR1736466 ... ERR1736469 ERR1736470
## colData names(1): genotype

6 Problem 6: Run DESeq2 and compare significance thresholds

Objective: Fit the DESeq2 model and compare results at alpha = 0.05 and alpha = 1e-5.

dds <- DESeq(dds)

res <- results(dds, contrast = c("genotype", "KO", "WT"), alpha = 0.05)

res_strict <- results(dds, contrast = c("genotype", "KO", "WT"), alpha = 1e-5)

cat("Results summary at alpha = 0.05:\n")
## Results summary at alpha = 0.05:
summary(res, alpha = 0.05)
## 
## out of 14548 with nonzero total read count
## adjusted p-value < 0.05
## LFC > 0 (up)       : 2155, 15%
## LFC < 0 (down)     : 2068, 14%
## outliers [1]       : 32, 0.22%
## low counts [2]     : 1411, 9.7%
## (mean count < 3)
## [1] see 'cooksCutoff' argument of ?results
## [2] see 'independentFiltering' argument of ?results
cat("\nResults summary at alpha = 1e-5:\n")
## 
## Results summary at alpha = 1e-5:
summary(res_strict, alpha = 1e-5)
## 
## out of 14548 with nonzero total read count
## adjusted p-value < 1e-05
## LFC > 0 (up)       : 769, 5.3%
## LFC < 0 (down)     : 538, 3.7%
## outliers [1]       : 32, 0.22%
## low counts [2]     : 2257, 16%
## (mean count < 5)
## [1] see 'cooksCutoff' argument of ?results
## [2] see 'independentFiltering' argument of ?results

7 Problem 7: Convert results to a data frame and merge gene names

Objective: Combine DESeq2 results with Gene.ID and Gene.Name.

result_df <- as.data.frame(res)

result_df$Gene.ID <- rownames(result_df)

result_df <- merge(result_df,genes,by = "Gene.ID")

head(result_df)
##           Gene.ID   baseMean log2FoldChange      lfcSE       stat       pvalue
## 1 ENSG00000000003  568.53155    -0.29035333 0.12614348 -2.3017704 2.134813e-02
## 2 ENSG00000000419 3024.93792     0.08912822 0.07892131  1.1293302 2.587586e-01
## 3 ENSG00000000457   65.74124     0.30418198 0.31272292  0.9726885 3.307081e-01
## 4 ENSG00000000460  404.36962    -0.25642207 0.14721473 -1.7418234 8.153935e-02
## 5 ENSG00000001036 2665.73118    -0.04814284 0.09701088 -0.4962623 6.197093e-01
## 6 ENSG00000001084 2025.65061    -0.60795968 0.12152635 -5.0026984 5.653336e-07
##           padj Gene.Name
## 1 6.278438e-02    TSPAN6
## 2 4.091495e-01      DPM1
## 3 4.886054e-01     SCYL3
## 4 1.725175e-01  C1orf112
## 5 7.491966e-01     FUCA2
## 6 6.028231e-06      GCLC

8 Problem 8: Find the twenty most significant genes

Objective: Remove missing adjusted p-values and rank genes by padj (smallest first).

results_valid_padj <- result_df[!is.na(result_df$padj), , drop = FALSE]
results_valid_padj <- results_valid_padj[order(results_valid_padj$padj), , drop = FALSE]

top20 <- head(results_valid_padj[, c("Gene.ID", "Gene.Name",
                                 "log2FoldChange", "padj")], 20)
knitr::kable(top20, digits = 2,
             caption = "Top 20 genes ranked by adjusted p-value")
Top 20 genes ranked by adjusted p-value
Gene.ID Gene.Name log2FoldChange padj
6885 ENSG00000151632 AKR1C2 -6.16 0
2676 ENSG00000108691 CCL2 8.49 0
1561 ENSG00000096696 DSP -4.30 0
4076 ENSG00000124212 PTGIS -5.39 0
12769 ENSG00000225383 SFTA1P 4.75 0
7848 ENSG00000163017 ACTG2 3.00 0
3119 ENSG00000113389 NPR3 5.66 0
11085 ENSG00000187134 AKR1C1 -2.93 0
1123 ENSG00000079482 OPHN1 -3.56 0
7953 ENSG00000163584 RPL22L1 -2.47 0
2850 ENSG00000110852 CLEC2B -5.05 0
11945 ENSG00000198959 TGM2 2.56 0
4991 ENSG00000133816 MICAL2 2.36 0
8946 ENSG00000169429 CXCL8 -4.88 0
8420 ENSG00000166165 CKB 3.63 0
13013 ENSG00000231290 APCDD1L-DT -1.65 0
8605 ENSG00000167355 OR51B5 -6.79 0
2525 ENSG00000106976 DNM1 5.49 0
3310 ENSG00000115380 EFEMP1 -1.93 0
7481 ENSG00000159167 STC1 3.14 0

9 Problem 9: Examine four genes of interest

Objective: Extract THY1, SFMBT2, PASD1, and SNAI1 and interpret KO-versus-WT changes.

genes_check = c("THY1", "SFMBT2", "PASD1", "SNAI1")
result_df[result_df$Gene.Name %in% genes_check, ]
##               Gene.ID   baseMean log2FoldChange     lfcSE      stat
## 4078  ENSG00000124216   53.88646      -4.155055 0.4547677 -9.136655
## 7087  ENSG00000154096 3146.01390      12.371934 0.8406168 14.717685
## 8402  ENSG00000166049  251.22854      11.075918 1.1921449  9.290748
## 11913 ENSG00000198879  260.22214     -11.807534 1.2088618 -9.767480
##             pvalue         padj Gene.Name
## 4078  6.441369e-20 3.488188e-18     SNAI1
## 7087  4.963618e-49 1.182695e-46      THY1
## 8402  1.532072e-20 9.167944e-19     PASD1
## 11913 1.552653e-22 1.105844e-20    SFMBT2

10 Problem 10: Make a volcano plot

Objective: Visualize effect size (log2FoldChange) and statistical significance (padj).

volcano_df <- result_df[
  !is.na(result_df$log2FoldChange) &
    !is.na(result_df$padj) & is.finite(result_df$padj),
  , drop = FALSE
]
volcano_df$Label <- ifelse(
  is.na(volcano_df$Gene.Name) | volcano_df$Gene.Name == "",
  volcano_df$Gene.ID, volcano_df$Gene.Name
)

EnhancedVolcano(
  volcano_df,
  lab = volcano_df$Label,
  x = "log2FoldChange",
  y = "padj",
  pCutoff = 0.05,
  FCcutoff = 1,
  title = "KO versus WT: Differential gene expression",
  subtitle = "DESeq2 results (adjusted p-value cutoff: 0.05)",
  pointSize = 1.6,
  labSize = 3.0,
  legendPosition = "right"
)
## 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.

11 Problem 11: Diagnose a sample-mismatch bug

Objective: Simulate incorrectly paired sample labels, fit DESeq2 again, and compare top hits.

#Can do, but will take time.

12 Problem 12: Compare three low-count filtering strategies

Objective: Run DESeq2 with total-count filters > 0, > 10, and > 50 and compare findings.

#Can do, but will take time.

Interpretation: More stringent count filtering retains fewer genes. Removing weakly informative genes can change the number of tests, DESeq2’s independent filtering, and adjusted p-values. A different significant-gene count is not, by itself, evidence that one filter is biologically preferable; the choice should be justified by coverage and study design.

13 Problem 13: Export a collaborator-ready results report

Objective: Create a clean CSV with the requested result columns plus a three-category expression-direction label.

final_report <- result_df[, c("Gene.ID", "Gene.Name", "baseMean",
                                  "log2FoldChange", "lfcSE", "stat",
                                  "pvalue", "padj")]
final_report$Direction <- ifelse(
  !is.na(final_report$padj) & final_report$padj < 0.05 &
    !is.na(final_report$log2FoldChange) &
    final_report$log2FoldChange > 0,
  "Up_in_KO",
  ifelse(!is.na(final_report$padj) & final_report$padj < 0.05 &
           !is.na(final_report$log2FoldChange) &
           final_report$log2FoldChange < 0,
         "Down_in_KO", "Not_Significant")
)

# Export beside the R Markdown document (the knitting working directory).
write.csv(final_report, "DESeq2_WT_vs_KO_results.csv", row.names = FALSE)

cat("Direction counts:\n")
## Direction counts:
print(table(final_report$Direction))
## 
##      Down_in_KO Not_Significant        Up_in_KO 
##            2068           10325            2155
knitr::kable(head(final_report, 10), digits = 5,
             caption = "Preview of the collaborator-ready CSV")
Preview of the collaborator-ready CSV
Gene.ID Gene.Name baseMean log2FoldChange lfcSE stat pvalue padj Direction
ENSG00000000003 TSPAN6 568.53155 -0.29035 0.12614 -2.30177 0.02135 0.06278 Not_Significant
ENSG00000000419 DPM1 3024.93792 0.08913 0.07892 1.12933 0.25876 0.40915 Not_Significant
ENSG00000000457 SCYL3 65.74124 0.30418 0.31272 0.97269 0.33071 0.48861 Not_Significant
ENSG00000000460 C1orf112 404.36962 -0.25642 0.14721 -1.74182 0.08154 0.17252 Not_Significant
ENSG00000001036 FUCA2 2665.73118 -0.04814 0.09701 -0.49626 0.61971 0.74920 Not_Significant
ENSG00000001084 GCLC 2025.65061 -0.60796 0.12153 -5.00270 0.00000 0.00001 Down_in_KO
ENSG00000001167 NFYA 1042.36797 -0.15674 0.12461 -1.25784 0.20845 0.34933 Not_Significant
ENSG00000001460 STPG1 118.52593 9.99199 1.20429 8.29701 0.00000 0.00000 Up_in_KO
ENSG00000001461 NIPAL3 518.97013 1.37899 0.14409 9.57042 0.00000 0.00000 Up_in_KO
ENSG00000001497 LAS1L 813.21547 0.25522 0.10188 2.50522 0.01224 0.03989 Up_in_KO 
cat("\nCSV exported as: DESeq2_WT_vs_KO_results.csv\n")
## 
## CSV exported as: DESeq2_WT_vs_KO_results.csv

Data source: Expression Atlas E-MTAB-5244.