EdgeR execise
Maryam
2026-07-21
1.Packages:
## Loading required package: limma##
## Attaching package: 'dplyr'## The following objects are masked from 'package:stats':
##
## filter, lag## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, union## 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 object is masked from 'package:dplyr':
##
## explain## 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 object is masked from 'package:dplyr':
##
## combine## The following object is masked from 'package:limma':
##
## plotMA## 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## 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")'.## Loading required package: IRanges## Loading required package: S4Vectors##
## Attaching package: 'S4Vectors'## The following objects are masked from 'package:dplyr':
##
## first, rename## The following object is masked from 'package:utils':
##
## findMatches## The following objects are masked from 'package:base':
##
## expand.grid, I, unname##
## Attaching package: 'IRanges'## The following objects are masked from 'package:dplyr':
##
## collapse, desc, slice## The following object is masked from 'package:grDevices':
##
## windows##
## Attaching package: 'AnnotationDbi'## The following object is masked from 'package:dplyr':
##
## select## 2.Raw RNA-seq count data and grouping the data:
We should use Unnormalized counts, but I couldn’t find it so lets go on with the normalized one this time. and it’s better to grouping it with meaningful names.
3.making DGEList:
## Formal class 'DGEList' [package "edgeR"] with 1 slot
## ..@ .Data:List of 2
## .. ..$ : num [1:24035, 1:12] 2976.8 59.8 21.2 40.1 0 ...
## .. .. ..- attr(*, "dimnames")=List of 2
## .. .. .. ..$ : chr [1:24035] "ENSMUSG00000000001" "ENSMUSG00000000028" "ENSMUSG00000000031" "ENSMUSG00000000037" ...
## .. .. .. ..$ : chr [1:12] "Mmus_C57.6J_KDN_FLT_Rep1_M23" "Mmus_C57.6J_KDN_FLT_Rep2_M24" "Mmus_C57.6J_KDN_FLT_Rep3_M25" "Mmus_C57.6J_KDN_FLT_Rep4_M26" ...
## .. ..$ :'data.frame': 12 obs. of 3 variables:
## .. .. ..$ group : Factor w/ 2 levels "control","treatment": 1 1 1 1 1 1 2 2 2 2 ...
## .. .. ..$ lib.size : num [1:12] 40266365 40740336 37739541 39196969 36820645 ...
## .. .. ..$ norm.factors: num [1:12] 1 1 1 1 1 1 1 1 1 1 ...
## ..$ names: chr [1:2] "counts" "samples"4.Removing low expressed Genes:
filterByExpr() is the current recommended approach (adapts the threshold to library sizes and group sizes automatically), rather than a manual CPM cutoff.manual rule (cpm > 2 in >= 4 samples) is not wrong, but filterByExpr is now the standard and more robust.
# It is old way
#keep <- rowSums(
# cpm(dgeGlm)>2
#) >= 4
#dgeGlm <- dgeGlm[keep,]
# Its is new way of filtering low expressed genes
keep <- filterByExpr(dgeGlm, group = group)
dgeGlm <- dgeGlm[keep, , keep.lib.sizes = FALSE]
table(keep)## keep
## FALSE TRUE
## 6268 177675.Normalization
calcNormFactors() computes TMM normalization factors to correct for RNA composition differences between libraries. Without this step, edgeR uses raw library sizes only, which can bias logFC and p-values, especially when a few highly expressed genes differ strongly between groups.
## calcNormFactors has been renamed to normLibSizes## group lib.size norm.factors
## Mmus_C57.6J_KDN_FLT_Rep1_M23 control 40244279 0.9372231
## Mmus_C57.6J_KDN_FLT_Rep2_M24 control 40715743 0.9287846
## Mmus_C57.6J_KDN_FLT_Rep3_M25 control 37713767 0.9990904
## Mmus_C57.6J_KDN_FLT_Rep4_M26 control 39171133 0.9571221
## Mmus_C57.6J_KDN_FLT_Rep5_M27 control 36794493 1.0176533
## Mmus_C57.6J_KDN_FLT_Rep6_M28 control 36447848 1.0293697
## Mmus_C57.6J_KDN_GC_Rep1_M33 treatment 36561141 1.0279594
## Mmus_C57.6J_KDN_GC_Rep2_M34 treatment 36525208 1.0353444
## Mmus_C57.6J_KDN_GC_Rep3_M35 treatment 39550907 0.9527398
## Mmus_C57.6J_KDN_GC_Rep4_M36 treatment 36322742 1.0361061
## Mmus_C57.6J_KDN_GC_Rep5_M37 treatment 37640633 1.0006113
## Mmus_C57.6J_KDN_GC_Rep6_M38 treatment 34634568 1.09092256.Exploratory data analysis
MDS plot = quick way to see whether samples cluster by group as expected, and to catch potential outliers/mislabeled samples BEFORE testing for DE.
plotMDS(dgeGlm, col = as.numeric(group), main = "MDS Plot")
legend("topright", legend = levels(group), col = 1:length(levels(group)), pch = 16)
# 7.making design matrix
## (Intercept) grouptreatment
## Mmus_C57.6J_KDN_FLT_Rep1_M23 1 0
## Mmus_C57.6J_KDN_FLT_Rep2_M24 1 0
## Mmus_C57.6J_KDN_FLT_Rep3_M25 1 0
## Mmus_C57.6J_KDN_FLT_Rep4_M26 1 0
## Mmus_C57.6J_KDN_FLT_Rep5_M27 1 0
## Mmus_C57.6J_KDN_FLT_Rep6_M28 1 0
## Mmus_C57.6J_KDN_GC_Rep1_M33 1 1
## Mmus_C57.6J_KDN_GC_Rep2_M34 1 1
## Mmus_C57.6J_KDN_GC_Rep3_M35 1 1
## Mmus_C57.6J_KDN_GC_Rep4_M36 1 1
## Mmus_C57.6J_KDN_GC_Rep5_M37 1 1
## Mmus_C57.6J_KDN_GC_Rep6_M38 1 1
## attr(,"assign")
## [1] 0 1
## attr(,"contrasts")
## attr(,"contrasts")$group
## [1] "contr.treatment"8.Estimate Dispersion:
EstimateDisp() replaces the three legacy calls (estimateGLMCommonDisp / estimateGLMTrendedDisp / estimateGLMTagwiseDisp) with a single, more numerically stable function that does all three (common -> trended -> tagwise) internally. Functionally equivalent to your original three lines, just the modern/recommended API.
#dgeGlmComDisp <- estimateGLMCommonDisp(dgeGlm,design,verbose = TRUE)
#dgeGlmTrendedDisp <- estimateGLMTrendedDisp(dgeGlmComDisp,design)
#dgeGlmTagDisp <- estimateGLMTagwiseDisp(dgeGlmTrendedDisp,design)
dgeGlmDisp <- estimateDisp(dgeGlm,design)
plotBCV(dgeGlmDisp)
9. Differential Expression
#Generalized Linear Model for each genes
#fit <- glmFit(dgeGlmTagDisp,design)
#Likelihood Ratio Test
#lrt <- glmLRT(fit,coef = 2)
#results
#ttGlm <- topTags(lrt,n = Inf)
fit <- glmQLFit(dgeGlmDisp,design)
qlf <- glmQLFTest(fit,coef = 2)
ttGlm <- topTags(qlf, n = Inf)
head(ttGlm)## Coefficient: grouptreatment
## logFC logCPM F PValue FDR
## ENSMUSG00000026077 -1.4395242 3.595296 95.30252 1.073182e-07 0.001906722
## ENSMUSG00000002250 -0.9003315 5.279324 76.06173 4.310147e-07 0.003828919
## ENSMUSG00000007872 0.8221989 5.520618 61.21486 1.575672e-06 0.007799390
## ENSMUSG00000021775 0.8892720 6.231754 60.09245 1.755927e-06 0.007799390
## ENSMUSG00000040998 -0.4974178 8.601906 52.73331 3.733239e-06 0.012990495
## ENSMUSG00000027111 -0.4593931 7.744927 51.25792 4.386952e-06 0.01299049510. Selection of significant genes
## [1] 88 5tagsGlm <- rownames(DGEs)
#plotSmear(lrt,de.tags = tagsGlm)
#hist2 <- deGlm[deGlm$FDR < 0.1,]
{plotMD(qlf,status = decideTests(qlf, p.value = 0.05), main = "MD Plot")
abline(h = c(-1,1), col = "blue")}
11. Gene IDs conversion with mapIDs
ttGlm2 <- as.data.frame(ttGlm$table)
ttGlm2$symbol = mapIds(org.Mm.eg.db,
keys = row.names(ttGlm2),
column = "SYMBOL",
keytype = "ENSEMBL",
multiVals = "first"
)## 'select()' returned 1:many mapping between keys and columnsttGlm2$entrez = mapIds(org.Mm.eg.db,
keys = row.names(ttGlm2),
column = "ENTREZID",
keytype = "ENSEMBL",
multiVals = "first")## 'select()' returned 1:many mapping between keys and columnsttGlm2$name = mapIds(org.Mm.eg.db,
keys = row.names(ttGlm2),
column = "GENENAME",
keytype = "ENSEMBL",
multiVals = "first")## 'select()' returned 1:many mapping between keys and columns## logFC logCPM F PValue FDR symbol
## ENSMUSG00000026077 -1.4395242 3.595296 95.30252 1.073182e-07 0.001906722 Npas2
## ENSMUSG00000002250 -0.9003315 5.279324 76.06173 4.310147e-07 0.003828919 Ppard
## ENSMUSG00000007872 0.8221989 5.520618 61.21486 1.575672e-06 0.007799390 Id3
## ENSMUSG00000021775 0.8892720 6.231754 60.09245 1.755927e-06 0.007799390 Nr1d2
## ENSMUSG00000040998 -0.4974178 8.601906 52.73331 3.733239e-06 0.012990495 Npnt
## ENSMUSG00000027111 -0.4593931 7.744927 51.25792 4.386952e-06 0.012990495 Itga6
## entrez name
## ENSMUSG00000026077 18143 neuronal PAS domain protein 2
## ENSMUSG00000002250 19015 peroxisome proliferator activator receptor delta
## ENSMUSG00000007872 15903 inhibitor of DNA binding 3
## ENSMUSG00000021775 353187 nuclear receptor subfamily 1, group D, member 2
## ENSMUSG00000040998 114249 nephronectin
## ENSMUSG00000027111 16403 integrin alpha 612. Pathway Enrichment (KEGG)
## ## ##############################################################################
## 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(gageData)
data(kegg.sets.mm)
data(sigmet.idx.mm)
#preparing data for gage
ttGlm2_kegg <- ttGlm2[!is.na(ttGlm2$entrez),]
foldchanges <- ttGlm2_kegg$logFC
names(foldchanges) <- ttGlm2_kegg$entrez
kegg.mm.sigmet <- kegg.sets.mm[sigmet.idx.mm]
keggres = gage(
foldchanges,
gsets = kegg.mm.sigmet,
same.dir = TRUE
)
lapply(keggres, head)## $greater
## p.geomean stat.mean
## mmu03010 Ribosome 1.435044e-08 5.897486
## mmu04350 TGF-beta signaling pathway 7.063114e-03 2.483894
## mmu00982 Drug metabolism - cytochrome P450 1.351800e-02 2.242716
## mmu00980 Metabolism of xenobiotics by cytochrome P450 1.694353e-02 2.154660
## mmu00830 Retinol metabolism 2.586426e-02 1.977727
## mmu04360 Axon guidance 3.297250e-02 1.847423
## p.val q.val
## mmu03010 Ribosome 1.435044e-08 2.281720e-06
## mmu04350 TGF-beta signaling pathway 7.063114e-03 5.615176e-01
## mmu00982 Drug metabolism - cytochrome P450 1.351800e-02 6.735054e-01
## mmu00980 Metabolism of xenobiotics by cytochrome P450 1.694353e-02 6.735054e-01
## mmu00830 Retinol metabolism 2.586426e-02 7.529994e-01
## mmu04360 Axon guidance 3.297250e-02 7.529994e-01
## set.size exp1
## mmu03010 Ribosome 95 1.435044e-08
## mmu04350 TGF-beta signaling pathway 74 7.063114e-03
## mmu00982 Drug metabolism - cytochrome P450 54 1.351800e-02
## mmu00980 Metabolism of xenobiotics by cytochrome P450 47 1.694353e-02
## mmu00830 Retinol metabolism 40 2.586426e-02
## mmu04360 Axon guidance 122 3.297250e-02
##
## $less
## p.geomean stat.mean p.val q.val
## mmu04110 Cell cycle 0.003920960 -2.682110 0.003920960 0.4666428
## mmu04145 Phagosome 0.008798041 -2.389143 0.008798041 0.4666428
## mmu00910 Nitrogen metabolism 0.008804581 -2.497676 0.008804581 0.4666428
## mmu00100 Steroid biosynthesis 0.019812285 -2.160624 0.019812285 0.6820620
## mmu04115 p53 signaling pathway 0.021448489 -2.046456 0.021448489 0.6820620
## mmu03440 Homologous recombination 0.031887125 -1.894921 0.031887125 0.8063396
## set.size exp1
## mmu04110 Cell cycle 117 0.003920960
## mmu04145 Phagosome 135 0.008798041
## mmu00910 Nitrogen metabolism 18 0.008804581
## mmu00100 Steroid biosynthesis 16 0.019812285
## mmu04115 p53 signaling pathway 62 0.021448489
## mmu03440 Homologous recombination 27 0.031887125
##
## $stats
## stat.mean exp1
## mmu03010 Ribosome 5.897486 5.897486
## mmu04350 TGF-beta signaling pathway 2.483894 2.483894
## mmu00982 Drug metabolism - cytochrome P450 2.242716 2.242716
## mmu00980 Metabolism of xenobiotics by cytochrome P450 2.154660 2.154660
## mmu00830 Retinol metabolism 1.977727 1.977727
## mmu04360 Axon guidance 1.847423 1.84742313.Drawing KEGG pathways
# Top 5 up-regulated pathways
top_up = greaters %>%
rownames_to_column("id") %>%
filter(!is.na(p.geomean)) %>%
slice_head(n = 5) %>%
pull(id)
# Top 5 down-regulated pathways
top_down <- lessers %>%
rownames_to_column("id") %>%
filter(!is.na(p.geomean))%>%
slice_head(n = 5)%>%
pull(id)
# Extract the KEGG pathway IDs (first 8 characters, e.g. "mmu04110")
ids_up <- substr(top_up,1,8)
ids_up## [1] "mmu03010" "mmu04350" "mmu00982" "mmu00980" "mmu00830"## [1] "mmu04110" "mmu04145" "mmu00910" "mmu00100" "mmu04115"# pathview() wants the numeric part only ("04110" not "mmu04110")
ids_up_num <- sub("mmu","",ids_up)
ids_down_num <- sub("mmu","",ids_down)
out_dir <- file.path(getwd(), "pathway_plots")
dir.create(out_dir, showWarnings = FALSE)
old_wd <- setwd(out_dir) # pathview writes files to the working directory
invisible(lapply(ids_up_num, function(pid)
tryCatch(pathview(gene.data = foldChanges, pathway.id = pid, species = "mmu"),
error = function(e) message("Failed: ", pid, " - ", e$message))))## Failed: 03010 - object 'foldChanges' not found## Failed: 04350 - object 'foldChanges' not found## Failed: 00982 - object 'foldChanges' not found## Failed: 00980 - object 'foldChanges' not found## Failed: 00830 - object 'foldChanges' not foundinvisible(lapply(ids_down_num, function(pid)
tryCatch(pathview(gene.data = foldChanges, pathway.id = pid, species = "mmu"),
error = function(e) message("Failed: ", pid, " - ", e$message))))## Failed: 04110 - object 'foldChanges' not found## Failed: 04145 - object 'foldChanges' not found## Failed: 00910 - object 'foldChanges' not found## Failed: 00100 - object 'foldChanges' not found## Failed: 04115 - object 'foldChanges' not found## Warning: package 'imager' was built under R version 4.6.1## Loading required package: magrittr##
## Attaching package: 'imager'## The following object is masked from 'package:magrittr':
##
## add## The following objects are masked from 'package:IRanges':
##
## resize, width## The following object is masked from 'package:S4Vectors':
##
## width## The following object is masked from 'package:Biobase':
##
## channel## The following object is masked from 'package:BiocGenerics':
##
## width## The following object is masked from 'package:dplyr':
##
## where## The following objects are masked from 'package:stats':
##
## convolve, spectrum## The following object is masked from 'package:graphics':
##
## frame## The following object is masked from 'package:base':
##
## save.imagepathway_files <- list.files(path = out_dir, pattern = "*.pathview.png$",
full.names = TRUE)
if (length(pathway_files) > 0) {
all_images <- lapply(pathway_files, load.image)
knitr::include_graphics(pathway_files)
} else {
message("No pathview PNG files found in: ", out_dir)
}
