單細胞GSEA分析需要的文件有兩個：
1. 單細胞基因表達變化數(shù)據(jù)（兩列街佑，一列是geneid/gene symbol波闹，一列是logFC）（文件格式：.rnk）
2. 目標基因集（文件格式：.gmx或.gmt）一行是一個term/基因集浙滤，第一列是基因集的名稱，后面是基因分類等，再后面是geneid/gene symbol（要和單細胞矩陣的對應(yīng)，否則不能識別）

1. 單細胞基因表達變化數(shù)據(jù)的準備

• 1.1 下載單細胞數(shù)據(jù)（pbmc5k）

wget http://cf.10xgenomics.com/samples/cell-exp/3.0.2/5k_pbmc_v3/5k_pbmc_v3_filtered_feature_bc_matrix.tar.gz
tar xvzf 5k_pbmc_v3_filtered_feature_bc_matrix.tar.gz

• 1.2 參考Seurat標準流程對數(shù)據(jù)進行預(yù)處理得到分群

library(tidyverse)
library(Seurat)
# Load the PBMC dataset
pbmc.data <- Read10X(data.dir = "filtered_feature_bc_matrix/")
# Initialize the Seurat object with the raw (non-normalized data).
pbmc <- CreateSeuratObject(counts = pbmc.data, project = "pbmc5k", min.cells = 3, min.features = 200)
#質(zhì)控
pbmc[["percent.mt"]] <- PercentageFeatureSet(pbmc, pattern = "^MT-")
#取子集
pbmc <- subset(pbmc, subset = nFeature_RNA > 200 & nFeature_RNA < 5000 & percent.mt < 25)
#標準化和歸一化
pbmc <- NormalizeData(pbmc, normalization.method = "LogNormalize", scale.factor = 10000)
pbmc <- FindVariableFeatures(pbmc, selection.method = "vst", nfeatures = 2000)
all.genes <- rownames(pbmc)
pbmc <- ScaleData(pbmc, vars.to.regress = "percent.mt")
#PCA
pbmc <- RunPCA(pbmc, features = VariableFeatures(object = pbmc), verbose = FALSE)
#聚類
pbmc <- FindNeighbors(pbmc, dims = 1:20)
pbmc <- FindClusters(pbmc, resolution = 0.5)
#TSNE和UMAP降維
pbmc <- RunUMAP(pbmc, dims = 1:20)
pbmc<- RunTSNE(pbmc, dims = 1:20)
DimPlot(pbmc, reduction = "umap", label = TRUE)
#查找marker基因
pbmc.markers <- FindAllMarkers(pbmc, only.pos = TRUE, min.pct = 0.25, logfc.threshold = 0.25)
saveRDS(pbmc, "pbmc_5k_v3.rds")

分為13個亞群

• 1.3 使用wilcoxauc()計算每個cluster的差異基因

pbmc<- readRDS("pbmc_5k_v3.rds")
library(presto)
pbmc.genes <- wilcoxauc(pbmc, 'seurat_clusters')
head(pbmc.genes)
#     feature group     avgExpr         logFC statistic       auc
# 1 AL627309.1     0 0.004801016 -0.0048337160   2113594 0.4966881
# 2 AL627309.3     0 0.000000000 -0.0004923881   2126401 0.4996978
# 3 AL669831.5     0 0.058225088 -0.0101451241   2081447 0.4891337
# 4     FAM87B     0 0.000000000 -0.0017590952   2121900 0.4986401
# 5  LINC00115     0 0.033568380 -0.0054911605   2090280 0.4912094
# 6     FAM41C     0 0.015927696 -0.0162813395   2077756 0.4882664
#          pval         padj    pct_in    pct_out
# 1 0.0451711892 0.0817186691 0.5443235 1.20882442
# 2 0.3781096732 0.4754137751 0.0000000 0.06044122
# 3 0.0182067030 0.0364774662 6.8429238 9.21728619
# 4 0.0612494406 0.1063124181 0.0000000 0.27198549
# 5 0.0148011270 0.0301102065 3.7325039 5.59081293
# 6 0.0001622764 0.0004571717 2.0217729 4.38198852
dplyr::count(pbmc.genes, group)# 查看每個cluster中有多少基因（與矩陣的基因數(shù)一致）
#    group     n
# 1      0 18791
# 2      1 18791
# 3     10 18791
# 4     11 18791
# 5     12 18791
# 6      2 18791
# 7      3 18791
# 8      4 18791
# 9      5 18791
# 10     6 18791
# 11     7 18791
# 12     8 18791
# 13     9 18791

• 1.4 選取group0的差異基因進行GSEA分析件舵，將cluster0的 差異基因整理為GSEA分析需要的數(shù)據(jù)格式（也可輸入注釋好的細胞群的差異基因）

# 僅選擇fgsea的feature和auc列
cluster0.genes<- pbmc.genes %>% dplyr::filter(group == "0") %>% arrange(desc(auc)) %>% dplyr::select(feature, auc)
ranks<- deframe(cluster0.genes)
head(ranks)
# RPL30     RPS3A     RPS13    RPL35A     RPL32     RPL34 
# 0.9376551 0.9358956 0.9250745 0.9196189 0.9185743 0.9132980 

2. 選擇自己數(shù)據(jù)的物種以及要做的GSEA的數(shù)據(jù)庫類型卸察，準備目標基因集的輸入文件

為了進行基因集富集分析脯厨，首先需要注釋基因集。一個重要的來源是Broad Institute開發(fā)的MsigDB坑质。
網(wǎng)址：https://www.gsea-msigdb.org/gsea/msigdb/

library(msigdbr)
library(fgsea)
library(dplyr)
library(ggplot2)
##查看物種的數(shù)據(jù) msigdbr_show_species()合武；#我們使用C7免疫基因集
m_df<- msigdbr(species = "Homo sapiens", category = "C7")
head(m_df)

##將m_df的基因與通路取出并改成一個通路對應(yīng)相應(yīng)基因的格式
fgsea_sets<- m_df %>% split(x = .$gene_symbol, f = .$gs_name)
#以gs_name為factor對gene_symbol進行分類临梗，統(tǒng)計落在每個gs_name中的gene_symbol的個數(shù)，并生成list稼跳。
summary(fgsea_sets)

3. 使用fgsea進行基因集富集并繪圖

• 3.1 富集分析

fgseaRes<- fgsea(fgsea_sets, stats = ranks, nperm = 1000)
#nperm設(shè)置的是permutation次數(shù)

• 3.2 整理數(shù)據(jù)：

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

• 3.3 繪圖

應(yīng)用標準化富集分數(shù)繪制barplot

# 顯示top20信號通路
ggplot(fgseaResTidy %>% filter(padj < 0.008) %>% head(n= 20), aes(reorder(pathway, NES), NES)) +
  geom_col(aes(fill= NES < 7.5)) +
  coord_flip() +
  labs(x="Pathway", y="Normalized Enrichment Score",
       title="Hallmark pathways NES from GSEA") +
  theme_minimal() ####以7.5進行繪圖填色

GSEA圖

plotEnrichment(fgsea_sets[["GSE10325_CD4_TCELL_VS_MYELOID_UP"]],
               ranks) + labs(title="GSE10325 CD4 TCELL VS MYELOID UP")

上圖中盟庞，最上面的綠線是遍歷排好序的基因列表以計算ES值的過程。遍歷基因集中的基因時汤善，當(dāng)基因出現(xiàn)在目標基因（橫軸的基因）中則加分什猖，反之減分。加減分值由基因與表型的相關(guān)性決定红淡。當(dāng)分值積累到最大就是富集分數(shù)不狮。
X軸是實驗中的所有基因（在這里大約為20,000）。每個黑條是該基因集中的基因（途徑）在旱，我們可以知道基因在排序列表中的位置摇零。
如果基因集位于預(yù)先排列的基因列表的頂部，則通過某種度量計算出富集分數(shù)（enrichment score桶蝎，ES）驻仅，ES為正。如果基因集位于預(yù)先排列的基因列表的底部登渣，則ES為負噪服。
從以上圖中可以看出多數(shù)基因都落在了峰值之前（綠線峰值）的核心基因集中，表明基因在該數(shù)據(jù)集中的顯著富集胜茧，也就是高表達芯咧。

4. 常規(guī)GSEA分析

GSEA（軟件）分析需要準備的文件：

• 基因表達量表（.gct / .txt）文件
• 樣本分類關(guān)系表（.cls）文件
• 基因集合（.gmt）文件

• 4.1 基因表達量表文件

新建excel，將數(shù)據(jù)整理成下面的格式竹揍，存儲為.txt文件敬飒。

第一列是NAME，也就是基因名芬位；第二列是Description无拗，全是gene即可，不需更改昧碉；后面的列是樣本的表達量英染，需要是數(shù)值類型。

• 4.2 樣本分類關(guān)系表（.cls）文件

新建excel
第一行分別是：樣本數(shù)被饿、幾個分組、1(固定不需要改)
第二行分別是：#狭握、分組1闪金、分組2
第三行是每個組的樣本，0代表CON組，1代表CASE組

將文件保存成.txt格式哎垦，隨后更改后綴為.cls （一定要先保存成.txt格式囱嫩，若先存成.xlsx格式再改成.cls，讀入GSEA軟件會報錯）

• 4.3 基因集合（.gmt）文件
  人的文件可以直接從GSEA官網(wǎng)上下載
  小鼠的文件（以KEGG上下載的數(shù)據(jù)集做演示）：
  首先漏设，在KEGG數(shù)據(jù)庫中下載小鼠的KEGG通路數(shù)據(jù)庫文件

點擊mmu

點擊Brite hierachy

這樣就得到了下載.keg文件

先在Linux下解析keg文件

#解析信息
perl -alne '{if(/^C/){/PATH:mmu(\d+)/;$kegg=$1}else{print "$kegg\t$F[1]" if /^D/ and $kegg;}}' mmu00001.keg > mmu00001_kegg2gene.txt
grep '^C' mmu00001.keg | grep "mmu" | sed 's/^C    //g' | sed 's/\([0-9]*\)\s/\1\t/' > mmu00001_kegg2description.txt

再用R制作.gmt文件

#轉(zhuǎn)換ID同時制作.gmt文件
#安裝轉(zhuǎn)換ID需要用到的各個R包
install.packages("data.table")
BiocManager::install("org.Mm.eg.db")
BiocManager::install("clusterProfiler")

#調(diào)用各個R包
library(data.table)
library(org.Mm.eg.db)
library(clusterProfiler)

#讀入數(shù)據(jù)
path2gene_file<-"mmu00001_kegg2gene.txt"
tmp=read.table(path2gene_file,sep="\t",colClasses=c('character'))
dt <- tmp
dt$geneid <- rownames(dt)

# transform id  
map_dt <- bitr(dt$V2, fromType = "ENTREZID",toType = c( "SYMBOL"),OrgDb = org.Mm.eg.db)
dt_merge <- merge(map_dt,dt, by.y = "V2", by.x = "ENTREZID")

# 可以把轉(zhuǎn)換ID后的表輸出
#write.table(dt_merge, file = "example42.txt", sep = "\t", col.names = T, row.names = T, quote = F)

# first column is kegg ID, second column is entrez ID
GeneID2kegg_list<- tapply(tmp[,1],as.factor(tmp[,2]),function(x) x)
kegg2symbol_list<- tapply(dt_merge$SYMBOL,as.factor(dt_merge$V1),function(x) x)
View(kegg2symbol_list)

#輸出數(shù)據(jù)墨闲，但描述列為NA
write.gmt <- function(geneSet=kegg2symbol_list,gmt_file='kegg2symbol.gmt'){
  sink( gmt_file )
  for (i in 1:length(geneSet)){
    cat(names(geneSet)[i])
    cat('\tNA\t')
    cat(paste(geneSet[[i]],collapse = '\t'))
    cat('\n')
  }
  sink()
}

write.gmt()

在上面的gmt文件里添加自定義基因集：

gene50 <- matrix$mygene[1:50]
GeneID2kegg_list$gene50 <- gene50

然后重復(fù)上面的輸出gmt文件。如果只想要自定義基因集郑口，可以在加入自定義基因集后刪掉其他基因集鸳碧。

• 4.4 排序文件
  前面我們得到了GSEA軟件需要用的三個文件，但GSEA軟件的圖默認是.png犬性，且不夠清晰杆兵。如果想使用R畫圖，除了需要上面的.gmt文件還需要排序文件

rm(list = ls())  
options(stringsAsFactors = F)
# library(org.Hs.eg.db)
library(org.Mm.eg.db)
library(clusterProfiler)
library(enrichplot)
library(tidyverse)
library(ggstatsplot)

#讀入差異分析結(jié)果
DEG_DESeq2_all <- read_excel("DEG_DESeq2_all.xlsx")

## 物種設(shè)置
organism = 'hsa'    #  人類'hsa' 小鼠'mmu'   
OrgDb = 'org.Hs.eg.db'#人類"org.Hs.eg.db" 小鼠"org.Mm.eg.db"

#### 按照需要可選擇不同的DEG方法數(shù)據(jù)集 ####
need_DEG <- DEG_DESeq2_all
need_DEG <- need_DEG[,c(2,4,7)] #選擇SYMBOL, log2FoldChange和pvalue（湊成數(shù)據(jù)框）

colnames(need_DEG) <- c('SYMBOL','log2FoldChange','pvalue')

##### 創(chuàng)建gsea分析的geneList（包含從大到小排列的log2FoldChange和ENTREZID信息）####
#轉(zhuǎn)化id  
df <- bitr(need_DEG$SYMBOL, fromType = "SYMBOL",toType =  "ENTREZID", OrgDb = 'org.Hs.eg.db') #人數(shù)據(jù)庫org.Hs.eg.db 小鼠org.Mm.eg.db
need_DEG <- merge(need_DEG, df, by='SYMBOL')  #按照SYMBOL合并注釋信息
need_DEG <- filter(need_DEG,need_DEG$log2FoldChange!='Inf')
need_DEG <- filter(need_DEG,need_DEG$log2FoldChange!='=Inf')
geneList <- need_DEG$log2FoldChange
#geneList <- as.numeric(geneList)
names(geneList) <- need_DEG$ENTREZID
geneList <- sort(geneList, decreasing = T)   #從大到小排序

這樣得到的是以ENTREZID命名的"numeric"仔夺，接下來就可以進行常規(guī)GSEA分析

##### gsea富集 ####
KEGG_kk_entrez <- gseKEGG(geneList     = geneList,
                          organism     = 'mmu', #人hsa 鼠mmu
                          pvalueCutoff = 0.25)  #實際為padj閾值,可調(diào)整 
KEGG_kk <- DOSE::setReadable(KEGG_kk_entrez, 
                             OrgDb=OrgDb,
                             keyType='ENTREZID')#轉(zhuǎn)化id             

GO_kk_entrez <- gseGO(geneList     = geneList,
                      ont          = "ALL",  # "BP"琐脏、"MF"和"CC"或"ALL"
                      OrgDb        = 'org.Mm.eg.db',#人類org.Hs.eg.db 鼠org.Mm.eg.db
                      keyType      = "ENTREZID",
                      pvalueCutoff = 0.25)   #實際為padj閾值可調(diào)整
GO_kk <- DOSE::setReadable(GO_kk_entrez, 
                           OrgDb=OrgDb,
                           keyType='ENTREZID')#轉(zhuǎn)化id 

#save(KEGG_kk_entrez, GO_kk_entrez, file = "GSEA_result.RData")


###可視化
##選取富集結(jié)果
kk_gse <- KEGG_kk
kk_gse_entrez <- KEGG_kk_entrez

###條件篩選 
#一般認為|NES|>1，NOM pvalue<0.05缸兔，F(xiàn)DR（padj）<0.25的通路是顯著富集的
kk_gse_cut <- kk_gse[kk_gse$pvalue<0.05 & kk_gse$p.adjust<0.25 & abs(kk_gse$NES)>1]
kk_gse_cut_down <- kk_gse_cut[kk_gse_cut$NES < 0,]
kk_gse_cut_up <- kk_gse_cut[kk_gse_cut$NES > 0,]

#選擇展現(xiàn)NES前幾個通路 
down_gsea <- kk_gse_cut_down[tail(order(kk_gse_cut_down$NES,decreasing = T),10),]
up_gsea <- kk_gse_cut_up[head(order(kk_gse_cut_up$NES,decreasing = T),10),]
diff_gsea <- kk_gse_cut[head(order(abs(kk_gse_cut$NES),decreasing = T),10),]


#### 經(jīng)典的GSEA圖 
up_gsea$Description
i=2
gseap1 <- gseaplot2(kk_gse,
                    up_gsea$ID[i],#富集的ID編號
                    title = up_gsea$Description[i],#標題
                    color = "red", #GSEA線條顏色
                    base_size = 20,#基礎(chǔ)字體大小
                    rel_heights = c(1.5, 0.5, 1),#副圖的相對高度
                    subplots = 1:3,   #要顯示哪些副圖 如subplots=c(1,3) #只要第一和第三個圖
                    ES_geom = "line", #enrichment score用線還是用點"dot"
                    pvalue_table = T) #顯示pvalue等信息
gseap1
ggsave(gseap1, filename = 'GSEA_up_1.pdf', width =10, height =8)

#### 合并 GSEA通路 
gseap2 <- gseaplot2(kk_gse,
                    up_gsea$ID,#富集的ID編號
                    title = "UP_GSEA_all",#標題
                    color = "red",#GSEA線條顏色
                    base_size = 20,#基礎(chǔ)字體大小
                    rel_heights = c(1.5, 0.5, 1),#副圖的相對高度
                    subplots = 1:3, #要顯示哪些副圖 如subplots=c(1,3) #只要第一和第三個圖
                    ES_geom = "line",#enrichment score用線還是用點"dot"
                    pvalue_table = T) #顯示pvalue等信息
gseap2 
ggsave(gseap2, filename = "GSEA_up_all.pdf",width =12,height =12)

但如果想使用4.3中的自定義gmt文件日裙，則需要注意，如果生成的gmt文件下使用的是基因SYMBOL惰蜜，則上面構(gòu)建排序文件時不需將SYMBOL轉(zhuǎn)換成ENTREZID昂拂。
在差異分析結(jié)果中取出'SYMBOL'和'log2FoldChange'后，將SYMBOL去重復(fù)抛猖，再單獨取出log2FoldChange的值格侯，使用SYMBOL去對log2FoldChange命名，然后按log2FoldChange排序即可财著。

need_DEG=need_DEG[duplicated(need_DEG$SYMBOL)==FALSE,]
geneList <- need_DEG$log2FoldChange
names(geneList) <- a$SYMBOL
geneList <- sort(geneList, decreasing = T)

隨后就可以使用fgsea進行分析了

#pathway是讀入的.gmt文件
fgseaRes <- fgsea(pathway, geneList, minSize=15, maxSize=500, nperm=1000) 
fgseaResTidy <- fgseaRes %>%  as_tibble() %>% arrange(desc(NES))
fgseaResTidy %>% dplyr::select(-leadingEdge, -ES, -nMoreExtreme) %>% arrange(padj) %>% head()
plotEnrichment(pathway[["gene50"]],geneList) + labs(title="gene50")

?該教程中的一些小細節(jié)：

1. wilcoxauc函數(shù)
這個函數(shù)來自presto包联四，可以基于高斯近似法計算Wilcoxon p值和auROC〕沤蹋可以輸入Dense matrix/data.frame朝墩，Sparse matrix比如dgCMatrix，Seurat V3對象和SingleCellExperiment對象伟姐。

2. dplyr：
參考：dplyr包的函數(shù)及用法
這里使用到了dplyr包中的4個函數(shù)：
filter函數(shù)：針對行進行操作收苏，提取一個或多個分組變量中的某個觀測
arrange函數(shù)：針對某個變量對矩陣進行排序，默認升序
select函數(shù)：針對列進行操作愤兵，提取指定的一或多列
count函數(shù)：對數(shù)據(jù)集中某個分組變量的各個觀測值的數(shù)量進行統(tǒng)計

3. fgsea函數(shù)：
fgsea能夠快速對預(yù)選基因集進行GSEA富集分析鹿霸，預(yù)選基因集可以是自己設(shè)定，一般使用MSigdb數(shù)據(jù)庫（同樣由提出GSEA方法團隊提供）秆乳。
fgsea（）函數(shù)需要一個基因集列表以及對應(yīng)值懦鼠，主要是基因名和AUC(ROC曲線下方的面積大小，簡單說就是代表準確性，準確性越高葛闷，AUC值越大)憋槐，其中基因集中的基因名要與數(shù)據(jù)集（pathway）中的基因名相對應(yīng)双藕。
fgsea包中的plotEnrichment函數(shù)用于GSEA圖的繪制淑趾。
fgsea函數(shù)詳細參考：https://stephenturner.github.io/deseq-to-fgsea/

4. msigdbr：
參考：基因功能富集方法和基因注釋數(shù)據(jù)庫介紹