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:
rawData <- read.csv("GLDS-102_rna_seq_Normalized_Counts.csv",row.names = 1)
group <- factor(c("1","1","1","1","1","1","2","2","2","2","2","2"))
3.making DGEList:
dgeGlm <- DGEList(
counts = rawData,
group = group
)
str(dgeGlm)
## 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 "1","2": 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:
keep <- rowSums(
cpm(dgeGlm)>2
) >= 4
dgeGlm <- dgeGlm[keep,]
5.making design matrix
design <- model.matrix(~group)
design
## (Intercept) group2
## 1 1 0
## 2 1 0
## 3 1 0
## 4 1 0
## 5 1 0
## 6 1 0
## 7 1 1
## 8 1 1
## 9 1 1
## 10 1 1
## 11 1 1
## 12 1 1
## attr(,"assign")
## [1] 0 1
## attr(,"contrasts")
## attr(,"contrasts")$group
## [1] "contr.treatment"
6.Estimate Dispersion:
dgeGlmComDisp <- estimateGLMCommonDisp(dgeGlm,design,verbose = TRUE)
## Disp = 0.02648 , BCV = 0.1627
dgeGlmTrendedDisp <- estimateGLMTrendedDisp(dgeGlmComDisp,design)
dgeGlmTagDisp <- estimateGLMTagwiseDisp(dgeGlmTrendedDisp,design)
plotBCV(dgeGlmTagDisp)
7. 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)
head(ttGlm)
## Coefficient: group2
## logFC logCPM LR PValue FDR
## ENSMUSG00000026077 -1.3621448 3.576328 79.72497 4.303282e-19 5.858058e-15
## ENSMUSG00000007872 0.8863833 5.525394 76.76969 1.921032e-18 1.307551e-14
## ENSMUSG00000021775 0.9537055 6.241961 63.05270 2.012493e-15 9.132021e-12
## ENSMUSG00000002250 -0.8282806 5.266124 62.46754 2.708731e-15 9.218490e-12
## ENSMUSG00000059824 2.2593713 4.561763 57.98741 2.638001e-14 7.182221e-11
## ENSMUSG00000074715 -1.9895168 3.805293 47.03450 6.974820e-12 1.582470e-08
8. Selection of significant genes
DGEs <- dgeGlm[
dgeGlm$FDR < 0.1,
]
summary(DGEs)
## Length Class Mode
## counts 0 -none- numeric
## samples 3 data.frame list
tagsGlm <- rownames(DGEs)
plotSmear(lrt,de.tags = tagsGlm)
## Warning in plot.xy(xy.coords(x, y), type = type, ...): "panel.first" is not a
## graphical parameter
hist2 <- dgeGlm[dgeGlm$FDR < 0.1,]
#Saving
write.csv(hist2,"Mouse_EdgeR_Results_Zero_vs_1_Exercise.csv")
9. 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 columns
ttGlm2$entrez = mapIds(org.Mm.eg.db,
keys = row.names(ttGlm2),
column = "ENTREZID",
keytype = "ENSEMBL",
multiVals = "first")
## 'select()' returned 1:many mapping between keys and columns
ttGlm2$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 LR PValue FDR
## ENSMUSG00000026077 -1.3621448 3.576328 79.72497 4.303282e-19 5.858058e-15
## ENSMUSG00000007872 0.8863833 5.525394 76.76969 1.921032e-18 1.307551e-14
## ENSMUSG00000021775 0.9537055 6.241961 63.05270 2.012493e-15 9.132021e-12
## ENSMUSG00000002250 -0.8282806 5.266124 62.46754 2.708731e-15 9.218490e-12
## ENSMUSG00000059824 2.2593713 4.561763 57.98741 2.638001e-14 7.182221e-11
## ENSMUSG00000074715 -1.9895168 3.805293 47.03450 6.974820e-12 1.582470e-08
## symbol entrez
## ENSMUSG00000026077 Npas2 18143
## ENSMUSG00000007872 Id3 15903
## ENSMUSG00000021775 Nr1d2 353187
## ENSMUSG00000002250 Ppard 19015
## ENSMUSG00000059824 Dbp 13170
## ENSMUSG00000074715 Ccl28 56838
## name
## ENSMUSG00000026077 neuronal PAS domain protein 2
## ENSMUSG00000007872 inhibitor of DNA binding 3
## ENSMUSG00000021775 nuclear receptor subfamily 1, group D, member 2
## ENSMUSG00000002250 peroxisome proliferator activator receptor delta
## ENSMUSG00000059824 D site albumin promoter binding protein
## ENSMUSG00000074715 C-C motif chemokine ligand 28
10. 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(go.subs.mm)
data(sigmet.idx.mm)
#preparing data for gage
foldchanges <- ttGlm2$logFC
names(foldchanges) <- ttGlm2$entrez
data(kegg.sets.mm)
kegg.mm <- kegg.sets.mm
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 7.695685e-09 6.112711
## mmu04350 TGF-beta signaling pathway 3.671996e-03 2.724347
## mmu00982 Drug metabolism - cytochrome P450 9.942013e-03 2.373913
## mmu04330 Notch signaling pathway 2.119323e-02 2.061627
## mmu00980 Metabolism of xenobiotics by cytochrome P450 2.570064e-02 1.982056
## mmu00830 Retinol metabolism 2.672841e-02 1.972303
## p.val q.val
## mmu03010 Ribosome 7.695685e-09 1.223614e-06
## mmu04350 TGF-beta signaling pathway 3.671996e-03 2.919237e-01
## mmu00982 Drug metabolism - cytochrome P450 9.942013e-03 5.269267e-01
## mmu04330 Notch signaling pathway 2.119323e-02 7.083030e-01
## mmu00980 Metabolism of xenobiotics by cytochrome P450 2.570064e-02 7.083030e-01
## mmu00830 Retinol metabolism 2.672841e-02 7.083030e-01
## set.size exp1
## mmu03010 Ribosome 88 7.695685e-09
## mmu04350 TGF-beta signaling pathway 65 3.671996e-03
## mmu00982 Drug metabolism - cytochrome P450 43 9.942013e-03
## mmu04330 Notch signaling pathway 44 2.119323e-02
## mmu00980 Metabolism of xenobiotics by cytochrome P450 36 2.570064e-02
## mmu00830 Retinol metabolism 30 2.672841e-02
##
## $less
## p.geomean stat.mean
## mmu04145 Phagosome 0.002100896 -2.890969
## mmu00910 Nitrogen metabolism 0.003869883 -2.843207
## mmu04640 Hematopoietic cell lineage 0.004541940 -2.677574
## mmu04110 Cell cycle 0.005802394 -2.545883
## mmu00100 Steroid biosynthesis 0.006459181 -2.658249
## mmu04622 RIG-I-like receptor signaling pathway 0.015171792 -2.201562
## p.val q.val set.size
## mmu04145 Phagosome 0.002100896 0.2054020 119
## mmu00910 Nitrogen metabolism 0.003869883 0.2054020 17
## mmu04640 Hematopoietic cell lineage 0.004541940 0.2054020 41
## mmu04110 Cell cycle 0.005802394 0.2054020 108
## mmu00100 Steroid biosynthesis 0.006459181 0.2054020 15
## mmu04622 RIG-I-like receptor signaling pathway 0.015171792 0.3543823 48
## exp1
## mmu04145 Phagosome 0.002100896
## mmu00910 Nitrogen metabolism 0.003869883
## mmu04640 Hematopoietic cell lineage 0.004541940
## mmu04110 Cell cycle 0.005802394
## mmu00100 Steroid biosynthesis 0.006459181
## mmu04622 RIG-I-like receptor signaling pathway 0.015171792
##
## $stats
## stat.mean exp1
## mmu03010 Ribosome 6.112711 6.112711
## mmu04350 TGF-beta signaling pathway 2.724347 2.724347
## mmu00982 Drug metabolism - cytochrome P450 2.373913 2.373913
## mmu04330 Notch signaling pathway 2.061627 2.061627
## mmu00980 Metabolism of xenobiotics by cytochrome P450 1.982056 1.982056
## mmu00830 Retinol metabolism 1.972303 1.972303
greaters <- keggres$greater
lessers <- keggres$less
11.Drawing KEGG pathways
keggrespathways = data.frame(id = rownames(greaters), greaters) %>%
tibble::as.tibble() %>%
filter(row_number() <=5 )%>%
.$id %>%
as.character()
## Warning: `as.tibble()` was deprecated in tibble 2.0.0.
## ℹ Please use `as_tibble()` instead.
## ℹ The signature and semantics have changed, see `?as_tibble`.
## This warning is displayed once per session.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.
## [1] "mmu03010 Ribosome"
## [2] "mmu04350 TGF-beta signaling pathway"
## [3] "mmu00982 Drug metabolism - cytochrome P450"
## [4] "mmu04330 Notch signaling pathway"
## [5] "mmu00980 Metabolism of xenobiotics by cytochrome P450"
keggresids = substr(keggrespathways, start = 1, stop = 8)
keggresids
## [1] "mmu03010" "mmu04350" "mmu00982" "mmu04330" "mmu00980"
plot_pathway = function(pid) pathview(gene.data = foldchanges, pathway.id = pid, species = "mmu", new.signature = FALSE)
#plot multiple pathway
keggresids3 <- sub("mmu", "", keggresids)
tmp = sapply(keggresids3, 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
## Info: Writing image file mmu03010.pathview.png
## 'select()' returned 1:1 mapping between keys and columns
## Info: Working in directory H:/maryam/R-Youtube
## Info: Writing image file mmu04350.pathview.png
## 'select()' returned 1:1 mapping between keys and columns
## Info: Working in directory H:/maryam/R-Youtube
## Info: Writing image file mmu00982.pathview.png
## 'select()' returned 1:1 mapping between keys and columns
## Info: Working in directory H:/maryam/R-Youtube
## Info: Writing image file mmu04330.pathview.png
## 'select()' returned 1:1 mapping between keys and columns
## Info: Working in directory H:/maryam/R-Youtube
## Info: Writing image file mmu00980.pathview.png
keggresspathways2 = data.frame(id=rownames(lessers),lessers) %>%
tibble::as.tibble()%>%
filter(row_number() <=5) %>%
.$id %>%
as.character()
keggresspathways2
## [1] "mmu04145 Phagosome" "mmu00910 Nitrogen metabolism"
## [3] "mmu04640 Hematopoietic cell lineage" "mmu04110 Cell cycle"
## [5] "mmu00100 Steroid biosynthesis"
keggresids2 = substr(keggresspathways2, start = 1 , stop = 8)
keggresids2
## [1] "mmu04145" "mmu00910" "mmu04640" "mmu04110" "mmu00100"
plot_pathway2 <- function(pid) pathview(gene.data = foldchanges, pathway.id = pid, species = "mmu", new.signature = FALSE)
keggresids4 <- sub("mmu","",keggresids2)
tmp2 = sapply(keggresids4, 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
## Info: Writing image file mmu04145.pathview.png
## 'select()' returned 1:1 mapping between keys and columns
## Info: Working in directory H:/maryam/R-Youtube
## Info: Writing image file mmu00910.pathview.png
## 'select()' returned 1:1 mapping between keys and columns
## Info: Working in directory H:/maryam/R-Youtube
## Info: Writing image file mmu04640.pathview.png
## 'select()' returned 1:1 mapping between keys and columns
## Info: Working in directory H:/maryam/R-Youtube
## Info: Writing image file mmu04110.pathview.png
## 'select()' returned 1:1 mapping between keys and columns
## Info: Working in directory H:/maryam/R-Youtube
## Info: Writing image file mmu00100.pathview.png
## 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.image
filenames <- list.files(path ="H:/maryam/R-Youtube", pattern = ".*pathview.png" )
all_images <- lapply(filenames, load.image)
knitr:: include_graphics(filenames)
