library("edgeR")
library("dplyr")
library("AnnotationDbi")
library("org.Mm.eg.db")
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"))
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"
str(group)
## Factor w/ 2 levels "1","2": 1 1 1 1 1 1 2 2 2 2 ...
keep <- rowSums(cpm(dgeGlm)>2) >=4
dgeGlm <- dgeGlm[keep, ]
names(dgeGlm)
## [1] "counts" "samples"
dgeGlm[["samples"]]
## group lib.size norm.factors
## Mmus_C57.6J_KDN_FLT_Rep1_M23 1 40266365 1
## Mmus_C57.6J_KDN_FLT_Rep2_M24 1 40740336 1
## Mmus_C57.6J_KDN_FLT_Rep3_M25 1 37739541 1
## Mmus_C57.6J_KDN_FLT_Rep4_M26 1 39196969 1
## Mmus_C57.6J_KDN_FLT_Rep5_M27 1 36820645 1
## Mmus_C57.6J_KDN_FLT_Rep6_M28 1 36475638 1
## Mmus_C57.6J_KDN_GC_Rep1_M33 2 36585661 1
## Mmus_C57.6J_KDN_GC_Rep2_M34 2 36550490 1
## Mmus_C57.6J_KDN_GC_Rep3_M35 2 39573830 1
## Mmus_C57.6J_KDN_GC_Rep4_M36 2 36346170 1
## Mmus_C57.6J_KDN_GC_Rep5_M37 2 37663039 1
## Mmus_C57.6J_KDN_GC_Rep6_M38 2 34662629 1
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"
dgeGlmComDisp <- estimateGLMCommonDisp(dgeGlm,design,verbose = TRUE)
## Disp = 0.02648 , BCV = 0.1627
dgeGlmTrendDisp <- estimateGLMTrendedDisp(dgeGlmComDisp, design)
dgeGlmTagDisp <- estimateGLMTagwiseDisp(dgeGlmTrendDisp,design)
plotBCV(dgeGlmTagDisp)
fit <- glmFit(dgeGlmTagDisp,design)
colnames(coef(fit))
## [1] "(Intercept)" "group2"
lrt <- glmLRT(fit, coef = 2) # LRT مخفف Likelihood Ratio Test
ttGlm <- topTags(lrt,n= Inf)
head(ttGlm$table)
## 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
class(ttGlm)
## [1] "TopTags"
## attr(,"package")
## [1] "edgeR"
#summary(deGlm <- decideTestsDGE(lrt,p = 0.05, adjust="fdr")) #it doesn't work in this version of edgeR
deGlm <- ttGlm$table
DGEs <- deGlm[
deGlm$FDR < 0.1 #&
#abs(deGlm$logFC) > 1
,
]
summary(DGEs)
## logFC logCPM LR PValue
## Min. :-2.1769 Min. : 0.6987 Min. : 9.916 Min. :0.000e+00
## 1st Qu.:-0.3112 1st Qu.: 4.2196 1st Qu.:11.113 1st Qu.:4.104e-05
## Median : 0.3044 Median : 5.2921 Median :13.011 Median :3.097e-04
## Mean : 0.1525 Mean : 5.3090 Mean :16.390 Mean :4.909e-04
## 3rd Qu.: 0.4561 3rd Qu.: 6.5851 3rd Qu.:16.828 3rd Qu.:8.572e-04
## Max. : 2.2594 Max. :11.2980 Max. :79.725 Max. :1.638e-03
## FDR
## Min. :0.000000
## 1st Qu.:0.009882
## Median :0.037644
## Mean :0.040870
## 3rd Qu.:0.069187
## Max. :0.099999
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
hits2 <- deGlm[deGlm$FDR < 0.1, ]
write.csv(hits2, "Mouse_EdgeR_Results_Zero_vs_1.csv")
columns(org.Mm.eg.db)
## [1] "ACCNUM" "ALIAS" "ENSEMBL" "ENSEMBLPROT" "ENSEMBLTRANS"
## [6] "ENTREZID" "ENZYME" "EVIDENCE" "EVIDENCEALL" "GENENAME"
## [11] "GENETYPE" "GO" "GOALL" "IPI" "MGI"
## [16] "ONTOLOGY" "ONTOLOGYALL" "PATH" "PFAM" "PMID"
## [21] "PROSITE" "REFSEQ" "SYMBOL" "UNIPROT"
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
head(ttGlm2)
## 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
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(gage)
##
library(gageData)
data(kegg.sets.mm)
data(go.subs.mm)
data(sigmet.idx.mm)
foldchanges <- ttGlm2$logFC
names(foldchanges) <- ttGlm2$entrez
# kegg.mm = kegg.gsets("mouse", id.type = "entrez") # the server doesn't work
data(kegg.sets.mm)
kegg.mm <- kegg.sets.mm
#kegg.mm.sigmet = kegg.mm$kg.sets[kegg.mm$sigmet.idx] #It doesn't work cause it's for old version
kegg.mm.sigmet <- kegg.sets.mm[sigmet.idx.mm]
#Lets get the results
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
keggrespathways = data.frame(id = rownames(keggres$greater), keggres$greater) %>%
tibble::as_tibble() %>%
filter(row_number()<=5) %>%
.$id %>%
as.character()
keggrespathways
## [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
keggresids2 <- sub("mmu", "", keggresids)
tmp = sapply(keggresids2[1], 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
keggrespathways = data.frame(id = rownames(keggres$less), keggres$less) %>%
tibble::as_tibble() %>%
filter(row_number()<=5) %>%
.$id %>%
as.character()
keggrespathways
## [1] "mmu04145 Phagosome" "mmu00910 Nitrogen metabolism"
## [3] "mmu04640 Hematopoietic cell lineage" "mmu04110 Cell cycle"
## [5] "mmu00100 Steroid biosynthesis"
keggresids = substr(keggrespathways, start = 1, stop = 8)
keggresids
## [1] "mmu04145" "mmu00910" "mmu04640" "mmu04110" "mmu00100"
plot_pathway = function(pid) pathview(gene.data = foldchanges, pathway.id = pid, species = "mmu", new.signature = FALSE)
#plot multiple pathway
keggresids2 <- sub("mmu", "", keggresids)
tmp = sapply(keggresids2[1], 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
library(imager)
## 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)
