library(Seurat)
## Registered S3 method overwritten by 'spatstat.geom':
## method from
## print.boxx cli
## Attaching SeuratObject
library(reactome.db)
## Loading required package: AnnotationDbi
## Loading required package: stats4
## Loading required package: BiocGenerics
## Loading required package: parallel
##
## Attaching package: 'BiocGenerics'
## The following objects are masked from 'package:parallel':
##
## clusterApply, clusterApplyLB, clusterCall, clusterEvalQ,
## clusterExport, clusterMap, parApply, parCapply, parLapply,
## parLapplyLB, parRapply, parSapply, parSapplyLB
## The following objects are masked from 'package:stats':
##
## IQR, mad, sd, var, xtabs
## The following objects are masked from 'package:base':
##
## anyDuplicated, append, as.data.frame, basename, cbind, colnames,
## dirname, do.call, duplicated, eval, evalq, Filter, Find, get, grep,
## grepl, intersect, is.unsorted, lapply, Map, mapply, match, mget,
## order, paste, pmax, pmax.int, pmin, pmin.int, Position, rank,
## rbind, Reduce, rownames, sapply, setdiff, sort, table, tapply,
## union, 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 object is masked from 'package:base':
##
## expand.grid
library(ReactomePA)
##
## ReactomePA v1.34.0 For help: https://guangchuangyu.github.io/ReactomePA
##
## If you use ReactomePA in published research, please cite:
## Guangchuang Yu, Qing-Yu He. ReactomePA: an R/Bioconductor package for reactome pathway analysis and visualization. Molecular BioSystems 2016, 12(2):477-479
library(ReactomeGSA)
library(ReactomeGSA.data)
## Loading required package: limma
##
## Attaching package: 'limma'
## The following object is masked from 'package:BiocGenerics':
##
## plotMA
## Loading required package: edgeR
library(ggplot2)
pbmc <- readRDS("/mnt/nectar_volume/home/eraz0001/KELLY 2020/E11.5/Final_15_clusters.RDS")
library(ReactomeGSA)
ReactomeResults <- analyse_sc_clusters(pbmc, verbose = TRUE)
## Calculating average cluster expression...
## Converting expression data to string... (This may take a moment)
## Conversion complete
## Submitting request to Reactome API...
## Compressing request data...
## Reactome Analysis submitted succesfully
## Mapping identifiers...
## Performing gene set analysis using ssGSEA
## Analysing dataset 'Seurat' using ssGSEA
## Retrieving result...
ReactomeResults
## ReactomeAnalysisResult object
## Reactome Release: 78
## Results:
## - Seurat:
## 1741 pathways
## 9697 fold changes for genes
## No Reactome visualizations available
## ReactomeAnalysisResult
pathway_expression <- pathways(ReactomeResults)
colnames(pathway_expression) <- gsub("\\.Seurat", "", colnames(pathway_expression))
pathway_expression[1:20,] #you can choose your number of results here
## Name
## R-HSA-1059683 Interleukin-6 signaling
## R-HSA-109606 Intrinsic Pathway for Apoptosis
## R-HSA-109703 PKB-mediated events
## R-HSA-109704 PI3K Cascade
## R-HSA-110056 MAPK3 (ERK1) activation
## R-HSA-110312 Translesion synthesis by REV1
## R-HSA-110313 Translesion synthesis by Y family DNA polymerases bypasses lesions on DNA template
## R-HSA-110314 Recognition of DNA damage by PCNA-containing replication complex
## R-HSA-110320 Translesion Synthesis by POLH
## R-HSA-110328 Recognition and association of DNA glycosylase with site containing an affected pyrimidine
## R-HSA-110329 Cleavage of the damaged pyrimidine
## R-HSA-110330 Recognition and association of DNA glycosylase with site containing an affected purine
## R-HSA-110331 Cleavage of the damaged purine
## R-HSA-110357 Displacement of DNA glycosylase by APEX1
## R-HSA-110362 POLB-Dependent Long Patch Base Excision Repair
## R-HSA-110373 Resolution of AP sites via the multiple-nucleotide patch replacement pathway
## R-HSA-110381 Resolution of AP sites via the single-nucleotide replacement pathway
## R-HSA-111367 SLBP independent Processing of Histone Pre-mRNAs
## R-HSA-111446 Activation of BIM and translocation to mitochondria
## R-HSA-111447 Activation of BAD and translocation to mitochondria
## Chondrocytes DCh Endothelial_Cells Fibroblasts Hoxd13
## R-HSA-1059683 0.05479127 0.05007248 0.09365249 0.05369036 0.05488407
## R-HSA-109606 0.11346100 0.12536901 0.12821101 0.12444170 0.12345441
## R-HSA-109703 0.02206575 0.06213229 0.06229419 -0.01645497 0.03492748
## R-HSA-109704 0.04686066 0.03525082 0.04560987 0.02679795 0.02512239
## R-HSA-110056 0.07126449 0.10086349 0.13256362 0.10765450 0.09933605
## R-HSA-110312 0.12857521 0.13743074 0.13943199 0.15148118 0.13568893
## R-HSA-110313 0.12595859 0.13385880 0.12876995 0.13943208 0.13195110
## R-HSA-110314 0.15704390 0.16299796 0.15958918 0.17105583 0.16700481
## R-HSA-110320 0.14207824 0.14654451 0.14587102 0.15155373 0.14767426
## R-HSA-110328 0.16449886 0.16945003 0.15934902 0.16617762 0.17084387
## R-HSA-110329 0.16449886 0.16945003 0.15934902 0.16617762 0.17084387
## R-HSA-110330 0.16332859 0.17141061 0.15698979 0.16455967 0.17475776
## R-HSA-110331 0.16332859 0.17141061 0.15698979 0.16455967 0.17475776
## R-HSA-110357 0.12318028 0.10431029 0.12655799 0.14171906 0.10268606
## R-HSA-110362 0.16837690 0.18174567 0.15952228 0.17113251 0.19313621
## R-HSA-110373 0.14342771 0.15273373 0.14952602 0.15869078 0.16435399
## R-HSA-110381 0.22524221 0.22168281 0.21008974 0.23938980 0.22484495
## R-HSA-111367 0.22067809 0.25783292 0.23422126 0.23372038 0.26951142
## R-HSA-111446 0.18437808 0.18297707 0.17330806 0.18594487 0.17043102
## R-HSA-111447 0.11595909 0.12514280 0.12665578 0.13417839 0.11601806
## HtCh Immune_cells Osteoblast Osteoblasts Osteocytes
## R-HSA-1059683 0.05432962 0.08945233 0.05256261 0.05213726 0.052964394
## R-HSA-109606 0.12341748 0.11696063 0.12506132 0.12107789 0.119324259
## R-HSA-109703 0.03494138 0.02067297 0.04109390 0.07063454 -0.001872312
## R-HSA-109704 0.03340717 0.02425637 0.03382052 0.03641556 0.024012218
## R-HSA-110056 0.10182775 0.10892491 0.09802541 0.08987214 0.099407203
## R-HSA-110312 0.13350916 0.13755787 0.13288484 0.13657826 0.143595313
## R-HSA-110313 0.13183685 0.13072737 0.13212776 0.13366992 0.131477026
## R-HSA-110314 0.16011733 0.15969585 0.16172550 0.16470106 0.168655916
## R-HSA-110320 0.14413394 0.14587930 0.14666407 0.14989513 0.146178439
## R-HSA-110328 0.16697481 0.16711425 0.16797280 0.16711095 0.172084333
## R-HSA-110329 0.16697481 0.16711425 0.16797280 0.16711095 0.172084333
## R-HSA-110330 0.17183194 0.17211633 0.16765285 0.16623558 0.171859981
## R-HSA-110331 0.17183194 0.17211633 0.16765285 0.16623558 0.171859981
## R-HSA-110357 0.09499640 0.07578329 0.11876292 0.11676171 0.133232180
## R-HSA-110362 0.18435947 0.17391125 0.18242884 0.18668692 0.174513839
## R-HSA-110373 0.15267014 0.16388032 0.15400759 0.15667835 0.159517924
## R-HSA-110381 0.21607315 0.21243002 0.21610450 0.21783414 0.220451027
## R-HSA-111367 0.26379227 0.23865114 0.25090827 0.25039390 0.242837454
## R-HSA-111446 0.18496181 0.14866817 0.17976850 0.17825642 0.183436413
## R-HSA-111447 0.11639849 0.10201890 0.12033999 0.11752817 0.124655756
## Terminal_HTCh
## R-HSA-1059683 0.063144195
## R-HSA-109606 0.118196601
## R-HSA-109703 -0.005252028
## R-HSA-109704 0.046423952
## R-HSA-110056 0.078530117
## R-HSA-110312 0.121279136
## R-HSA-110313 0.120953608
## R-HSA-110314 0.145993453
## R-HSA-110320 0.134099989
## R-HSA-110328 0.159176366
## R-HSA-110329 0.159176366
## R-HSA-110330 0.154303680
## R-HSA-110331 0.154303680
## R-HSA-110357 0.134751174
## R-HSA-110362 0.159379047
## R-HSA-110373 0.127910699
## R-HSA-110381 0.216991400
## R-HSA-111367 0.233113915
## R-HSA-111446 0.195852525
## R-HSA-111447 0.116505634
max_difference <- do.call(rbind, apply(pathway_expression, 1, function(row) {
values <- as.numeric(row[2:length(row)])
return(data.frame(name = row[1], min = min(values), max = max(values)))
}))
max_difference$diff <- max_difference$max - max_difference$min
max_difference <- max_difference[order(max_difference$diff, decreasing = T), ]
head(max_difference)
## name min
## R-HSA-140180 COX reactions -0.44827590
## R-HSA-211994 Sterols are 12-hydroxylated by CYP8B1 -0.31045130
## R-HSA-190374 FGFR1c and Klotho ligand binding and activation -0.07386122
## R-HSA-9636003 NEIL3-mediated resolution of ICLs -0.48382290
## R-HSA-2142696 Synthesis of Hepoxilins (HX) and Trioxilins (TrX) -0.04150702
## R-HSA-8964540 Alanine metabolism -0.21424010
## max diff
## R-HSA-140180 0.31864623 0.7669221
## R-HSA-211994 0.33375905 0.6442103
## R-HSA-190374 0.48041720 0.5542784
## R-HSA-9636003 0.06545338 0.5492763
## R-HSA-2142696 0.43433376 0.4758408
## R-HSA-8964540 0.25446999 0.4687101
plot_gsva_pathway(ReactomeResults, pathway_id = rownames(max_difference)[1])
plot_gsva_heatmap(ReactomeResults, max_pathways = 15, margins = c(6,20))
gsva_result = ReactomeResults
plot_gsva_pca(ReactomeResults)
for (i in seq(1, 12, by = 1)){
aa <- plot_gsva_pathway(ReactomeResults, pathway_id = rownames(max_difference)[i])
print(aa)
}
relevant_pathways <- c("R-HSA-140180", "R-HSA-9673163", "R-HSA-190371", "R-HSA-9636003")
plot_gsva_heatmap(ReactomeResults,
pathway_ids = relevant_pathways, # limit to these pathways
margins = c(6,30), # adapt the figure margins in heatmap.2
dendrogram = "col", # only plot column dendrogram
scale = "row", # scale for each pathway
key = FALSE, # don't display the color key
lwid=c(0.1,4))
