gsea_presto

Downstream analysis of scRNAseq data

Gene Set enrichment (GSEA) analysis

library(Seurat)
library(devtools)
options(buildtools.check = function(action) TRUE )
install_github('immunogenomics/presto')

BiocManager::install('githubinstall')
library(githubinstall)
githubinstall('immunogenomics/presto')

devtools::install_local("idmap2.zip")

library(tidyverse)
library(presto)

setwd("D:/Data")
load("seu.rda")

pbmc<- seu

nc<-length(unique(pbmc@meta.data $seurat_clusters)) clusters<-as.vector(unique(pbmc@meta.data$ seurat_clusters))

pbmc.genes <- wilcoxauc(pbmc, 'seurat_clusters')
head(pbmc.genes)

we have all the genes for each cluster

dplyr::count(pbmc.genes, group)

BiocManager::install(c("msigdbr", "fgsea"))

library(msigdbr)
library(fgsea)
library(dplyr)
library(ggplot2)

msigdbr_show_species()
m_df<- msigdbr(species = "Homo sapiens", category = "H")
head(m_df)
fgsea_sets<- m_df %>% split(x = . $gene_symbol, f = .$ gs_name)

p<-list()
nes<-matrix(0,50,nc)
p.me<-matrix(0,50,nc)

colnames(nes)<-clusters
colnames(p.me)<-clusters

for (i in 1:nc){

ci<-clusters[i]
print (ci)

select only the feature and auc columns for fgsea, which statistics to use is an open question

cluster0.genes<- pbmc.genes %>%
dplyr::filter(group == ci) %>%
arrange(desc(auc)) %>%
dplyr::select(feature, auc)
ranks<- deframe(cluster0.genes)
head(ranks)
fgseaRes<- fgsea(fgsea_sets, stats = ranks, nperm = 1000)

nes[,i]<-fgseaRes $NES rownames(nes)<-fgseaRes$ pathway
p.me[,i]<-fgseaRes $padj rownames(p.me)<-fgseaRes$ pathway

fgseaResTidy <- fgseaRes %>%
as_tibble() %>%
arrange(desc(NES))
fgseaResTidy %>%
dplyr::select(-leadingEdge, -ES, -nMoreExtreme) %>%
arrange(padj) %>%
head()

```

plot a barplot for with the normalized Enrichment score

}

p.me<--log10(p.me)

pdf("presto_hallmarker_all.pdf")
breaks<-seq(-10,10,1)
legend_breaks<-seq(-10,10,2)

pheatmap(nes, scale="none",#fontsize_row=8,
cluster_rows=TRUE, cluster_cols=TRUE,

     color = c(colorRampPalette(colors=c("navy","white"))(10),colorRampPalette(colors=c("white","firebrick3"))(10)),
     legend_breaks=legend_breaks,
     breaks=breaks,
     fontsize_row=7,
     fontsize_legend=1.0)

pheatmap(nes, scale="row",#fontsize_row=8,
cluster_rows=TRUE, cluster_cols=TRUE,

     color = c(colorRampPalette(colors=c("navy","white"))(5),colorRampPalette(colors=c("white","firebrick3"))(5)),
     
     fontsize_row=7,
     fontsize_legend=1.0)

breaks<-c(0,1.3,2,3)
legend_breaks<-breaks

pheatmap(p.me, scale="none",#fontsize_row=8,
cluster_rows=TRUE, cluster_cols=TRUE,

     clustering_distance_rows = "euclidean",
     clustering_distance_cols = "euclidean",
     color = c(colorRampPalette(colors=c("white","firebrick3"))(3)),
     legend_breaks=legend_breaks,
     breaks=breaks,
     fontsize_row=7,
     fontsize_legend=1.0)

dev.off()