1. Packages

library(DESeq2)

## Loading required package: S4Vectors

## Loading required package: stats4

## Loading required package: BiocGenerics

## Loading required package: generics

## Warning: package 'generics' was built under R version 4.6.1

## 
## Attaching package: 'generics'

## The following objects are masked from 'package:base':
## 
##     as.difftime, as.factor, as.ordered, intersect, is.element, setdiff,
##     setequal, union

## 
## Attaching package: 'BiocGenerics'

## The following objects are masked from 'package:stats':
## 
##     IQR, mad, sd, var, xtabs

## The following objects are masked from 'package:base':
## 
##     anyDuplicated, aperm, append, as.data.frame, basename, cbind,
##     colnames, dirname, do.call, duplicated, eval, evalq, Filter, Find,
##     get, grep, grepl, is.unsorted, lapply, Map, mapply, match, mget,
##     order, paste, pmax, pmax.int, pmin, pmin.int, Position, rank,
##     rbind, Reduce, rownames, sapply, saveRDS, table, tapply, unique,
##     unsplit, which.max, which.min

## 
## Attaching package: 'S4Vectors'

## The following object is masked from 'package:utils':
## 
##     findMatches

## The following objects are masked from 'package:base':
## 
##     expand.grid, I, unname

## Loading required package: IRanges

## 
## Attaching package: 'IRanges'

## The following object is masked from 'package:grDevices':
## 
##     windows

## Loading required package: GenomicRanges

## Loading required package: Seqinfo

## Loading required package: SummarizedExperiment

## Loading required package: MatrixGenerics

## Loading required package: matrixStats

## Warning: package 'matrixStats' was built under R version 4.6.1

## 
## Attaching package: 'MatrixGenerics'

## The following objects are masked from 'package:matrixStats':
## 
##     colAlls, colAnyNAs, colAnys, colAvgsPerRowSet, colCollapse,
##     colCounts, colCummaxs, colCummins, colCumprods, colCumsums,
##     colDiffs, colIQRDiffs, colIQRs, colLogSumExps, colMadDiffs,
##     colMads, colMaxs, colMeans2, colMedians, colMins, colOrderStats,
##     colProds, colQuantiles, colRanges, colRanks, colSdDiffs, colSds,
##     colSums2, colTabulates, colVarDiffs, colVars, colWeightedMads,
##     colWeightedMeans, colWeightedMedians, colWeightedSds,
##     colWeightedVars, rowAlls, rowAnyNAs, rowAnys, rowAvgsPerColSet,
##     rowCollapse, rowCounts, rowCummaxs, rowCummins, rowCumprods,
##     rowCumsums, rowDiffs, rowIQRDiffs, rowIQRs, rowLogSumExps,
##     rowMadDiffs, rowMads, rowMaxs, rowMeans2, rowMedians, rowMins,
##     rowOrderStats, rowProds, rowQuantiles, rowRanges, rowRanks,
##     rowSdDiffs, rowSds, rowSums2, rowTabulates, rowVarDiffs, rowVars,
##     rowWeightedMads, rowWeightedMeans, rowWeightedMedians,
##     rowWeightedSds, rowWeightedVars

## Loading required package: Biobase

## Welcome to Bioconductor
## 
##     Vignettes contain introductory material; view with
##     'browseVignettes()'. To cite Bioconductor, see
##     'citation("Biobase")', and for packages 'citation("pkgname")'.

## 
## Attaching package: 'Biobase'

## The following object is masked from 'package:MatrixGenerics':
## 
##     rowMedians

## The following objects are masked from 'package:matrixStats':
## 
##     anyMissing, rowMedians

library(apeglm)
library(dplyr)

## 
## Attaching package: 'dplyr'

## The following object is masked from 'package:Biobase':
## 
##     combine

## The following object is masked from 'package:matrixStats':
## 
##     count

## The following objects are masked from 'package:GenomicRanges':
## 
##     intersect, setdiff, union

## The following object is masked from 'package:Seqinfo':
## 
##     intersect

## The following objects are masked from 'package:IRanges':
## 
##     collapse, desc, intersect, setdiff, slice, union

## The following objects are masked from 'package:S4Vectors':
## 
##     first, intersect, rename, setdiff, setequal, union

## The following objects are masked from 'package:BiocGenerics':
## 
##     combine, intersect, setdiff, setequal, union

## The following object is masked from 'package:generics':
## 
##     explain

## The following objects are masked from 'package:stats':
## 
##     filter, lag

## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union

library(org.Mm.eg.db)

## Loading required package: AnnotationDbi

## 
## Attaching package: 'AnnotationDbi'

## The following object is masked from 'package:dplyr':
## 
##     select

##

library(AnnotationDbi)
library(clusterProfiler)

##

## clusterProfiler v4.20.0 Learn more at https://yulab-smu.top/contribution-knowledge-mining/
## 
## Please cite:
## 
## T Wu, E Hu, S Xu, M Chen, P Guo, Z Dai, T Feng, L Zhou, W Tang, L Zhan,
## X Fu, S Liu, X Bo, and G Yu. clusterProfiler 4.0: A universal
## enrichment tool for interpreting omics data. The Innovation. 2021,
## 2(3):100141

## 
## Attaching package: 'clusterProfiler'

## The following object is masked from 'package:AnnotationDbi':
## 
##     select

## The following object is masked from 'package:IRanges':
## 
##     slice

## The following object is masked from 'package:S4Vectors':
## 
##     rename

## The following object is masked from 'package:stats':
## 
##     filter

library(pathview)

##

## ##############################################################################
## Pathview is an open source software package distributed under GNU General
## Public License version 3 (GPLv3). Details of GPLv3 is available at
## http://www.gnu.org/licenses/gpl-3.0.html. Particullary, users are required to
## formally cite the original Pathview paper (not just mention it) in publications
## or products. For details, do citation("pathview") within R.
## 
## The pathview downloads and uses KEGG data. Non-academic uses may require a KEGG
## license agreement (details at http://www.kegg.jp/kegg/legal.html).
## ##############################################################################

library(enrichplot)

## enrichplot v1.32.0 Learn more at https://yulab-smu.top/contribution-knowledge-mining/
## 
## Please cite:
## 
## Guangchuang Yu, Li-Gen Wang, and Qing-Yu He. ChIPseeker: an
## R/Bioconductor package for ChIP peak annotation, comparison and
## visualization. Bioinformatics. 2015, 31(14):2382-2383

2. Import count matrix and sample metadata

countData <- read.csv("GLDS-102_rna_seq_Unnormalized_Counts.csv", row.names = 1)
colData   <- read.csv("GLDS-102_rna_seq_SampleTable1.csv", sep = ",", row.names = 1)

3. Prepare metadata (experimental design)

colData$group <- factor(colData$group, levels = c("0", "1"))

4. Quality check and sample matching

stopifnot(all(rownames(colData) %in% colnames(countData)))
stopifnot(all(colnames(countData) %in% rownames(colData)))

countData <- countData[, rownames(colData)]

identical(colnames(countData), rownames(colData))

## [1] TRUE

5. Prepare count matrix

countData <- countData %>%
  mutate(across(where(is.numeric), as.integer))


rownames(countData) <- sub("\\..*$", "", rownames(countData))

6. Create DESeq2 dataset

dds <- DESeqDataSetFromMatrix(
  countData = countData,
  colData   = colData,
  design    = ~ group
)

7. Filter low expression genes

min_group_size <- min(table(colData$group))

keep <- rowSums(counts(dds) >= 10) >= min_group_size
dds  <- dds[keep, ]

8. Normalization and statistical modeling

dds <- DESeq(dds)

## estimating size factors

## estimating dispersions

## gene-wise dispersion estimates

## mean-dispersion relationship

## final dispersion estimates

## fitting model and testing

## -- replacing outliers and refitting for 9 genes
## -- DESeq argument 'minReplicatesForReplace' = 7 
## -- original counts are preserved in counts(dds)

## estimating dispersions

## fitting model and testing

plotDispEsts(dds)

9. Extract differential expression results

resultsNames(dds)

## [1] "Intercept"    "group_1_vs_0"

res <- results(dds, contrast = c("group", "1", "0"), alpha = 0.05)
summary(res)

## 
## out of 18165 with nonzero total read count
## adjusted p-value < 0.05
## LFC > 0 (up)       : 79, 0.43%
## LFC < 0 (down)     : 143, 0.79%
## outliers [1]       : 9, 0.05%
## low counts [2]     : 2817, 16%
## (mean count < 29)
## [1] see 'cooksCutoff' argument of ?results
## [2] see 'independentFiltering' argument of ?results

res_shrunk <- lfcShrink(dds, coef = "group_1_vs_0", type = "apeglm")

## using 'apeglm' for LFC shrinkage. If used in published research, please cite:
##     Zhu, A., Ibrahim, J.G., Love, M.I. (2018) Heavy-tailed prior distributions for
##     sequence count data: removing the noise and preserving large differences.
##     Bioinformatics. https://doi.org/10.1093/bioinformatics/bty895

## Warning in nbinomGLM(x = x, Y = YNZ, size = size, weights = weightsNZ, offset =
## offsetNZ, : the line search routine failed, unable to sufficiently decrease the
## function value

res_df <- as.data.frame(res)
write.csv(res_df, file = "Mouse_DESeq2_Results.csv", row.names = TRUE)

10. Visualization

# PCA Plot
vsd <- vst(dds, blind = FALSE)
pdf("PCA_Plot.pdf")
plotPCA(vsd, intgroup = "group")

## using ntop=500 top features by variance

dev.off()

## png 
##   2

# heatmap 
sampleDists <- dist(t(assay(vsd)))
pdf("Sample_Distance_Heatmap.pdf")
heatmap(as.matrix(sampleDists))
dev.off()

## png 
##   2

# MA Plot 
pdf("MA_Plot.pdf")
plotMA(res_shrunk, ylim = c(-2, 2))
dev.off()

## png 
##   2

11. Functional annotation

res_df$SYMBOL <- mapIds(
  org.Mm.eg.db,
  keys      = rownames(res_df),
  column    = "SYMBOL",
  keytype   = "ENSEMBL",
  multiVals = "first"
)

## 'select()' returned 1:many mapping between keys and columns

res_df$ENTREZID <- mapIds(
  org.Mm.eg.db,
  keys      = rownames(res_df),
  column    = "ENTREZID",
  keytype   = "ENSEMBL",
  multiVals = "first"
)

## 'select()' returned 1:many mapping between keys and columns

res_df$GENENAME <- mapIds(
  org.Mm.eg.db,
  keys      = rownames(res_df),
  column    = "GENENAME",
  keytype   = "ENSEMBL",
  multiVals = "first"
)

## 'select()' returned 1:many mapping between keys and columns

res_df <- res_df[order(res_df$padj), ]
write.csv(res_df, "DESeq2_final_result.csv", row.names = TRUE)

12. Filter significant DEGs

DEGs <- res_df %>%
  filter(
    !is.na(padj),
    !is.na(log2FoldChange),
    padj < 0.05
  )

write.csv(DEGs, "DESeq2_DEGs_file.csv", row.names = TRUE)

13. Pathway analysis

geneList <- na.omit(DEGs$ENTREZID)
universe <- na.omit(res_df$ENTREZID)

# GO Enrichment analysis
ego_BP <- enrichGO(
  gene          = geneList,
  universe      = universe,
  OrgDb         = org.Mm.eg.db,
  keyType       = "ENTREZID",
  ont           = "BP",
  pAdjustMethod = "BH",
  pvalueCutoff  = 0.05,
  qvalueCutoff  = 0.05,
  readable      = TRUE
)


head(ego_BP)

##                    ID
## GO:0033490 GO:0033490
## GO:0033489 GO:0033489
## GO:0036196 GO:0036196
## GO:0036197 GO:0036197
## GO:0010142 GO:0010142
## GO:1902767 GO:1902767
##                                                              Description
## GO:0033490              cholesterol biosynthetic process via lathosterol
## GO:0033489              cholesterol biosynthetic process via desmosterol
## GO:0036196                                  zymosterol metabolic process
## GO:0036197                               zymosterol biosynthetic process
## GO:0010142 farnesyl diphosphate biosynthetic process, mevalonate pathway
## GO:1902767                isoprenoid biosynthetic process via mevalonate
##            GeneRatio  BgRatio RichFactor FoldEnrichment   zScore       pvalue
## GO:0033490    10/206 20/14390  0.5000000       34.92718 18.29717 4.737170e-14
## GO:0033489    10/206 21/14390  0.4761905       33.26399 17.83052 8.931676e-14
## GO:0036196     9/206 16/14390  0.5625000       39.29308 18.46890 2.223835e-13
## GO:0036197     9/206 16/14390  0.5625000       39.29308 18.46890 2.223835e-13
## GO:0010142     7/206 10/14390  0.7000000       48.89806 18.25946 1.287796e-11
## GO:1902767     7/206 10/14390  0.7000000       48.89806 18.25946 1.287796e-11
##                p.adjust       qvalue
## GO:0033490 1.244629e-10 4.432124e-11
## GO:0033489 1.244629e-10 4.432124e-11
## GO:0036196 1.549457e-10 5.517617e-11
## GO:0036197 1.549457e-10 5.517617e-11
## GO:0010142 5.981812e-09 2.130123e-09
## GO:1902767 5.981812e-09 2.130123e-09
##                                                          geneID Count
## GO:0033490 Mvd/Dhcr7/Idi1/Fdps/Acat2/Msmo1/Sqle/Mvk/Nsdhl/Hmgcr    10
## GO:0033489 Mvd/Dhcr7/Idi1/Fdps/Acat2/Msmo1/Sqle/Mvk/Nsdhl/Hmgcr    10
## GO:0036196       Mvd/Idi1/Fdps/Acat2/Msmo1/Sqle/Mvk/Nsdhl/Hmgcr     9
## GO:0036197       Mvd/Idi1/Fdps/Acat2/Msmo1/Sqle/Mvk/Nsdhl/Hmgcr     9
## GO:0010142                 Mvd/Idi1/Fdps/Acat2/Mvk/Hmgcs2/Hmgcr     7
## GO:1902767                 Mvd/Idi1/Fdps/Acat2/Mvk/Hmgcs2/Hmgcr     7

write.csv(as.data.frame(ego_BP), file = "egoBP.csv", row.names = FALSE)

dotplot(ego_BP, showCategory = 20)

barplot(ego_BP, showCategory = 20)

# KEGG Enrichment
ekegg <- enrichKEGG(
  gene         = geneList,
  universe     = universe,
  organism     = "mmu",
  pvalueCutoff = 0.05
)

## Reading KEGG annotation online: "https://rest.kegg.jp/link/mmu/pathway"...

## Reading KEGG annotation online: "https://rest.kegg.jp/list/pathway/mmu"...

head(ekegg)

##                    category                              subcategory       ID
## mmu04710 Organismal Systems                 Environmental adaptation mmu04710
## mmu00900         Metabolism Metabolism of terpenoids and polyketides mmu00900
## mmu00100         Metabolism                         Lipid metabolism mmu00100
##                              Description GeneRatio BgRatio RichFactor
## mmu04710                Circadian rhythm    10/135 34/6777  0.2941176
## mmu00900 Terpenoid backbone biosynthesis     8/135 23/6777  0.3478261
## mmu00100            Steroid biosynthesis     4/135 18/6777  0.2222222
##          FoldEnrichment    zScore       pvalue     p.adjust       qvalue
## mmu04710       14.76471 11.470543 6.148393e-10 1.463318e-07 1.463318e-07
## mmu00900       17.46087 11.273021 7.679842e-09 9.139012e-07 9.139012e-07
## mmu00100       11.15556  6.150393 3.709734e-04 2.943056e-02 2.943056e-02
##                                                                 geneID Count
## mmu04710 13170/12953/353187/11865/18143/79362/18628/18627/18030/217166    10
## mmu00900           56075/192156/319554/110196/110460/17855/15360/15357     8
## mmu00100                                       13360/66234/20775/18194     4

write.csv(as.data.frame(ekegg), "KEGG_Pathway.csv", row.names = FALSE)
dotplot(ekegg, showCategory = 20)

foldChanges <- res_df$log2FoldChange
names(foldChanges) <- res_df$ENTREZID
foldChanges <- foldChanges[!is.na(names(foldChanges))]


top_pathways <- head(ekegg$ID, 20)

pathview_dir <- "pathview_output"
dir.create(pathview_dir, showWarnings = FALSE)
old_wd <- getwd()
setwd(pathview_dir)



pv.out <- lapply(
  top_pathways,
  function(pid) {
    pathview(
      gene.data  = foldChanges,
      pathway.id = pid,
      species    = "mmu"
    )
  }
)

## 'select()' returned 1:1 mapping between keys and columns

## Info: Working in directory H:/maryam/R-Youtube/DESeq2_pipline/pathview_output

## Info: Writing image file mmu04710.pathview.png

## 'select()' returned 1:1 mapping between keys and columns

## Info: Working in directory H:/maryam/R-Youtube/DESeq2_pipline/pathview_output

## Info: Writing image file mmu00900.pathview.png

## 'select()' returned 1:1 mapping between keys and columns

## Info: Working in directory H:/maryam/R-Youtube/DESeq2_pipline/pathview_output

## Info: Writing image file mmu00100.pathview.png

setwd(old_wd)


# For reporting
filenames <- list.files(path =pathview_dir, pattern = ".*pathview.png", full.names = TRUE )

knitr:: include_graphics(filenames)

DESeq2 Pipeline

Maryam

2026-09-18

1. Packages

2. Import count matrix and sample metadata

3. Prepare metadata (experimental design)

4. Quality check and sample matching

5. Prepare count matrix

6. Create DESeq2 dataset

7. Filter low expression genes

8. Normalization and statistical modeling

9. Extract differential expression results

10. Visualization

11. Functional annotation

12. Filter significant DEGs

13. Pathway analysis