简 介
软件包功能允许用户从存储的数据或其他来源输入单细胞RNA-SEQ 计数和任何基因集路径。富集计算本身使用两种方法:
1)gsva R包和RNA的泊松分布;
2)UCell包。
在单细胞RNA测序的背景下进基因集富集分析(GSEA)。使用 raw count information, Seurat objects, or SingleCellExperiment format 作为输入，用户可以对单细胞数据进行分析并可视化GSEA结果。
我们使用这篇文字的数据集进行测试。
软件包安装
if (!requireNamespace("BiocManager", quietly = TRUE))
install.packages("BiocManager")
options(BioC_mirror="https://mirrors.tuna.tsinghua.edu.cn/bioconductor")
BiocManager::install("dittoSeq")
devtools::install_github("dtm2451/dittoSeq")
BiocManager::install("escape")
devtools::install_github("ncborcherding/escape@dev")
devtools::install_github("ncborcherding/escape")数据读取
准备好的数据集：https://ncborcherding.github.io/vignettes/escape_seurat_ex.rds 直接下载即可。该数据是一个更大的、已发表的皮肤T细胞淋巴瘤项目(CTCL)的一个子集，该项目包含1141个非恶性或正常T细胞(N)和847个恶性T细胞(T)。我们已经将这些数据整合并聚集成8个不同的集群。将更多地了解细胞类型之间的差异。
直接读取.rds文件，并查看单细胞数据分布：
library(escape)
library(Seurat)
library(dittoSeq)
library(ggplot2)
library(grDevices)
seurat_ex <- readRDS("escape_seurat_ex.rds")
seurat_ex <- UpdateSeuratObject(seurat_ex)
DimPlot(seurat_ex, label = T) + NoLegend()
查看样本分类：
colorblind_vector <- colorRampPalette(rev(c("#0D0887FF", "#47039FFF", "#7301A8FF",
"#9C179EFF", "#BD3786FF", "#D8576BFF", "#ED7953FF", "#FA9E3BFF", "#FDC926FF",
"#F0F921FF")))
DimPlot(seurat_ex, group.by = "Type") + scale_color_manual(values = colorblind_vector(2))
设置数据集
Option 1: Molecular Signture Database
进行基因集富集分析的第一步是确定我们想要使用的基因集。函数getGeneSets()允许用户从GSEABase GeneSetCollection 对象列表中选择整个或多个库。可以通过将参数library设置为library/libraries of interest来对内置 Molecular Signature Database 中的基因集集合进行此操作。
此外: -通过设置基因，可以从选定的文库中分离出单个通路/基因集。gene.sets =感兴趣的基因集合的名称。
—可通过subcategory参数选择单个库的子类别。
如果单细胞数据的测序是在“智人”以外的物种上执行的，请确保使用getGeneSets()中的物种参数，以获得正确的基因命名法。
gene.sets <- getGeneSets(library = "H")
gene.sets
## GeneSetCollection
## names: HALLMARK_ADIPOGENESIS, HALLMARK_ALLOGRAFT_REJECTION, ..., HALLMARK_XENOBIOTIC_METABOLISM (50 total)
## unique identifiers: ABCA1, ABCB8, ..., XDH (4383 total)
## types in collection:
## geneIdType: NullIdentifier (1 total)
## collectionType: NullCollection (1 total)Option 2: Built-In gene sets
data("escape.gene.sets", package="escape")
gene.sets <- escape.gene.setsOption 3: Define personal gene sets
gene.sets <- list(Tcell_signature = c("CD2","CD3E","CD3D"),
Myeloid_signature = c("SPI1","FCER1G","CSF1R"))实例操作
下一步是对RNA计数数据进行富集。函数enrichIt()既可以处理原始计数数据的矩阵，也可以直接从singlecellexexperiment或Seurat对象中提取数据。的基因。函数中的set参数是GeneSetCollection，可以从getGeneSets()生成，也可以从用户生成。浓缩分数将在所有单个细胞中计算，组是打破浓缩的n个大小，而核心是在浓缩计算期间并行执行的核心数量。
enrichIt()可以使用两种不同的方法来使用方法参数进行量化- Barbie等人2009年描述的“ssGSEA”方法或Andreatta和Carmona 2021年描述的“UCell”方法。
为了防止对单细胞对象计数数据中表示的基因数量较少的基因集进行富集计算时出现问题，建立了min.size来去除少于指示基因的基因集。
ES <- enrichIt(obj = seurat_ex, gene.sets = gene.sets, groups = 1000, cores = 1,
min.size = NULL)
## [1] "Using sets of 1000 cells. Running 2 times."
## Setting parallel calculations through a SnowParam back-end
## with workers=14 and tasks=100.
## Estimating ssGSEA scores for 50 gene sets.
## [1] "Calculating ranks..."
## [1] "Calculating absolute values from ranks..."
## Setting parallel calculations through a SnowParam back-end
## with workers=14 and tasks=100.
## Estimating ssGSEA scores for 50 gene sets.
## [1] "Calculating ranks..."
## [1] "Calculating absolute values from ranks..."
seurat_ex <- AddMetaData(seurat_ex, ES)结果可视化
绘制热图
所有富集结果
seurat_ex@meta.data$active.idents <- seurat_ex@active.ident
table(seurat_ex@meta.data$Type)
##
## N T
## 1141 847
dittoHeatmap(seurat_ex, genes = NULL, metas = names(ES), heatmap.colors = rev(colorblind_vector(50)),
annot.by = c("active.idents", "Type"), cluster_cols = TRUE, fontsize = 7)
筛选特点基因集
dittoHeatmap(seurat_ex, genes = NULL, metas = c("HALLMARK_APOPTOSIS", "HALLMARK_DNA_REPAIR",
"HALLMARK_P53_PATHWAY"), heatmap.colors = rev(colorblind_vector(50)), annot.by = c("active.idents",
"Type"), cluster_cols = TRUE, fontsize = 7)
绘制小提琴图
dittoPlot(seurat_ex, "HALLMARK_DNA_REPAIR", group.by = "Type") + scale_fill_manual(values = colorblind_vector(2))
绘制密度富集图
dittoScatterHex(seurat_ex, x.var = "HALLMARK_DNA_REPAIR", y.var = "HALLMARK_MTORC1_SIGNALING",
do.contour = TRUE) + theme_classic() + scale_fill_gradientn(colors = rev(colorblind_vector(11))) +
geom_vline(xintercept = 0, lty = 2) + geom_hline(yintercept = 0, lty = 2)
分组绘制密度富集图
dittoScatterHex(seurat_ex, x.var = "HALLMARK_DNA_REPAIR", y.var = "HALLMARK_MTORC1_SIGNALING",
do.contour = TRUE, split.by = "Type") + theme_classic() + scale_fill_gradientn(colors = rev(colorblind_vector(11))) +
geom_vline(xintercept = 0, lty = 2) + geom_hline(yintercept = 0, lty = 2)
绘制山脊富集图
ES2 <- data.frame(seurat_ex[[]], Idents(seurat_ex))
colnames(ES2)[ncol(ES2)] <- "cluster"
ridgeEnrichment(ES2, gene.set = "HALLMARK_DNA_REPAIR", group = "Type", add.rug = TRUE)
绘制细胞类型山脊富集图
ridgeEnrichment(ES2, gene.set = "HALLMARK_DNA_REPAIR", group = "cluster", facet = "Type",
add.rug = TRUE)
绘制分组小提琴图
与上面一样，我们可以通过使用ES2对象调用SplitEnrichment()来探索“HALLMARK_DNA_REPAIR”在类型之间的分布。我们指定split = " Type " "，这将分离小提琴图本身的肿瘤细胞和正常细胞。
splitEnrichment(ES2, split = "Type", gene.set = "HALLMARK_DNA_REPAIR")
也可以通过指定x轴=“cluster”来探索集群分布。
splitEnrichment(ES2, x.axis = "cluster", split = "Type", gene.set = "HALLMARK_DNA_REPAIR")## Warning: Groups with fewer than two data points have been dropped.
## Groups with fewer than two data points have been dropped.
## Groups with fewer than two data points have been dropped.
绘制Enrichment Plots
enrichmentPlot(seurat_ex, gene.set = "HALLMARK_DNA_REPAIR", gene.sets = gene.sets,
group = "Type") + scale_color_manual(values = colorblind_vector(5)[c(1, 4)])
分组免疫浸润差异
output <- getSignificance(ES2, group = "Type", fit = "T.test")
T.statistic p.value FDR median.N median.T
nCount_RNA -16.790749 2.187225e-58 7.655289e-57 5073.0000 7330.000
nFeature_RNA -24.099722 2.027481e-110 9.529161e-109 1487.0000 2163.000
HALLMARK_ADIPOGENESIS -11.233782 2.341109e-28 3.979886e-27 1059.5467 1293.156
HALLMARK_ALLOGRAFT_REJECTION -12.024954 4.444757e-32 8.000562e-31 2235.5484 2462.286
HALLMARK_ANDROGEN_RESPONSE -14.716000 3.056918e-46 7.642294e-45 1715.6572 1963.711
HALLMARK_ANGIOGENESIS -9.952924 1.000529e-22 1.400741e-21 234.0722 544.154Reference
Borcherding, N., Vishwakarma, A., Voigt, A.P. et al. Mapping the immune environment in clear cell renal carcinoma by single-cell genomics. Commun Biol 4, 122 (2021). https://doi.org/10.1038/s42003-020-01625-6.
单细胞生信分析教程
桓峰基因公众号推出单细胞生信分析教程并配有视频在线教程，目前整理出来的相关教程目录如下：
SCS【4】单细胞转录组数据可视化分析 (Seurat 4.0)
SCS【6】单细胞转录组之细胞类型自动注释 (SingleR)
SCS【7】单细胞转录组之轨迹分析 (Monocle 3) 聚类、分类和计数细胞
SCS【8】单细胞转录组之筛选标记基因 (Monocle 3)
SCS【9】单细胞转录组之构建细胞轨迹 (Monocle 3)
SCS【10】单细胞转录组之差异表达分析 (Monocle 3)
SCS【11】单细胞ATAC-seq 可视化分析 (Cicero)
SCS【12】单细胞转录组之评估不同单细胞亚群的分化潜能 (Cytotrace)
SCS【13】单细胞转录组之识别细胞对“基因集”的响应 (AUCell)
SCS【15】细胞交互：受体-配体及其相互作用的细胞通讯数据库 (CellPhoneDB)
SCS【16】从肿瘤单细胞RNA-Seq数据中推断拷贝数变化 (inferCNV)
SCS【17】从单细胞转录组推断肿瘤的CNV和亚克隆 (copyKAT)
SCS【18】细胞交互：受体-配体及其相互作用的细胞通讯数据库 (iTALK)
SCS【21】单细胞空间转录组可视化 (Seurat V5)
SCS【22】单细胞转录组之 RNA 速度估计 (Velocyto.R)
SCS【24】单细胞数据量化代谢的计算方法 (scMetabolism)
SCS【25】单细胞细胞间通信第一部分细胞通讯可视化（CellChat）
SCS【26】单细胞细胞间通信第二部分通信网络的系统分析（CellChat）
SCS【28】单细胞转录组加权基因共表达网络分析(hdWGCNA)
SCS【29】单细胞基因富集分析 (singleseqgset)
SCS【30】单细胞空间转录组学数据库(STOmics DB)
SCS【31】减少障碍，加速单细胞研究数据库(Single Cell PORTAL)
SCS【32】基于scRNA-seq数据中推断单细胞的eQTLs (eQTLsingle)
SCS【33】单细胞转录之全自动超快速的细胞类型鉴定 (ScType)
SCS【34】单细胞/T细胞/抗体免疫库数据分析(immunarch)
SCS【35】单细胞转录组之去除双细胞 (DoubletFinder)
SCS【36】单细胞转录组之k-近邻图差异丰度测试(miloR)
利用这个软件包实现了快速计算单细胞免疫浸润值，有需求的老师可以联系桓峰基因，关注桓峰基因公众号，轻松学生信，高效发文章！
号外号外，桓峰基因单细胞生信分析免费培训课程即将开始，快来报名吧！
桓峰基因，铸造成功的您！
未来桓峰基因公众号将不间断的推出单细胞系列生信分析教程，
敬请期待！！
桓峰基因官网正式上线，请大家多多关注，还有很多不足之处，大家多多指正！
http://www.kyohogene.com/
桓峰基因和投必得合作，文章润色优惠85折，需要文章润色的老师可以直接到网站输入领取桓峰基因专属优惠券码：KYOHOGENE，然后上传，付款时选择桓峰基因优惠券即可享受85折优惠哦！https://www.topeditsci.com/
794
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
