SCS【29】单细胞基因富集分析 (singleseqgset)

最新推荐文章于 2024-05-31 09:30:00 发布

桓峰基因

最新推荐文章于 2024-05-31 09:30:00 发布

阅读量1.1k

点赞数 1

文章标签：信息可视化

本文链接：https://blog.csdn.net/weixin_41368414/article/details/131388874

版权

桓峰基因公众号推出单细胞生信分析教程并配有视频在线教程，目前整理出来的相关教程目录如下：

Topic 6. 克隆进化之 Canopy

Topic 7. 克隆进化之 Cardelino

Topic 8. 克隆进化之 RobustClone

SCS【1】今天开启单细胞之旅，述说单细胞测序的前世今生

SCS【2】单细胞转录组之 cellranger

SCS【3】单细胞转录组数据 GEO 下载及读取

SCS【4】单细胞转录组数据可视化分析 (Seurat 4.0)

SCS【5】单细胞转录组数据可视化分析 (scater)

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【14】单细胞调节网络推理和聚类 (SCENIC)

SCS【15】细胞交互：受体-配体及其相互作用的细胞通讯数据库 (CellPhoneDB)

SCS【16】从肿瘤单细胞RNA-Seq数据中推断拷贝数变化 (inferCNV)

SCS【17】从单细胞转录组推断肿瘤的CNV和亚克隆 (copyKAT)

SCS【18】细胞交互：受体-配体及其相互作用的细胞通讯数据库 (iTALK)

SCS【19】单细胞自动注释细胞类型 (Symphony)

SCS【20】单细胞数据估计组织中细胞类型(Music)

SCS【21】单细胞空间转录组可视化 (Seurat V5)

SCS【22】单细胞转录组之 RNA 速度估计 (Velocyto.R)

SCS【23】单细胞转录组之数据整合 (Harmony)

SCS【24】单细胞数据量化代谢的计算方法 (scMetabolism)

SCS【25】单细胞细胞间通信第一部分细胞通讯可视化（CellChat）

SCS【26】单细胞细胞间通信第二部分通信网络的系统分析（CellChat）

SCS【27】单细胞转录组之识别标记基因（scran）

SCS【28】单细胞转录组加权基因共表达网络分析(hdWGCNA)

简介

基因集富集分析是一种普遍的工具，用于从基因组数据集中提取有生物学意义的信息。有许多不同风格的工具可用于基因集富集分析，但最常见的是 Subramanian 等人的开创性工作发表在 PNAS 2005。他们的关键概念在任何基因集富集分析是比较一个指标(如对数倍的变化在组之间的基因表达)的基因集内与基因集外的基因。基因集可以从许多地方获得，包括 Broad 的分子特征数据库 (MSigDB) 和 The Gene Ontology Resource。选择最能回答手头生物学问题的基因组是很重要的。例如，你通常不会使用来自Broad的C7(免疫基因集)的基因集来回答有关神经元发育的问题(除非你有一个非常好的理由!)。我在这里开发的基因集富集方法来自于 Di Wu 和 Gordon K Smyth 的成果发表在 Nucleic Acids Research 2012。Wu和Smyth开发了一种相关调整平均秩基因集测试(CAMERA)。CAMERA是一种竞争性基因集富集测试，控制基因集内基因间的相关性。之所以选择使用CAMERA，是因为它不依赖于基因独立性的假设(因为我们知道基因通常具有结构化的共表达模式)，也不依赖于基因标签的排列或测试样本的重新采样。CAMERA 通过产生基于基因间相关程度的方差膨胀因子来控制基因间相关性，并将这种方差膨胀纳入 Wilcoxon 秩和检验。

演示了 singleseqgset 在模拟数据上的标准用法，主要有以下步骤:

1).使用R软件包 splatter 模拟表达数据;

2).使用 msigdbr下载感兴趣的基因集;

3).将特定的基因集添加到模拟数据中;

4).使用标准的Seurat工作流(v.2.3.4)处理数据;

5).使用singleseqgset执行基因集富集分析;

6).将结果绘制在热图中。

注：我们改为Seurat使用V4版，在利用Seurat时修改了里面涉及的函数及代码。另外这个软件包里面有一个bug，在我们桓峰基因公众号已经给出解决方法，有报错的可以参考:Debug 1. 单细胞富集分析 (singleseqgset) 报错怎么解？

软件包安装

这里我们需要安装几个软件包，主要包括 singleseqgset， splatter，heatmap3，msigdbr，若没有的话，自己安装即可。

library(devtools)
install_github("arc85/singleseqgset")
BiocManager::install("splatter")

例子实操

1).使用R软件包 splatter 模拟表达数据;

首先，加载所有必要的包，并模拟与 splatter 使用的数据。创建了一些由4个簇组成的模拟数据，簇之间有不同比例的基因差异表达。

suppressMessages({
    library(splatter)
    library(Seurat)
    options(Seurat.object.assay.version = "v4")
    library(msigdbr)
    library(singleseqgset)
    library(heatmap3)
})

# Create parameters and simulate data
sc.params <- newSplatParams(nGenes = 1000, batchCells = 5000)
sim.groups <- splatSimulate(params = sc.params, method = "groups", group.prob = c(0.4,
    0.2, 0.3, 0.1), de.prob = c(0.2, 0.2, 0.1, 0.3), verbose = F)
sim.groups  #Check out the SingleCellExperiment object with simulated dataset
## class: SingleCellExperiment 
## dim: 1000 5000 
## metadata(1): Params
## assays(6): BatchCellMeans BaseCellMeans ... TrueCounts counts
## rownames(1000): Gene1 Gene2 ... Gene999 Gene1000
## rowData names(8): Gene BaseGeneMean ... DEFacGroup3 DEFacGroup4
## colnames(5000): Cell1 Cell2 ... Cell4999 Cell5000
## colData names(4): Cell Batch Group ExpLibSize
## reducedDimNames(0):
## mainExpName: NULL
## altExpNames(0):

colData(sim.groups)  #The colData column 'Group' contains the groups of simulated cells
## DataFrame with 5000 rows and 4 columns
##                 Cell       Batch    Group ExpLibSize
##          <character> <character> <factor>  <numeric>
## Cell1          Cell1      Batch1   Group3    62905.5
## Cell2          Cell2      Batch1   Group4    48401.0
## Cell3          Cell3      Batch1   Group1    64592.9
## Cell4          Cell4      Batch1   Group1    56475.9
## Cell5          Cell5      Batch1   Group3    82686.9
## ...              ...         ...      ...        ...
## Cell4996    Cell4996      Batch1   Group2    59076.6
## Cell4997    Cell4997      Batch1   Group4    84844.4
## Cell4998    Cell4998      Batch1   Group1    63239.6
## Cell4999    Cell4999      Batch1   Group3    75500.1
## Cell5000    Cell5000      Batch1   Group3    55475.8
# We will pull the simulated counts and groups from the sim.groups object
sim.counts <- assays(sim.groups)$counts
dim(sim.counts)
## [1] 1000 5000
groups <- colData(sim.groups)$Group
names(groups) <- rownames(colData(sim.groups))
table(groups)
## groups
## Group1 Group2 Group3 Group4 
##   1960   1008   1494    538

2).使用 msigdbr下载感兴趣的基因集;

我们将使用misigdbr软件包下载人类贺曼基因集。或者，我们可以直接从Broad的网站:http://software.broadinstitute.org/gsea/msigdb/index.jsp下载基因集，并使用GSEABase软件包将它们读取到R中。

注意:必须选择合适的基因注释样式(即基因符号或entrez基因id;这里我们只是基因符号)和物种(这里我们使用人类)为您的项目。否则，您的基因集中的基因名称将与您的数据中的基因名称不匹配，并且所有基因集将获得“0”富集分数。

library(msigdbr)
h.human <- msigdbr(species = "Homo sapiens", category = "H")

h.names <- unique(h.human$gs_name)

h.sets <- vector("list", length = length(h.names))
names(h.sets) <- h.names

library(dplyr)
for (i in names(h.sets)) {
    h.sets[[i]] <- pull(h.human[h.human$gs_name == i, "gene_symbol"])
}

3).将特定的基因集添加到模拟数据中;

我们将指定5个基因集在每个集群中独占表达。我们将随机选择这20组基因，每个基因将从零膨胀的泊松分布中抽取，成功概率等于50%，lambda值为10。这将模拟中度辍学和相对较低的基因表达计数。我们已经成功地修改了模拟计数矩阵，以包含特定集群中感兴趣的基因集。

# Add gene sets to simulated clusters
sets.to.use <- sample(seq(1, 50, 1), 20, replace = F)
sets.and.groups <- data.frame(set = sets.to.use, group = paste("Group", rep(1:4,
    each = 5), sep = ""))

for (i in 1:nrow(sets.and.groups)) {

    set.to.use <- sets.and.groups[i, "set"]
    group.to.use <- sets.and.groups[i, "group"]
    gene.set.length <- length(h.sets[[set.to.use]])
    blank.matrix <- matrix(0, ncol = ncol(sim.counts), nrow = gene.set.length)
    rownames(blank.matrix) <- h.sets[[sets.to.use[i]]]
    matching <- rownames(blank.matrix) %in% rownames(sim.counts)

    if (any(matching == TRUE)) {

        matching.genes <- rownames(blank.matrix)[matching]
        nonmatching.genes <- setdiff(rownames(blank.matrix), matching.genes)
        blank.matrix <- blank.matrix[!matching, ]
        sim.counts <- rbind(sim.counts, blank.matrix)

    } else {

        sim.counts <- rbind(sim.counts, blank.matrix)
        matching.genes <- rownames(blank.matrix)

    }

    group.cells <- colnames(sim.counts)[groups == group.to.use]

    for (z in group.cells) {

        if (any(matching == TRUE)) {

            sim.counts[matching.genes, z] <- ifelse(rbinom(length(matching.genes),
                size = 1, prob = 0.5) > 0, 0, rpois(length(matching.genes), lambda = 10))
            sim.counts[nonmatching.genes, z] <- ifelse(rbinom(length(nonmatching.genes),
                size = 1, prob = 0.5) > 0, 0, rpois(length(nonmatching.genes), lambda = 10))

        } else {

            sim.counts[matching.genes, z] <- ifelse(rbinom(length(matching.genes),
                size = 1, prob = 0.5) > 0, 0, rpois(length(matching.genes), lambda = 10))

        }
    }

}
dim(sim.counts)
## [1] 3547 5000
# Here, we will check out the sum of expression for the first gene set
len.set1 <- length(h.sets[[sets.to.use[[1]]]])
plot(apply(sim.counts[1001:(1000 + len.set1), ], 2, sum), xlim = c(0, 5000), xlab = "Cells",
    ylab = "Sum of Gene Set 1 Expression")

4).使用标准的Seurat工作流(v.2.3.4)处理数据;

源代码使用的是Seurat工作流(v.2.3.4)处理数据，而我们为了方便将使用标准的Seurat工作流(v4)来可视化我们的模拟数据。

# Pass simulated data to Seurat
sim.counts = sim.counts[!duplicated(rownames(sim.counts)), ]
ser <- CreateSeuratObject(counts = sim.counts)

# Normal Seurat workflow
ser <- NormalizeData(ser)
ser = FindVariableFeatures(ser)
ser@assays$RNA@meta.data$var.features <- rownames(ser@assays$RNA$counts)
ser <- ScaleData(ser, features = VariableFeatures(ser))
# We will pretend like the 1000 simulated genes are the 'variable genes', and
# we will skip FindVariableGenes from Seurat
ser <- RunPCA(ser, features = VariableFeatures(ser))

ElbowPlot(ser)

看起来前4组捕获了大部分的差异

# Looks like the first 4 PCs capture most of the variance We will include the
# first 5

ser <- RunTSNE(ser, dims.use = 1:5)

# Add the simulated data IDs to the Seurat object

data.to.add <- colData(sim.groups)$Group
names(data.to.add) <- colnames(ser@assays$RNA$counts)
ser <- AddMetaData(ser, metadata = data.to.add, col.name = "real_id")
ser <- SetIdent(ser, value = "real_id")

# We will skip clustering with Seurat because we already know the true clusters
# IDs

TSNEPlot(ser, group.by = "real_id")

DimPlot(ser, reduction = "tsne", group.by = "real_id")

5).使用singleseqgset执行基因集富集分析;

首先，我们将计算基因集富集测试的度量，在这种情况下，所有基因簇之间的对数倍变化。我们选择使用已根据库大小进行了校正的日志规范化数据。

# Use singleseqgset to perform gene set enrichment analysis
source("cluster_logfc.R")
logfc.data <- logFC(cluster.ids = as.vector(ser@meta.data$real_id), expr.mat = ser@assays$RNA$data)
names(logfc.data)
## [1] "cluster.cells"  "log.fc.cluster"
source("wmw_gsea.R")
gse.res <- wmw_gsea(expr.mat = ser@assays$RNA$data, cluster.cells = logfc.data[[1]],
    log.fc.cluster = logfc.data[[2]], gene.sets = h.sets)
## [1] "Removed 1 rows with all z-scores equal to zero."
names(gse.res)
## [1] "GSEA_statistics" "GSEA_p_values"
res.stats <- gse.res[["GSEA_statistics"]]
res.pvals <- gse.res[["GSEA_p_values"]]

res.pvals <- apply(res.pvals, 2, p.adjust, method = "fdr")  #Correct for multiple comparisons

res.stats[order(res.stats[, 1], decreasing = TRUE)[1:10], ]  #Top gene sets enriched by z scores
##                                       Group1     Group2     Group3     Group4
## HALLMARK_APICAL_JUNCTION           10.029068 -8.8738897 -7.7777472 -14.888549
## HALLMARK_ADIPOGENESIS               9.649492 -6.9360684 -6.7605911 -10.268786
## HALLMARK_P53_PATHWAY                9.371631 -6.1210394 -7.8104146  -6.825038
## HALLMARK_COAGULATION                8.891660 -6.0562046 -7.0869384  -6.410326
## HALLMARK_TNFA_SIGNALING_VIA_NFKB    8.812246 -6.0724976 -6.0367430  -7.170197
## HALLMARK_COMPLEMENT                 4.210233 -4.9892184 -2.3790894  -6.511253
## HALLMARK_INFLAMMATORY_RESPONSE      1.482738 -3.8412838  0.1291806  -3.670790
## HALLMARK_KRAS_SIGNALING_UP          1.439934 -2.1742014 -2.0405954  -2.236894
## HALLMARK_OXIDATIVE_PHOSPHORYLATION  1.400599 -3.4582284 -0.2656490  -4.504055
## HALLMARK_TGF_BETA_SIGNALING         1.378856 -0.7791616 -2.4287152  -2.424578

res.pvals[order(res.stats[, 1], decreasing = TRUE)[1:10], ]  #Top gene sets by p values
##                                          Group1       Group2       Group3
## HALLMARK_APICAL_JUNCTION           5.565623e-22 6.922181e-18 5.167918e-14
## HALLMARK_ADIPOGENESIS              1.210294e-20 2.822138e-11 6.734063e-11
## HALLMARK_P53_PATHWAY               1.166531e-19 5.061531e-09 4.655014e-14
## HALLMARK_COAGULATION               7.374795e-18 6.208333e-09 8.398073e-12
## HALLMARK_TNFA_SIGNALING_VIA_NFKB   1.202099e-17 6.170863e-09 7.005034e-09
## HALLMARK_COMPLEMENT                8.333511e-05 2.475484e-06 4.475883e-02
## HALLMARK_INFLAMMATORY_RESPONSE     3.076844e-01 3.998155e-04 9.769672e-01
## HALLMARK_KRAS_SIGNALING_UP         3.165010e-01 6.140669e-02 9.196646e-02
## HALLMARK_OXIDATIVE_PHOSPHORYLATION 3.165010e-01 1.665204e-03 9.008132e-01
## HALLMARK_TGF_BETA_SIGNALING        3.165010e-01 5.476498e-01 4.124828e-02
##                                          Group4
## HALLMARK_APICAL_JUNCTION           1.916723e-48
## HALLMARK_ADIPOGENESIS              2.386853e-23
## HALLMARK_P53_PATHWAY               4.785764e-11
## HALLMARK_COAGULATION               5.929350e-10
## HALLMARK_TNFA_SIGNALING_VIA_NFKB   4.586997e-12
## HALLMARK_COMPLEMENT                3.319812e-10
## HALLMARK_INFLAMMATORY_RESPONSE     6.235955e-04
## HALLMARK_KRAS_SIGNALING_UP         3.997969e-02
## HALLMARK_OXIDATIVE_PHOSPHORYLATION 2.041737e-05
## HALLMARK_TGF_BETA_SIGNALING        2.503279e-02

names(h.sets)[sets.to.use[1:5]]  #Compare to the simulate sets we created
## [1] "HALLMARK_P53_PATHWAY"             "HALLMARK_TNFA_SIGNALING_VIA_NFKB"
## [3] "HALLMARK_ADIPOGENESIS"            "HALLMARK_APICAL_JUNCTION"        
## [5] "HALLMARK_COAGULATION"

6).将结果绘制在热图中

恭喜你，你已经用singleseqgset完成了你的第一个基因集富集分析!

# Plot the z scores with heatmap3
heatmap3(res.stats, Colv = NA, cexRow = 0.5, cexCol = 1, scale = "row")

Reference

Subramanian A, Tamayo P, Mootha VK, et al. Gene set enrichment analysis: a knowledge-based approach for interpreting genome-wide expression profiles. Proc Natl Acad Sci U S A. 2005;102(43):15545-15550. doi:10.1073/pnas.0506580102
Wu D, Smyth GK. Camera: a competitive gene set test accounting for inter-gene correlation. Nucleic Acids Res. 2012;40(17):e133. doi:10.1093/nar/gks461

桓峰基因，铸造成功的您！

未来桓峰基因公众号将不间断的推出单细胞系列生信分析教程，

敬请期待！！

桓峰基因官网正式上线，请大家多多关注，还有很多不足之处，大家多多指正！

http://www.kyohogene.com/

桓峰基因和投必得合作，文章润色优惠85折，需要文章润色的老师可以直接到网站输入领取桓峰基因专属优惠券码：KYOHOGENE，然后上传，付款时选择桓峰基因优惠券即可享受85折优惠哦！https://www.topeditsci.com/