Gene set co-regulation analysis tutorial
This vignette describes GESECA (gene set co-regulation analysis): a method to identify gene sets that have high gene correlation. We will show how GESECA can be used to find regulated pathways in multi-conditional data, where there is no obvious contrast that can be used to rank genes for GSEA analysis. As examples we will consider a time course microarray experiment and a spatial transcriptomics dataset.
Overiew of GESECA method
GESECA takes as an input:
- E - gene expression matrix, where rows and columns correspond to genes and samples respectively.
- P - list of gene sets (i.e. hallmark gene sets).
Note: genes identifier type should be the same for both elements of P and for row names of matrix E.
By default, GESECA method performs centering for rows of the matrix E. So, after that, the gene values are assumed to have zero mean. Then for each gene set p in P let us introduce the gene set score in the following form:
This score was inspired by the variance of principal components from the principal component analysis (PCA). Therefore, the given score can be viewed in terms of explained variance by the gene set p. Geometrically, this can be considered as an embedding of samples into a one-dimensional space, given by a unit vector in which nonzero positions correspond to genes from gene set p.
In the case of row-centered matrix E the variance of highly correlated genes is summed up to a higher score. While the genes that are not correlated cancel each other and the total gene set variance is low. See the toy example:
Toy example of GESECA score calculation
Another major feature of the proposed score is that it does not require an explicit sample annotation or a contrast. As the result, GESECA can be applied to various types of sequencing technologies: RNA-seq, single-cell sequencing, spatial RNA-seq, etc.
To assess statistical significance for a given gene set p we
calculate an empirical P-value by using gene permutations. The
definition of the P-value is given by the following expression: P(random score ≥ score of p). The estimation
of the given P-value is done by sampling random gene sets with the same
size as p from the row names of matrix E. In practice,
the theoretical P-value can be extremely small, so we use the adaptive
multilevel Markov Chain Monte Carlo scheme, that we used previously in
fgseaMultilevel procedure. For more details, see the preprint.
Analysis of time course data
In the first example we will consider a time course data of Th2 activation from the dataset GSE200250.
First, let prepare the dataset. We load it from Gene Expression Omnibus, apply log and quantile normalization and filter lowly expressed genes.
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))
#> GSM6025497 GSM6025506 GSM6025498 GSM6025507 GSM6025499 GSM6025508
#> 20042 15.95540 15.98219 15.98219 16.00741 15.93798 15.96460
#> 20005 15.93039 15.96460 15.90505 15.91366 15.95540 15.72169
#> 20088 15.93798 15.91054 15.84485 15.87839 15.92647 15.91366
#> 20102 15.88375 15.78037 15.82486 15.87013 15.98219 15.94702
#> 20103 15.84885 15.86085 15.77502 15.84141 15.72169 15.85271
#> 20090 15.91054 15.80253 15.93798 15.74269 15.96460 15.90505
#> GSM6025500 GSM6025509 GSM6025501 GSM6025510 GSM6025502 GSM6025511
#> 20042 15.93039 16.00741 15.94702 16.00741 15.92647 15.98219
#> 20005 15.89942 15.94702 15.95540 15.79744 15.91366 15.80778
#> 20088 15.84485 15.92647 15.86573 15.86573 15.80778 15.88972
#> 20102 15.75168 15.90505 15.79744 15.78575 15.91054 15.89411
#> 20103 15.86085 15.82057 15.91054 15.91981 15.91981 15.70228
#> 20090 15.84885 15.88375 15.76562 15.82057 15.76024 15.82057
#> GSM6025503 GSM6025512 GSM6025504 GSM6025513 GSM6025505 GSM6025514
#> 20042 15.98219 16.00741 15.96460 16.00741 16.00741 16.00741
#> 20005 15.90505 15.84485 15.85590 15.78575 15.83803 15.89411
#> 20088 15.82486 15.83803 15.86573 15.93039 15.85271 15.87839
#> 20102 15.86085 15.75168 15.88972 15.94702 15.91054 15.88972
#> 20103 15.84141 15.93039 15.90505 15.89942 15.90505 15.91054
#> 20090 15.93798 15.83089 15.91366 15.73523 15.91981 15.84885Then we obtain the pathway list. Here we use Hallmarks collection from MSigDB database.
library(msigdbr)
pathwaysDF <- msigdbr("mouse", category="H")
pathways <- split(as.character(pathwaysDF$entrez_gene), pathwaysDF$gs_name)Now we can run GESECA analysis:
The resulting table contain GESECA scores and the corresponding P-values:
head(gesecaRes, 10)
#> pathway pctVar pval padj
#> <char> <num> <num> <num>
#> 1: HALLMARK_E2F_TARGETS 1.5886651 3.726293e-48 1.415991e-46
#> 2: HALLMARK_HYPOXIA 1.1041997 3.509752e-34 6.668528e-33
#> 3: HALLMARK_G2M_CHECKPOINT 1.0281398 1.078730e-32 1.366391e-31
#> 4: HALLMARK_MYC_TARGETS_V1 0.5788220 6.993558e-21 6.643880e-20
#> 5: HALLMARK_GLYCOLYSIS 0.5963190 1.773982e-20 1.348227e-19
#> 6: HALLMARK_MYC_TARGETS_V2 0.7244975 3.370652e-20 2.134746e-19
#> 7: HALLMARK_TNFA_SIGNALING_VIA_NFKB 0.5399796 5.436774e-19 2.951391e-18
#> 8: HALLMARK_MITOTIC_SPINDLE 0.3461497 4.304476e-13 2.044626e-12
#> 9: HALLMARK_INTERFERON_GAMMA_RESPONSE 0.2468646 6.688237e-10 2.541530e-09
#> 10: HALLMARK_IL2_STAT5_SIGNALING 0.2528643 6.688237e-10 2.541530e-09
#> log2err size
#> <num> <int>
#> 1: 1.8030940 181
#> 2: 1.5161076 149
#> 3: 1.4815676 175
#> 4: 1.1778933 186
#> 5: 1.1690700 153
#> 6: 1.1601796 50
#> 7: 1.1239150 156
#> 8: 0.9214260 165
#> 9: 0.8012156 157
#> 10: 0.8012156 177We can plot gene expression profile of HALLMARK_E2F_TARGETS pathway and see that these genes are strongly activated at 24 hours time point:
plotCoregulationProfile(pathway=pathways[["HALLMARK_E2F_TARGETS"]],
E=exprs(es), titles = es$title, conditions=es$`time point:ch1`)
Hypoxia genes have slightly different profile, getting activated around 48 hours:
plotCoregulationProfile(pathway=pathways[["HALLMARK_HYPOXIA"]],
E=exprs(es), titles = es$title, conditions=es$`time point:ch1`)
To get an overview of the top pathway patterns we can use
plotGesecaTable function:
plotGesecaTable(gesecaRes |> head(10), pathways, E=exprs(es), titles = es$title)
#> Warning in get_plot_component(plot, "guide-box"): Multiple components found;
#> returning the first one. To return all, use `return_all = TRUE`.
When the expression matrix contains many samples, a PCA-reduced expression matrix can be used instead of the full matrix to improve the performance. Let reduce the sample space from 18 to 10 dimensions, preserving as much gene variation as possible.
E <- t(base::scale(t(exprs(es)), scale=FALSE))
pcaRev <- prcomp(E, center=FALSE)
Ered <- pcaRev$x[, 1:10]
dim(Ered)
#> [1] 12000 10Now we can run GESECA on the reduced matrix, however we need to disable automatic centering, as we already have done it before the reduction.
set.seed(1)
gesecaResRed <- geseca(pathways, Ered, minSize = 15, maxSize = 500, center=FALSE)
head(gesecaResRed, 10)
#> pathway pctVar pval padj
#> <char> <num> <num> <num>
#> 1: HALLMARK_E2F_TARGETS 1.6274291 1.666912e-47 6.167575e-46
#> 2: HALLMARK_HYPOXIA 1.1242216 4.502652e-33 5.642121e-32
#> 3: HALLMARK_G2M_CHECKPOINT 1.0527938 4.574693e-33 5.642121e-32
#> 4: HALLMARK_MYC_TARGETS_V2 0.7422724 8.826326e-21 6.602863e-20
#> 5: HALLMARK_MYC_TARGETS_V1 0.5923285 8.922788e-21 6.602863e-20
#> 6: HALLMARK_GLYCOLYSIS 0.6088536 5.051598e-20 3.115152e-19
#> 7: HALLMARK_TNFA_SIGNALING_VIA_NFKB 0.5493606 5.551497e-18 2.934363e-17
#> 8: HALLMARK_MITOTIC_SPINDLE 0.3530595 2.868317e-12 1.326596e-11
#> 9: HALLMARK_IL2_STAT5_SIGNALING 0.2584534 6.767387e-10 2.782148e-09
#> 10: HALLMARK_INTERFERON_GAMMA_RESPONSE 0.2501224 1.838503e-09 6.802462e-09
#> log2err size
#> <num> <int>
#> 1: 1.7915725 181
#> 2: 1.4885397 149
#> 3: 1.4885397 175
#> 4: 1.1778933 50
#> 5: 1.1778933 186
#> 6: 1.1512205 153
#> 7: 1.0864405 156
#> 8: 0.8986712 165
#> 9: 0.8012156 177
#> 10: 0.7749390 157The scores and P-values are similar to the ones we obtained for the full matrix.
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()
Analysis of single-cell RNA-seq
Les us load neccesary libraries. We a going to use
Seurat package for working with singe cell data.
library(Seurat)
#> Loading required package: SeuratObject
#> Loading required package: sp
#>
#> Attaching package: 'sp'
#> The following object is masked from 'package:IRanges':
#>
#> %over%
#> 'SeuratObject' was built under R 4.4.0 but the current version is
#> 4.4.2; it is recomended that you reinstall 'SeuratObject' as the ABI
#> for R may have changed
#>
#> Attaching package: 'SeuratObject'
#> The following object is masked from 'package:IRanges':
#>
#> intersect
#> The following object is masked from 'package:S4Vectors':
#>
#> intersect
#> The following object is masked from 'package:BiocGenerics':
#>
#> intersect
#> The following objects are masked from 'package:base':
#>
#> intersect, tAs an example dataset we will use GSE116240 (https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE116240). The dataset features single cell RNA sequencing of aortic CD45+ cells and foam cells from atherosclerotic aorta and is extensively described in the corresponding publication (https://pubmed.ncbi.nlm.nih.gov/30359200/). We also thank the authors for providing the corresponding Seurat object used for the publication.
obj <- readRDS(url("https://alserglab.wustl.edu/files/fgsea/GSE116240.rds"))
obj
#> An object of class Seurat
#> 27998 features across 3781 samples within 1 assay
#> Active assay: RNA (27998 features, 3623 variable features)
#> 2 layers present: counts, data
#> 2 dimensional reductions calculated: pca, tsne
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()
We apply an appropriate normalization (note that we are using 10000 genes, which will be later used as a gene universe for the analysis):
obj <- SCTransform(obj, verbose = FALSE, variable.features.n = 10000)
#> `vst.flavor` is set to 'v2' but could not find glmGamPoi installed.
#> Please install the glmGamPoi package for much faster estimation.
#> --------------------------------------------
#> install.packages('BiocManager')
#> BiocManager::install('glmGamPoi')
#> --------------------------------------------
#> Falling back to native (slower) implementation.
#> Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
#> control$trace > : iteration limit reached
#> Warning in sqrt(1/i): NaNs produced
#> Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
#> control$trace > : iteration limit reached
#> Warning in sqrt(1/i): NaNs produced
#> Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
#> control$trace > : iteration limit reached
#> Warning in sqrt(1/i): NaNs produced
#> Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
#> control$trace > : iteration limit reached
#> Warning in sqrt(1/i): NaNs produced
#> Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
#> control$trace > : iteration limit reached
#> Warning in sqrt(1/i): NaNs produced
#> Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
#> control$trace > : iteration limit reached
#> Warning in sqrt(1/i): NaNs produced
#> Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
#> control$trace > : iteration limit reached
#> Warning in sqrt(1/i): NaNs produced
#> Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
#> control$trace > : iteration limit reached
#> Warning in sqrt(1/i): NaNs produced
#> Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
#> control$trace > : iteration limit reached
#> Warning in sqrt(1/i): NaNs produced
#> Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
#> control$trace > : iteration limit reached
#> Warning in sqrt(1/i): NaNs produced
#> Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
#> control$trace > : iteration limit reached
#> Warning in sqrt(1/i): NaNs produced
#> Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
#> control$trace > : iteration limit reached
#> Warning in sqrt(1/i): NaNs produced
#> Warning in glm.nb(formula = as.formula(new_formula), data = data): alternation
#> limit reached
#> Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
#> control$trace > : iteration limit reached
#> Warning in sqrt(1/i): NaNs produced
#> Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
#> control$trace > : iteration limit reached
#> Warning in sqrt(1/i): NaNs produced
#> Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
#> control$trace > : iteration limit reached
#> Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
#> control$trace > : iteration limit reached
#> Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
#> control$trace > : iteration limit reached
#> Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
#> control$trace > : iteration limit reached
#> Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
#> control$trace > : iteration limit reached
#> Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
#> control$trace > : iteration limit reached
#> Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
#> control$trace > : iteration limit reached
#> Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
#> control$trace > : iteration limit reached
#> Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
#> control$trace > : iteration limit reached
#> Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
#> control$trace > : iteration limit reached
#> Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
#> control$trace > : iteration limit reached
#> Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
#> control$trace > : iteration limit reached
#> Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
#> control$trace > : iteration limit reached
#> Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
#> control$trace > : iteration limit reached
#> Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
#> control$trace > : iteration limit reached
#> Warning in sqrt(1/i): NaNs produced
#> Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
#> control$trace > : iteration limit reached
#> Warning in sqrt(1/i): NaNs produced
#> Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
#> control$trace > : iteration limit reached
#> Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
#> control$trace > : iteration limit reached
#> Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
#> control$trace > : iteration limit reached
#> Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
#> control$trace > : iteration limit reached
#> Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
#> control$trace > : iteration limit reached
#> Warning in sqrt(1/i): NaNs produced
#> Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
#> control$trace > : iteration limit reached
#> Warning in sqrt(1/i): NaNs produced
#> Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
#> control$trace > : iteration limit reached
#> Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
#> control$trace > : iteration limit reached
#> Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
#> control$trace > : iteration limit reached
#> Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
#> control$trace > : iteration limit reached
#> Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
#> control$trace > : iteration limit reached
#> Warning in sqrt(1/i): NaNs produced
#> Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
#> control$trace > : iteration limit reached
#> Warning in sqrt(1/i): NaNs produced
#> Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
#> control$trace > : iteration limit reached
#> Warning in sqrt(1/i): NaNs produced
#> Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
#> control$trace > : iteration limit reached
#> Warning in sqrt(1/i): NaNs produced
#> Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
#> control$trace > : iteration limit reached
#> Warning in sqrt(1/i): NaNs produced
#> Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
#> control$trace > : iteration limit reached
#> Warning in sqrt(1/i): NaNs produced
#> Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
#> control$trace > : iteration limit reached
#> Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
#> control$trace > : iteration limit reached
#> Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
#> control$trace > : iteration limit reached
#> Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
#> control$trace > : iteration limit reached
#> Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
#> control$trace > : iteration limit reached
#> Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
#> control$trace > : iteration limit reached
#> Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
#> control$trace > : iteration limit reached
#> Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
#> control$trace > : iteration limit reached
#> Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
#> control$trace > : iteration limit reached
#> Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
#> control$trace > : iteration limit reached
#> Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
#> control$trace > : iteration limit reached
#> Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
#> control$trace > : iteration limit reached
#> Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
#> control$trace > : iteration limit reached
#> Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
#> control$trace > : iteration limit reached
#> Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
#> control$trace > : iteration limit reached
#> Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
#> control$trace > : iteration limit reached
#> Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
#> control$trace > : iteration limit reached
#> Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
#> control$trace > : iteration limit reached
#> Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
#> control$trace > : iteration limit reached
#> Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
#> control$trace > : iteration limit reached
#> Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
#> control$trace > : iteration limit reached
#> Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
#> control$trace > : iteration limit reached
#> Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
#> control$trace > : iteration limit reached
#> Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
#> control$trace > : iteration limit reached
#> Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
#> control$trace > : iteration limit reached
#> Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
#> control$trace > : iteration limit reached
#> Warning in glm.nb(formula = as.formula(new_formula), data = data): alternation
#> limit reached
#> Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
#> control$trace > : iteration limit reached
#> Warning in sqrt(1/i): NaNs produced
#> Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
#> control$trace > : iteration limit reached
#> Warning in sqrt(1/i): NaNs produced
#> Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
#> control$trace > : iteration limit reached
#> Warning in sqrt(1/i): NaNs produced
#> Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
#> control$trace > : iteration limit reached
#> Warning in sqrt(1/i): NaNs produced
#> Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
#> control$trace > : iteration limit reached
#> Warning in sqrt(1/i): NaNs produced
#> Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
#> control$trace > : iteration limit reached
#> Warning in sqrt(1/i): NaNs produced
#> Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
#> control$trace > : iteration limit reached
#> Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
#> control$trace > : iteration limit reached
#> Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
#> control$trace > : iteration limit reached
#> Warning in sqrt(1/i): NaNs produced
#> Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
#> control$trace > : iteration limit reached
#> Warning in sqrt(1/i): NaNs produced
#> Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
#> control$trace > : iteration limit reached
#> Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
#> control$trace > : iteration limit reached
#> Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
#> control$trace > : iteration limit reached
#> Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
#> control$trace > : iteration limit reached
#> Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
#> control$trace > : iteration limit reached
#> Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
#> control$trace > : iteration limit reached
#> Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
#> control$trace > : iteration limit reached
#> Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
#> control$trace > : iteration limit reached
#> Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
#> control$trace > : iteration limit reached
#> Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
#> control$trace > : iteration limit reached
#> Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
#> control$trace > : iteration limit reached
#> Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
#> control$trace > : iteration limit reached
#> Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
#> control$trace > : iteration limit reached
#> Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
#> control$trace > : iteration limit reached
#> Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
#> control$trace > : iteration limit reached
#> Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
#> control$trace > : iteration limit reached
#> Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
#> control$trace > : iteration limit reached
#> Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
#> control$trace > : iteration limit reached
#> Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
#> control$trace > : iteration limit reached
#> Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
#> control$trace > : iteration limit reached
#> Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
#> control$trace > : iteration limit reached
#> Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
#> control$trace > : iteration limit reached
#> Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
#> control$trace > : iteration limit reached
#> Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
#> control$trace > : iteration limit reached
#> Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
#> control$trace > : iteration limit reached
#> Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
#> control$trace > : iteration limit reached
#> Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
#> control$trace > : iteration limit reached
#> Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
#> control$trace > : iteration limit reached
#> Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
#> control$trace > : iteration limit reached
#> Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
#> control$trace > : iteration limit reached
#> Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
#> control$trace > : iteration limit reached
#> Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
#> control$trace > : iteration limit reached
#> Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
#> control$trace > : iteration limit reached
#> Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
#> control$trace > : iteration limit reached
#> Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
#> control$trace > : iteration limit reached
#> Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
#> control$trace > : iteration limit reached
#> Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
#> control$trace > : iteration limit reached
#> Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
#> control$trace > : iteration limit reached
#> Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
#> control$trace > : iteration limit reached
#> Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
#> control$trace > : iteration limit reached
#> Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
#> control$trace > : iteration limit reached
#> Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
#> control$trace > : iteration limit reached
#> Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
#> control$trace > : iteration limit reached
#> Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
#> control$trace > : iteration limit reached
#> Warning in glm.nb(formula = as.formula(new_formula), data = data): alternation
#> limit reached
#> Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
#> control$trace > : iteration limit reached
#> Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
#> control$trace > : iteration limit reached
#> Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
#> control$trace > : iteration limit reached
#> Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
#> control$trace > : iteration limit reached
#> Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
#> control$trace > : iteration limit reached
#> Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
#> control$trace > : iteration limit reached
#> Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
#> control$trace > : iteration limit reached
#> Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
#> control$trace > : iteration limit reached
#> Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
#> control$trace > : iteration limit reached
#> Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
#> control$trace > : iteration limit reached
#> Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
#> control$trace > : iteration limit reached
#> Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
#> control$trace > : iteration limit reached
#> Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
#> control$trace > : iteration limit reached
#> Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
#> control$trace > : iteration limit reached
#> Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
#> control$trace > : iteration limit reached
#> Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
#> control$trace > : iteration limit reached
#> Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
#> control$trace > : iteration limit reached
#> Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
#> control$trace > : iteration limit reached
#> Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
#> control$trace > : iteration limit reached
#> Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
#> control$trace > : iteration limit reached
#> Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
#> control$trace > : iteration limit reached
#> Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
#> control$trace > : iteration limit reached
#> Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
#> control$trace > : iteration limit reached
#> Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
#> control$trace > : iteration limit reached
#> Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
#> control$trace > : iteration limit reached
#> Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
#> control$trace > : iteration limit reached
#> Warning in glm.nb(formula = as.formula(new_formula), data = data): alternation
#> limit reached
#> Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
#> control$trace > : iteration limit reached
#> Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
#> control$trace > : iteration limit reached
#> Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
#> control$trace > : iteration limit reached
#> Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
#> control$trace > : iteration limit reached
#> Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
#> control$trace > : iteration limit reached
#> Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
#> control$trace > : iteration limit reached
#> Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
#> control$trace > : iteration limit reached
#> Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
#> control$trace > : iteration limit reached
#> Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
#> control$trace > : iteration limit reached
#> Warning in sqrt(1/i): NaNs produced
#> Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
#> control$trace > : iteration limit reached
#> Warning in sqrt(1/i): NaNs produced
#> Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
#> control$trace > : iteration limit reached
#> Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
#> control$trace > : iteration limit reached
#> Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
#> control$trace > : iteration limit reached
#> Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
#> control$trace > : iteration limit reached
#> Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
#> control$trace > : iteration limit reached
#> Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
#> control$trace > : iteration limit reached
#> Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
#> control$trace > : iteration limit reached
#> Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
#> control$trace > : iteration limit reached
#> Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
#> control$trace > : iteration limit reached
#> Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
#> control$trace > : iteration limit reached
#> Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
#> control$trace > : iteration limit reached
#> Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
#> control$trace > : iteration limit reached
#> Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
#> control$trace > : iteration limit reached
#> Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
#> control$trace > : iteration limit reached
#> Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
#> control$trace > : iteration limit reached
#> Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
#> control$trace > : iteration limit reached
#> Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
#> control$trace > : iteration limit reached
#> Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
#> control$trace > : iteration limit reached
#> Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
#> control$trace > : iteration limit reached
#> Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
#> control$trace > : iteration limit reached
#> Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
#> control$trace > : iteration limit reached
#> Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
#> control$trace > : iteration limit reached
#> Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
#> control$trace > : iteration limit reached
#> Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
#> control$trace > : iteration limit reached
#> Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
#> control$trace > : iteration limit reached
#> Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
#> control$trace > : iteration limit reached
#> Warning in glm.nb(formula = as.formula(new_formula), data = data): alternation
#> limit reachedTo speed up the analysis, instead of using the full transformed gene
expression matrix, we will consider only its first principal components.
Note that a “reverse” PCA should be done: the principal components
should correspond to linear combinations of the cells, not linear
combinations of the genes as in “normal” PCA. By default
SCTransform returns centered gene expression, so we can run
PCA directly.
length(VariableFeatures(obj)) # make sure it's a full gene universe of 10000 genes
#> [1] 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.loadingsFollowing the authors we are going to use KEGG pathway collection.
library(msigdbr)
pathwaysDF <- msigdbr("mouse", category="C2", subcategory = "CP:KEGG")
pathways <- split(pathwaysDF$gene_symbol, pathwaysDF$gs_name)Now we can run the analysis (we set center=FALSE because
we use the reduced matrix):
set.seed(1)
gesecaRes <- geseca(pathways, E, minSize = 5, maxSize = 500, center = FALSE, eps=1e-100)
head(gesecaRes, 10)
#> pathway pctVar pval
#> <char> <num> <num>
#> 1: KEGG_RIBOSOME 1.1840580 1.508018e-41
#> 2: KEGG_LEISHMANIA_INFECTION 0.4326062 3.910683e-19
#> 3: KEGG_NOD_LIKE_RECEPTOR_SIGNALING_PATHWAY 0.3995194 2.762109e-18
#> 4: KEGG_ASTHMA 1.0006135 3.020242e-18
#> 5: KEGG_GRAFT_VERSUS_HOST_DISEASE 0.6329171 4.517426e-18
#> 6: KEGG_ANTIGEN_PROCESSING_AND_PRESENTATION 0.4169731 3.583911e-17
#> 7: KEGG_TYPE_I_DIABETES_MELLITUS 0.5785945 7.540207e-17
#> 8: KEGG_ALLOGRAFT_REJECTION 0.5492550 1.820777e-16
#> 9: KEGG_LYSOSOME 0.2791601 8.867462e-16
#> 10: KEGG_CYTOKINE_CYTOKINE_RECEPTOR_INTERACTION 0.2401663 6.083392e-15
#> padj log2err size
#> <num> <num> <int>
#> 1: 1.100853e-39 1.6719967 61
#> 2: 1.427399e-17 1.1239150 57
#> 3: 5.511942e-17 1.0959293 50
#> 4: 5.511942e-17 1.0959293 12
#> 5: 6.595442e-17 1.0959293 26
#> 6: 4.360425e-16 1.0672100 49
#> 7: 7.863358e-16 1.0574636 26
#> 8: 1.661459e-15 1.0376962 24
#> 9: 7.192497e-15 1.0175448 102
#> 10: 4.440876e-14 0.9865463 133Now we can plot profiles of the top pathways (we need to specify the reduction as we are using the one from the publication):
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)
We can see that inflammatory pathways (e.g. KEGG_LEISHMANIA_INFECTION) are more associated with the non-foamy intimal macrophages, which was one of the main points of the Kim et al. Another pathway highlighted by the authors, KEGG_LYSOSOME, is specific to intimal foamy macrophages:
plotCoregulationProfileReduction(pathways$KEGG_LYSOSOME,
obj,
title=sprintf("KEGG_LYSOSOME (pval=%.2g)",
gesecaRes[match("KEGG_LYSOSOME", pathway), pval]),
reduction="tsne")
Analysis of spatial RNA-seq
Similarly to single cell RNA-seq, GESECA can be applied to spatial transcriptomics profiling. As an example we will use glioblastoma sample from Ravi et al (https://pubmed.ncbi.nlm.nih.gov/35700707/).
As for scRNA-seq, we apply normalization for 10000 genes and run a do a PCA reduction.
obj <- SCTransform(obj, assay = "Spatial", verbose = FALSE, variable.features.n = 10000)
#> `vst.flavor` is set to 'v2' but could not find glmGamPoi installed.
#> Please install the glmGamPoi package for much faster estimation.
#> --------------------------------------------
#> install.packages('BiocManager')
#> BiocManager::install('glmGamPoi')
#> --------------------------------------------
#> Falling back to native (slower) implementation.
#> Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
#> control$trace > : iteration limit reached
#> Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
#> control$trace > : iteration limit reached
#> Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
#> control$trace > : iteration limit reached
#> Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
#> control$trace > : iteration limit reached
#> Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
#> control$trace > : iteration limit reached
#> Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
#> control$trace > : iteration limit reached
#> Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
#> control$trace > : iteration limit reached
#> Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
#> control$trace > : iteration limit reached
#> Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
#> control$trace > : iteration limit reached
#> Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
#> control$trace > : iteration limit reached
#> Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
#> control$trace > : iteration limit reached
#> Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
#> control$trace > : iteration limit reached
#> Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
#> control$trace > : iteration limit reached
#> Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
#> control$trace > : iteration limit reached
#> Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
#> control$trace > : iteration limit reached
#> Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
#> control$trace > : iteration limit reached
#> Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
#> control$trace > : iteration limit reached
#> Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
#> control$trace > : iteration limit reached
#> Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
#> control$trace > : iteration limit reached
#> Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
#> control$trace > : iteration limit reached
#> Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
#> control$trace > : iteration limit reached
#> Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
#> control$trace > : iteration limit reached
#> Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
#> control$trace > : iteration limit reached
#> Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
#> control$trace > : iteration limit reached
#> Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
#> control$trace > : iteration limit reached
#> Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
#> control$trace > : iteration limit reached
#> Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
#> control$trace > : iteration limit reached
#> Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
#> control$trace > : iteration limit reached
#> Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
#> control$trace > : iteration limit reached
#> Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
#> control$trace > : iteration limit reached
#> Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
#> control$trace > : iteration limit reached
#> Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
#> control$trace > : iteration limit reached
#> Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
#> control$trace > : iteration limit reached
#> Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
#> control$trace > : iteration limit reached
#> Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
#> control$trace > : iteration limit reached
#> Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
#> control$trace > : iteration limit reached
#> Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
#> control$trace > : iteration limit reached
#> Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
#> control$trace > : iteration limit reached
#> Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
#> control$trace > : iteration limit reached
#> Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
#> control$trace > : iteration limit reached
#> Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
#> control$trace > : iteration limit reached
#> Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
#> control$trace > : iteration limit reached
#> Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
#> control$trace > : iteration limit reached
#> Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
#> control$trace > : iteration limit reached
#> Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
#> control$trace > : iteration limit reached
#> Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
#> control$trace > : iteration limit reached
#> Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
#> control$trace > : iteration limit reached
#> Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
#> control$trace > : iteration limit reached
#> Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
#> control$trace > : iteration limit reached
#> Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
#> control$trace > : iteration limit reached
#> Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
#> control$trace > : iteration limit reached
#> Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
#> control$trace > : iteration limit reached
#> Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
#> control$trace > : iteration limit reached
#> Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
#> control$trace > : iteration limit reached
#> Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
#> control$trace > : iteration limit reached
#> Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
#> control$trace > : iteration limit reached
#> Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
#> control$trace > : iteration limit reached
#> Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
#> control$trace > : iteration limit reached
#> Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
#> control$trace > : iteration limit reached
#> Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
#> control$trace > : iteration limit reached
#> Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
#> control$trace > : iteration limit reached
#> Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
#> control$trace > : iteration limit reached
#> Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
#> control$trace > : iteration limit reached
#> Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
#> control$trace > : iteration limit reached
#> Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
#> control$trace > : iteration limit reached
#> Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
#> control$trace > : iteration limit reached
#> Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
#> control$trace > : iteration limit reached
#> Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
#> control$trace > : iteration limit reached
#> Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
#> control$trace > : iteration limit reached
#> Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
#> control$trace > : iteration limit reached
#> Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
#> control$trace > : iteration limit reached
#> Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
#> control$trace > : iteration limit reached
#> Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
#> control$trace > : iteration limit reached
#> Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
#> control$trace > : iteration limit reached
#> Warning in glm.nb(formula = as.formula(new_formula), data = data): alternation
#> limit reached
#> Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
#> control$trace > : iteration limit reached
#> Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
#> control$trace > : iteration limit reached
#> Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
#> control$trace > : iteration limit reached
#> Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
#> control$trace > : iteration limit reached
#> Warning in glm.nb(formula = as.formula(new_formula), data = data): alternation
#> limit reached
#> Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
#> control$trace > : iteration limit reached
#> Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
#> control$trace > : iteration limit reached
#> Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
#> control$trace > : iteration limit reached
#> Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
#> control$trace > : iteration limit reached
#> Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
#> control$trace > : iteration limit reached
#> Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
#> control$trace > : iteration limit reached
#> Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
#> control$trace > : iteration limit reached
#> Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
#> control$trace > : iteration limit reached
#> Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
#> control$trace > : iteration limit reached
#> Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
#> control$trace > : iteration limit reached
#> Warning in glm.nb(formula = as.formula(new_formula), data = data): alternation
#> limit reached
#> Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
#> control$trace > : iteration limit reached
#> Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
#> control$trace > : iteration limit reached
#> Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
#> control$trace > : iteration limit reached
#> Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
#> control$trace > : iteration limit reached
#> Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
#> control$trace > : iteration limit reached
#> Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
#> control$trace > : iteration limit reached
#> Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
#> control$trace > : iteration limit reached
#> Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
#> control$trace > : iteration limit reached
#> Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
#> control$trace > : iteration limit reached
#> Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
#> control$trace > : iteration limit reached
#> Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
#> control$trace > : iteration limit reached
#> Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
#> control$trace > : iteration limit reached
#> Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
#> control$trace > : iteration limit reached
#> Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
#> control$trace > : iteration limit reached
#> Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
#> control$trace > : iteration limit reached
#> Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
#> control$trace > : iteration limit reached
#> Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
#> control$trace > : iteration limit reached
#> Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
#> control$trace > : iteration limit reached
#> Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
#> control$trace > : iteration limit reached
#> Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
#> control$trace > : iteration limit reached
#> Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
#> control$trace > : iteration limit reached
#> Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
#> control$trace > : iteration limit reached
#> Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
#> control$trace > : iteration limit reached
#> Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
#> control$trace > : iteration limit reached
#> Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
#> control$trace > : iteration limit reached
#> Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
#> control$trace > : iteration limit reached
#> Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
#> control$trace > : iteration limit reached
#> Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
#> control$trace > : iteration limit reached
#> Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
#> control$trace > : iteration limit reached
#> Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
#> control$trace > : iteration limit reached
#> Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
#> control$trace > : iteration limit reached
#> Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
#> control$trace > : iteration limit reached
#> Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
#> control$trace > : iteration limit reached
#> Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
#> control$trace > : iteration limit reached
#> Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
#> control$trace > : iteration limit reached
#> Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
#> control$trace > : iteration limit reached
#> Warning in glm.nb(formula = as.formula(new_formula), data = data): alternation
#> limit reached
#> Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
#> control$trace > : iteration limit reached
#> Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
#> control$trace > : iteration limit reached
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.loadingsWe will use HALLMARK pathways as the gene set collection.
library(msigdbr)
pathwaysDF <- msigdbr("human", category="H")
pathways <- split(pathwaysDF$gene_symbol, pathwaysDF$gs_name)Now we can run the analysis (remember, we set
center=FALSE because we use the reduced matrix):
set.seed(1)
gesecaRes <- geseca(pathways, E, minSize = 15, maxSize = 500, center = FALSE)
#> Warning in geseca(pathways, E, minSize = 15, maxSize = 500, center = FALSE):
#> For some pathways, in reality P-values are less than 1e-50. You can set the
#> `eps` argument to zero for better estimation.
head(gesecaRes, 10)
#> pathway pctVar pval
#> <char> <num> <num>
#> 1: HALLMARK_HYPOXIA 1.2775395 1.000000e-50
#> 2: HALLMARK_TNFA_SIGNALING_VIA_NFKB 0.5286872 1.013678e-26
#> 3: HALLMARK_EPITHELIAL_MESENCHYMAL_TRANSITION 0.4610028 7.936135e-24
#> 4: HALLMARK_INTERFERON_GAMMA_RESPONSE 0.4533018 4.841692e-23
#> 5: HALLMARK_INTERFERON_ALPHA_RESPONSE 0.4756488 1.146510e-21
#> 6: HALLMARK_OXIDATIVE_PHOSPHORYLATION 0.2678858 5.120417e-14
#> 7: HALLMARK_MTORC1_SIGNALING 0.2338037 4.271464e-12
#> 8: HALLMARK_GLYCOLYSIS 0.1878535 1.566967e-09
#> 9: HALLMARK_ALLOGRAFT_REJECTION 0.1649848 7.426823e-09
#> 10: HALLMARK_INFLAMMATORY_RESPONSE 0.1672958 1.466740e-08
#> padj log2err size
#> <num> <num> <int>
#> 1: 2.500000e-49 NA 150
#> 2: 1.267097e-25 1.3344975 150
#> 3: 6.613446e-23 1.2627399 146
#> 4: 3.026058e-22 1.2378967 150
#> 5: 5.732552e-21 1.2039752 77
#> 6: 2.133507e-13 0.9545416 184
#> 7: 1.525523e-11 0.8870750 180
#> 8: 4.896772e-09 0.7881868 144
#> 9: 2.063006e-08 0.7477397 98
#> 10: 3.666850e-08 0.7337620 112Finally, let us plot spatial expression of the top four pathways:
topPathways <- gesecaRes[, pathway] |> head(4)
titles <- sub("HALLMARK_", "", topPathways)
pt.size.factor <- 1.6
# Starting from Seurat version 5.1.0, the scaling method for spatial plots was changed.
# As mentioned here: https://github.com/satijalab/seurat/issues/9049
# This code provides a workaround to adapt to the new scaling behavior.
if (packageVersion("Seurat") >= "5.1.0"){
sfactors <- ScaleFactors(obj@images$slice1)
pt.size.factor <- 1.5 * sfactors$fiducial / sfactors$hires
}
ps <- plotCoregulationProfileSpatial(pathways[topPathways], obj,
title=titles,
pt.size.factor=pt.size.factor)
cowplot::plot_grid(plotlist=ps, ncol=2)
Consistent with the Ravi et al, we see a distinct hypoxic region (defined by HALLMARK_HYPOXIA) and a reactive immune region (defined by HALLMARK_INTERFERON_GAMMA_RESPONSE). Further, we can explore behavior of HALLMARK_OXIDATIVE_PHOSPHORYLATION pathway to see that oxidative metabolism is more characteristic to the “normal” tissue region:
plotCoregulationProfileSpatial(pathways$HALLMARK_OXIDATIVE_PHOSPHORYLATION,
obj,
pt.size.factor=pt.size.factor,
title=sprintf("HALLMARK_OXIDATIVE_PHOSPHORYLATION\n(pval=%.2g)",
gesecaRes[
match("HALLMARK_OXIDATIVE_PHOSPHORYLATION", pathway),
pval]))
#> [[1]]
Session info
sessionInfo()
#> R version 4.4.2 (2024-10-31)
#> Platform: x86_64-pc-linux-gnu
#> Running under: Ubuntu 24.04.1 LTS
#>
#> Matrix products: default
#> BLAS: /usr/lib/x86_64-linux-gnu/openblas-pthread/libblas.so.3
#> LAPACK: /usr/lib/x86_64-linux-gnu/openblas-pthread/libopenblasp-r0.3.26.so; LAPACK version 3.12.0
#>
#> locale:
#> [1] LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C
#> [3] LC_TIME=en_US.UTF-8 LC_COLLATE=C
#> [5] LC_MONETARY=en_US.UTF-8 LC_MESSAGES=en_US.UTF-8
#> [7] LC_PAPER=en_US.UTF-8 LC_NAME=C
#> [9] LC_ADDRESS=C LC_TELEPHONE=C
#> [11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C
#>
#> time zone: Etc/UTC
#> tzcode source: system (glibc)
#>
#> attached base packages:
#> [1] stats4 stats graphics grDevices utils datasets methods
#> [8] base
#>
#> other attached packages:
#> [1] Seurat_5.1.0 SeuratObject_5.0.2 sp_2.1-4
#> [4] msigdbr_7.5.1 limma_3.63.2 GEOquery_2.75.0
#> [7] org.Mm.eg.db_3.20.0 AnnotationDbi_1.69.0 IRanges_2.41.2
#> [10] S4Vectors_0.45.2 Biobase_2.67.0 BiocGenerics_0.53.3
#> [13] generics_0.1.3 BiocParallel_1.41.0 ggplot2_3.5.1
#> [16] data.table_1.16.4 fgsea_1.33.2 rmarkdown_2.29
#>
#> loaded via a namespace (and not attached):
#> [1] RcppAnnoy_0.0.22 splines_4.4.2
#> [3] later_1.4.1 tibble_3.2.1
#> [5] R.oo_1.27.0 polyclip_1.10-7
#> [7] XML_3.99-0.17 fastDummies_1.7.4
#> [9] lifecycle_1.0.4 httr2_1.0.7
#> [11] globals_0.16.3 lattice_0.22-6
#> [13] MASS_7.3-61 magrittr_2.0.3
#> [15] plotly_4.10.4 sass_0.4.9
#> [17] jquerylib_0.1.4 yaml_2.3.10
#> [19] httpuv_1.6.15 sctransform_0.4.1
#> [21] spam_2.11-0 spatstat.sparse_3.1-0
#> [23] reticulate_1.40.0 cowplot_1.1.3
#> [25] pbapply_1.7-2 DBI_1.2.3
#> [27] buildtools_1.0.0 RColorBrewer_1.1-3
#> [29] abind_1.4-8 zlibbioc_1.52.0
#> [31] Rtsne_0.17 GenomicRanges_1.59.1
#> [33] purrr_1.0.2 R.utils_2.12.3
#> [35] rappdirs_0.3.3 GenomeInfoDbData_1.2.13
#> [37] ggrepel_0.9.6 irlba_2.3.5.1
#> [39] listenv_0.9.1 spatstat.utils_3.1-1
#> [41] maketools_1.3.1 rentrez_1.2.3
#> [43] reactome.db_1.89.0 goftest_1.2-3
#> [45] RSpectra_0.16-2 spatstat.random_3.3-2
#> [47] fitdistrplus_1.2-1 parallelly_1.41.0
#> [49] leiden_0.4.3.1 codetools_0.2-20
#> [51] DelayedArray_0.33.3 xml2_1.3.6
#> [53] tidyselect_1.2.1 UCSC.utils_1.3.0
#> [55] farver_2.1.2 matrixStats_1.4.1
#> [57] spatstat.explore_3.3-3 jsonlite_1.8.9
#> [59] progressr_0.15.1 ggridges_0.5.6
#> [61] survival_3.8-3 tools_4.4.2
#> [63] ica_1.0-3 Rcpp_1.0.13-1
#> [65] glue_1.8.0 gridExtra_2.3
#> [67] SparseArray_1.7.2 xfun_0.49
#> [69] MatrixGenerics_1.19.0 GenomeInfoDb_1.43.2
#> [71] dplyr_1.1.4 withr_3.0.2
#> [73] fastmap_1.2.0 digest_0.6.37
#> [75] R6_2.5.1 mime_0.12
#> [77] colorspace_2.1-1 scattermore_1.2
#> [79] tensor_1.5 spatstat.data_3.1-4
#> [81] RSQLite_2.3.9 R.methodsS3_1.8.2
#> [83] tidyr_1.3.1 httr_1.4.7
#> [85] htmlwidgets_1.6.4 S4Arrays_1.7.1
#> [87] uwot_0.2.2 pkgconfig_2.0.3
#> [89] gtable_0.3.6 blob_1.2.4
#> [91] lmtest_0.9-40 XVector_0.47.0
#> [93] sys_3.4.3 htmltools_0.5.8.1
#> [95] dotCall64_1.2 scales_1.3.0
#> [97] png_0.1-8 spatstat.univar_3.1-1
#> [99] knitr_1.49 tzdb_0.4.0
#> [101] reshape2_1.4.4 nlme_3.1-166
#> [103] curl_6.0.1 cachem_1.1.0
#> [105] zoo_1.8-12 stringr_1.5.1
#> [107] KernSmooth_2.23-24 parallel_4.4.2
#> [109] miniUI_0.1.1.1 pillar_1.10.0
#> [111] grid_4.4.2 vctrs_0.6.5
#> [113] RANN_2.6.2 promises_1.3.2
#> [115] xtable_1.8-4 cluster_2.1.8
#> [117] evaluate_1.0.1 readr_2.1.5
#> [119] cli_3.6.3 compiler_4.4.2
#> [121] rlang_1.1.4 crayon_1.5.3
#> [123] future.apply_1.11.3 labeling_0.4.3
#> [125] plyr_1.8.9 stringi_1.8.4
#> [127] deldir_2.0-4 viridisLite_0.4.2
#> [129] babelgene_22.9 munsell_0.5.1
#> [131] Biostrings_2.75.3 lazyeval_0.2.2
#> [133] spatstat.geom_3.3-4 Matrix_1.7-1
#> [135] RcppHNSW_0.6.0 hms_1.1.3
#> [137] patchwork_1.3.0 bit64_4.5.2
#> [139] future_1.34.0 KEGGREST_1.47.0
#> [141] statmod_1.5.0 shiny_1.10.0
#> [143] SummarizedExperiment_1.37.0 ROCR_1.0-11
#> [145] igraph_2.1.2 memoise_2.0.1
#> [147] bslib_0.8.0 fastmatch_1.1-4
#> [149] bit_4.5.0.1