## ----include = FALSE---------------------------------------------------------- knitr::opts_chunk$set( collapse = TRUE, comment = "#>" ) ## ----echo=FALSE--------------------------------------------------------------- library(BiocParallel) register(SerialParam()) ## ----eval=FALSE--------------------------------------------------------------- # score <- sum(colSums(E[p, ])**2) / length(p) ## ----message=FALSE------------------------------------------------------------ library(GEOquery) library(limma) gse200250 <- getGEO("GSE200250", AnnotGPL = TRUE)[[1]] es <- gse200250 es <- es[, grep("Th2_", es$title)] es$time <- as.numeric(gsub(" hours", "", es$`time point:ch1`)) es <- es[, order(es$time)] exprs(es) <- normalizeBetweenArrays(log2(exprs(es)), method="quantile") es <- es[order(rowMeans(exprs(es)), decreasing=TRUE), ] es <- es[!duplicated(fData(es)$`Gene ID`), ] rownames(es) <- fData(es)$`Gene ID` es <- es[!grepl("///", rownames(es)), ] es <- es[rownames(es) != "", ] fData(es) <- fData(es)[, c("ID", "Gene ID", "Gene symbol")] es <- es[head(order(rowMeans(exprs(es)), decreasing=TRUE), 12000), ] head(exprs(es)) ## ----------------------------------------------------------------------------- library(msigdbr) pathwaysDF <- msigdbr("mouse", category="H") pathways <- split(as.character(pathwaysDF$entrez_gene), pathwaysDF$gs_name) ## ----message=FALSE------------------------------------------------------------ library(fgsea) set.seed(1) gesecaRes <- geseca(pathways, exprs(es), minSize = 15, maxSize = 500) ## ----------------------------------------------------------------------------- head(gesecaRes, 10) ## ----fig.width=10, fig.height=4, out.width="100%"----------------------------- plotCoregulationProfile(pathway=pathways[["HALLMARK_E2F_TARGETS"]], E=exprs(es), titles = es$title, conditions=es$`time point:ch1`) ## ----fig.width=10, fig.height=4, out.width="100%"----------------------------- plotCoregulationProfile(pathway=pathways[["HALLMARK_HYPOXIA"]], E=exprs(es), titles = es$title, conditions=es$`time point:ch1`) ## ----fig.width=10, fig.height=6, out.width="100%"----------------------------- plotGesecaTable(gesecaRes |> head(10), pathways, E=exprs(es), titles = es$title) ## ----------------------------------------------------------------------------- E <- t(base::scale(t(exprs(es)), scale=FALSE)) pcaRev <- prcomp(E, center=FALSE) Ered <- pcaRev$x[, 1:10] dim(Ered) ## ----------------------------------------------------------------------------- set.seed(1) gesecaResRed <- geseca(pathways, Ered, minSize = 15, maxSize = 500, center=FALSE) head(gesecaResRed, 10) ## ----fig.width=4, fig.height=4------------------------------------------------ library(ggplot2) ggplot(data=merge(gesecaRes[, list(pathway, logPvalFull=-log10(pval))], gesecaResRed[, list(pathway, logPvalRed=-log10(pval))])) + geom_point(aes(x=logPvalFull, y=logPvalRed)) + coord_fixed() + theme_classic() ## ----------------------------------------------------------------------------- library(Seurat) ## ----fig.width=8, fig.height=3.5---------------------------------------------- obj <- readRDS(url("https://ctlab.itmo.ru/files/software/fgsea/GSE116240.rds")) obj newIds <- c("0"="Adventitial MF", "3"="Adventitial MF", "5"="Adventitial MF", "1"="Intimal non-foamy MF", "2"="Intimal non-foamy MF", "4"="Intimal foamy MF", "7"="ISG+ MF", "8"="Proliferating cells", "9"="T-cells", "6"="cDC1", "10"="cDC2", "11"="Non-immune cells") obj <- RenameIdents(obj, newIds) DimPlot(obj) + ggplot2::coord_fixed() ## ----------------------------------------------------------------------------- obj <- SCTransform(obj, verbose = FALSE, variable.features.n = 10000) ## ----------------------------------------------------------------------------- length(VariableFeatures(obj)) # make sure it's a full gene universe of 10000 genes obj <- RunPCA(obj, assay = "SCT", verbose = FALSE, rev.pca = TRUE, reduction.name = "pca.rev", reduction.key="PCR_", npcs = 50) E <- obj@reductions$pca.rev@feature.loadings ## ----------------------------------------------------------------------------- library(msigdbr) pathwaysDF <- msigdbr("mouse", category="C2", subcategory = "CP:KEGG") pathways <- split(pathwaysDF$gene_symbol, pathwaysDF$gs_name) ## ----------------------------------------------------------------------------- set.seed(1) gesecaRes <- geseca(pathways, E, minSize = 5, maxSize = 500, center = FALSE, eps=1e-100) head(gesecaRes, 10) ## ----fig.width=12, fig.height=7, out.width="100%"----------------------------- topPathways <- gesecaRes[, pathway] |> head(4) titles <- sub("KEGG_", "", topPathways) ps <- plotCoregulationProfileReduction(pathways[topPathways], obj, title=titles, reduction="tsne") cowplot::plot_grid(plotlist=ps[1:4], ncol=2) ## ----fig.width=5, fig.height=3.5, out.width="50%"----------------------------- plotCoregulationProfileReduction(pathways$KEGG_LYSOSOME, obj, title=sprintf("KEGG_LYSOSOME (pval=%.2g)", gesecaRes[match("KEGG_LYSOSOME", pathway), pval]), reduction="tsne") ## ----message=FALSE------------------------------------------------------------ library(Seurat) obj <- readRDS(url("https://ctlab.itmo.ru/files/software/fgsea/275_T_seurat.rds")) ## ----------------------------------------------------------------------------- obj <- SCTransform(obj, assay = "Spatial", verbose = FALSE, variable.features.n = 10000) obj <- RunPCA(obj, assay = "SCT", verbose = FALSE, rev.pca = TRUE, reduction.name = "pca.rev", reduction.key="PCR_", npcs = 50) E <- obj@reductions$pca.rev@feature.loadings ## ----------------------------------------------------------------------------- library(msigdbr) pathwaysDF <- msigdbr("human", category="H") pathways <- split(pathwaysDF$gene_symbol, pathwaysDF$gs_name) ## ----------------------------------------------------------------------------- set.seed(1) gesecaRes <- geseca(pathways, E, minSize = 15, maxSize = 500, center = FALSE) head(gesecaRes, 10) ## ----fig.width=10, fig.height=7, out.width="100%"----------------------------- topPathways <- gesecaRes[, pathway] |> head(4) titles <- sub("HALLMARK_", "", topPathways) ps <- plotCoregulationProfileSpatial(pathways[topPathways], obj, title=titles) cowplot::plot_grid(plotlist=ps, ncol=2) ## ----fig.width=5, fig.height=3.5, out.width="50%"----------------------------- plotCoregulationProfileSpatial(pathways$HALLMARK_OXIDATIVE_PHOSPHORYLATION, obj, title=sprintf("HALLMARK_OXIDATIVE_PHOSPHORYLATION (pval=%.2g)", gesecaRes[ match("HALLMARK_OXIDATIVE_PHOSPHORYLATION", pathway), pval]))