单细胞之富集分析-1：单细胞GSEA分析流程
单细胞GSEA分析需要的文件有两个：
1. 单细胞基因表达变化数据（两列，一列是geneid/gene symbol，一列是logFC）（文件格式：.rnk）
2. 目标基因集（文件格式：.gmx或.gmt）一行是一个term/基因集，第一列是基因集的名称，后面是基因分类等，再后面是geneid/gene symbol（要和单细胞矩阵的对应，否则不能识别）
1. 单细胞基因表达变化数据的准备
- 1.1 下载单细胞数据（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标准流程对数据进行预处理得到分群
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)
#质控
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中有多少基因（与矩阵的基因数一致）
# 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分析需要的数据格式（也可输入注释好的细胞群的差异基因）
# 仅选择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. 选择自己数据的物种以及要做的GSEA的数据库类型，准备目标基因集的输入文件
为了进行基因集富集分析，首先需要注释基因集。一个重要的来源是Broad Institute开发的MsigDB。
网址：https://www.gsea-msigdb.org/gsea/msigdb/
library(msigdbr)
library(fgsea)
library(dplyr)
library(ggplot2)
##查看物种的数据 msigdbr_show_species()；#我们使用C7免疫基因集
m_df<- msigdbr(species = "Homo sapiens", category = "C7")
head(m_df)
##将m_df的基因与通路取出并改成一个通路对应相应基因的格式
fgsea_sets<- m_df %>% split(x = .$gene_symbol, f = .$gs_name)
#以gs_name为factor对gene_symbol进行分类，统计落在每个gs_name中的gene_symbol的个数，并生成list。
summary(fgsea_sets)
3. 使用fgsea进行基因集富集并绘图
- 3.1 富集分析
fgseaRes<- fgsea(fgsea_sets, stats = ranks, nperm = 1000)
#nperm设置的是permutation次数
- 3.2 整理数据：
fgseaResTidy <- fgseaRes %>% as_tibble() %>% arrange(desc(NES))
fgseaResTidy %>% dplyr::select(-leadingEdge, -ES, -nMoreExtreme) %>% arrange(padj) %>% head()
- 3.3 绘图
应用标准化富集分数绘制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值的过程。遍历基因集中的基因时，当基因出现在目标基因（横轴的基因）中则加分，反之减分。加减分值由基因与表型的相关性决定。当分值积累到最大就是富集分数。
X轴是实验中的所有基因（在这里大约为20,000）。每个黑条是该基因集中的基因（途径），我们可以知道基因在排序列表中的位置。
如果基因集位于预先排列的基因列表的顶部，则通过某种度量计算出富集分数（enrichment score，ES），ES为正。如果基因集位于预先排列的基因列表的底部，则ES为负。
从以上图中可以看出多数基因都落在了峰值之前（绿线峰值）的核心基因集中，表明基因在该数据集中的显著富集，也就是高表达。
✨该教程中的一些小细节：
1. wilcoxauc函数
这个函数来自presto包，可以基于高斯近似法计算Wilcoxon p值和auROC。可以输入Dense matrix/data.frame，Sparse matrix比如dgCMatrix，Seurat V3对象和SingleCellExperiment对象。
2. dplyr：
参考：dplyr包的函数及用法
这里使用到了dplyr包中的4个函数：
filter函数：针对行进行操作，提取一个或多个分组变量中的某个观测
arrange函数：针对某个变量对矩阵进行排序，默认升序
select函数：针对列进行操作，提取指定的一或多列
count函数：对数据集中某个分组变量的各个观测值的数量进行统计
3. fgsea函数：
fgsea能够快速对预选基因集进行GSEA富集分析，预选基因集可以是自己设定，一般使用MSigdb数据库（同样由提出GSEA方法团队提供）。
fgsea（）函数需要一个基因集列表以及对应值，主要是基因名和AUC(ROC曲线下方的面积大小，简单说就是代表准确性，准确性越高，AUC值越大)，其中基因集中的基因名要与数据集（pathway）中的基因名相对应。
fgsea包中的plotEnrichment函数用于GSEA图的绘制。
fgsea函数详细参考：https://stephenturner.github.io/deseq-to-fgsea/
4. msigdbr：
参考：基因功能富集方法和基因注释数据库介绍
