一个R包完成单细胞基因集富集分析 (全代码)

singleseqgset是用于单细胞RNA-seq数据的基因集富集分析的软件包。它使用简单的基础统计量(variance inflated Wilcoxon秩和检验)来确定不同cluster中感兴趣的基因集的富集。

Installation

library(devtools)
install_github("arc85/singleseqgset")
library(singleseqgset)

Introduction

基因集富集分析是一种无处不在的工具,用于从基因组数据集中提取具有生物学意义的信息。有许多不同类型的工具可用于基因集富集分析,以前的基因组富集分析工具的关键概念是比较组内基因与组外基因的差异(例如,各组之间基因表达的log fold change变化)。但最值得肯定的是Subramanian等人(PNAS 2005)(https://www.ncbi.nlm.nih.gov/pubmed/16199517)的开创性工作-Gene Set Enrichment Analysis(GSEA富集分析 - 界面操作)。

GSEA分析中已知功能的基因集可以从许多地方进行获取,包括Broad的Molecular Signatures Database(MSigDB)The Gene Ontology Resource,关键是选择最能回答当前生物学问题的基因集。举个栗子,通常您不会使用Broad的C7(免疫学基因集)中的基因集来回答有关神经元发育的问题(除非有充分的理由!)。

singleseqgset研究者在这里开发的基因集富集方法灵感来自Di WuGordon K SmythNucleic Acids Research 2012(https://www.ncbi.nlm.nih.gov/pubmed/16199517)上的工作,而WuSmyth在这个工作的基础上又开发了相关调整的MEan RAnk基因集测试(CAMERA)

CAMERA是一项竞争性基因集富集测试,可控制基因集内的基因间相关性。研究团队之所以选择使用CAMERA,是因为它不依赖于基因独立性的假设(因为我们知道基因通常具有结构化的共表达模式),也不依赖于基因标记的排列或测试样品的重采样。CAMERA通过基于基因间相关程度生成方差膨胀因子来控制基因间相关性,并将这种方差膨胀纳入Wilcoxon秩和检验中

此workflow演示了对模拟数据使用singleseqgset的标准用法。我们将执行以下步骤:

  1. 使用splatter R包模拟表达式数据;

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

  3. 将特定基因集添加到我们的模拟数据中(就是让我们的模拟数据可以富含某基因集);

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

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

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

Simulate data using splatter

首先,加载所有必要的程序包,并使用splatter模拟数据。

suppressMessages({
library(splatter)
library(Seurat)
library(msigdbr)
library(singleseqgset)
library(heatmap3)
})
#Splatter是用于模拟单细胞RNA测序count数据的软件包。
#创建参数并模拟数据
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.20,0.20,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):
## spikeNames(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
##          <factor> <character> <factor>        <numeric>
## Cell1       Cell1      Batch1   Group1 73324.2901799195
## Cell2       Cell2      Batch1   Group1 71115.1438815874
## Cell3       Cell3      Batch1   Group3  89689.371004738
## Cell4       Cell4      Batch1   Group3  44896.460028516
## Cell5       Cell5      Batch1   Group3 53956.0668720417
## ...           ...         ...      ...              ...
## Cell4996 Cell4996      Batch1   Group2 79041.7558615305
## Cell4997 Cell4997      Batch1   Group4 66684.3797711752
## Cell4998 Cell4998      Batch1   Group3 70407.9613848431
## Cell4999 Cell4999      Batch1   Group3 64645.3257427214
## Cell5000 Cell5000      Batch1   Group3 62266.9119469426
#我们将从sim.groups对象中提取模拟的计数和组别
sim.counts <- assays(sim.groups)$counts
groups <- colData(sim.groups)$Group
names(groups) <- rownames(colData(sim.groups))
table(groups)
## groups
## Group1 Group2 Group3 Group4
##   1977    992   1536    495

有次我们创建了一些模拟数据,这些数据由4个cluster组成,cluster中不同基因的表达比例不同。

Download gene sets of interest using msigdbr

在这里,我们将使用misigdbr软件包下载Hallmark Gene Sets。另外,我们可以直接从Broad网站(http://software.broadinstitute.org/gsea/msigdb/index.jsp)下载基因集,然后使用`GSEABase`软件包将其读入R。

注意:必须为您的项目选择适当的基因注释样式(gene symbols或entrez gene ID;在这里我们只是gene symbols)和物种(在这里我们使用human)。否则,您的基因组中的基因名称将与您的数据中的基因名称不匹配,并且所有基因组都将获得“0”富集得分。

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.namesfor (i in names(h.sets)) {h.sets[[i]] <- pull(h.human[h.human$gs_name==i,"gene_symbol"])
}

Add gene sets to simulated clusters

我们将分配5个基因集在每个cluster中专门表达(如果不进行分配,则可能出现无功能集在cluster中富集的情况。)。随机选择这20个集合,并从zero-inflatedPoisson分布中抽取每个基因,成功概率等于50%,λ值为10,这将模拟中等dropout和相对较低的基因表达量。

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))}}}#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")

图片

我们已经成功地修改了模拟计数矩阵,使其包含特定cluster中感兴趣的基因集。例如,我们知道Group1细胞应显示出以下基因的富集:

HALLMARK_APICAL_SURFACEHALLMARK_COMPLEMENTHALLMARK_DNA_REPAIRHALLMARK_IL2_STAT5_SIGNALINGHALLMARK_UNFOLDED_PROTEIN_RESPONSE

Seurat workflow on simulated data

#构建seurat object
ser <- CreateSeuratObject(raw.data=sim.counts)#归一化
ser <- NormalizeData(ser)
ser@var.genes <- rownames(ser@raw.data)[1:1000]
ser <- ScaleData(ser,genes.use=ser@var.genes)
## Scaling data matrix
#我们会假设1000个模拟基因是“可变基因”,
#我们将跳过Seurat的FindVariableGenes
ser <- RunPCA(ser,pc.genes=ser@var.genes,do.print=FALSE)PCElbowPlot(ser)

图片

#其中top4PCs代表了多数的差异
#我们将会选择top5 PCs;ser <- RunTSNE(ser,dims.use=1:5)#将模拟的dataID号加入Seurat object
data.to.add <- colData(sim.groups)$Group
names(data.to.add) <- ser@cell.names
ser <- AddMetaData(ser,metadata=data.to.add,col.name="real_id")
ser <- SetAllIdent(ser,id="real_id")#我们将会跳过使用Seurat聚类的步骤,因为我们已知道真正的聚类IDsTSNEPlot(ser,group.by="real_id")

图片

Use singleseqgset to perform gene set enrichment analysis

首先,我们将基于logFC计算基因集富集测试的矩阵。我们选择使用已针对library大小进行校正的log-normalized数据,并存储在Seurat的“@data”slot中。

logfc.data <- logFC(cluster.ids=ser@meta.data$real_id,expr.mat=ser@data)
names(logfc.data)
## [1] "cluster.cells"  "log.fc.cluster"

logFC函数返回一个长度为2的列表,其中包含每个cluster中的细胞以及cluster之间基因的对数倍变化。下一步需要此数据计算enrichment scoresp值。

gse.res <- wmw_gsea(expr.mat=ser@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"

此函数返回两个列表,第一个包含测试统计信息“GSEA_statistics”,第二个包含p值“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 comparisonsres.stats[order(res.stats[,1],decreasing=TRUE)[1:10],] #Top gene sets enriched by z scores
##                                         Group1      Group2    Group3
## HALLMARK_IL2_STAT5_SIGNALING       10.51698615 -7.55892580 -8.113191
## HALLMARK_DNA_REPAIR                10.47159479 -7.62264059 -7.326454
## HALLMARK_COMPLEMENT                10.17982979 -3.68698969 -8.347697
## HALLMARK_UNFOLDED_PROTEIN_RESPONSE  9.54278550 -8.98321228 -7.275756
## HALLMARK_APICAL_SURFACE             7.63384635 -5.70116613  0.000000
## HALLMARK_ANGIOGENESIS               1.32886963 -0.91447563  0.000000
## HALLMARK_HEDGEHOG_SIGNALING         0.77628519  0.66226233  0.000000
## HALLMARK_KRAS_SIGNALING_UP          0.59844933 -0.20261236 -2.365948
## HALLMARK_COAGULATION                0.57787404  9.03118414 -6.988972
## HALLMARK_INFLAMMATORY_RESPONSE      0.05206778 -0.08661822 -1.260939
##                                         Group4
## HALLMARK_IL2_STAT5_SIGNALING        -8.4400626
## HALLMARK_DNA_REPAIR                -10.3905166
## HALLMARK_COMPLEMENT                -12.5780500
## HALLMARK_UNFOLDED_PROTEIN_RESPONSE  -5.8693210
## HALLMARK_APICAL_SURFACE              0.0000000
## HALLMARK_ANGIOGENESIS               -0.6935659
## HALLMARK_HEDGEHOG_SIGNALING          0.3350711
## HALLMARK_KRAS_SIGNALING_UP          -1.1822606
## HALLMARK_COAGULATION                -4.3269207
## HALLMARK_INFLAMMATORY_RESPONSE      -3.3347607
res.pvals[order(res.stats[,1],decreasing=TRUE)[1:10],] #Top gene sets by p values
##                                          Group1       Group2       Group3
## HALLMARK_IL2_STAT5_SIGNALING       2.858190e-24 2.489271e-13 3.451510e-15
## HALLMARK_DNA_REPAIR                2.858190e-24 1.739767e-13 1.286641e-12
## HALLMARK_COMPLEMENT                3.985134e-23 6.949503e-04 6.821492e-16
## HALLMARK_UNFOLDED_PROTEIN_RESPONSE 1.362655e-20 2.577150e-18 1.687982e-12
## HALLMARK_APICAL_SURFACE            1.594960e-13 5.830539e-08 1.000000e+00
## HALLMARK_ANGIOGENESIS              3.003553e-01 5.697704e-01 1.000000e+00
## HALLMARK_HEDGEHOG_SIGNALING        5.642487e-01 7.318339e-01 1.000000e+00
## HALLMARK_KRAS_SIGNALING_UP         6.589077e-01 1.000000e+00 4.406075e-02
## HALLMARK_COAGULATION               6.589077e-01 2.080345e-18 1.130707e-11
## HALLMARK_INFLAMMATORY_RESPONSE     1.000000e+00 1.000000e+00 3.631134e-01
##                                          Group4
## HALLMARK_IL2_STAT5_SIGNALING       1.942653e-16
## HALLMARK_DNA_REPAIR                6.709712e-24
## HALLMARK_COMPLEMENT                1.366256e-34
## HALLMARK_UNFOLDED_PROTEIN_RESPONSE 2.144160e-08
## HALLMARK_APICAL_SURFACE            1.000000e+00
## HALLMARK_ANGIOGENESIS              6.462100e-01
## HALLMARK_HEDGEHOG_SIGNALING        7.856739e-01
## HALLMARK_KRAS_SIGNALING_UP         3.417063e-01
## HALLMARK_COAGULATION               4.939474e-05
## HALLMARK_INFLAMMATORY_RESPONSE     2.460747e-03
names(h.sets)[sets.to.use[1:5]] #Compare to the simulate sets we created
## [1] "HALLMARK_APICAL_SURFACE"
## [2] "HALLMARK_COMPLEMENT"
## [3] "HALLMARK_DNA_REPAIR"
## [4] "HALLMARK_IL2_STAT5_SIGNALING"
## [5] "HALLMARK_UNFOLDED_PROTEIN_RESPONSE"
#Plot the z scores with heatmap3
heatmap3(res.stats,Colv=NA,cexRow=0.5,cexCol=1,scale="row")

图片

恭喜,您已经使用singleseqgset执行了第一个基因集富集分析!

本文来自互联网用户投稿,该文观点仅代表作者本人,不代表本站立场。本站仅提供信息存储空间服务,不拥有所有权,不承担相关法律责任。如若转载,请注明出处:http://www.mzph.cn/pingmian/40579.shtml

如若内容造成侵权/违法违规/事实不符,请联系多彩编程网进行投诉反馈email:809451989@qq.com,一经查实,立即删除!

相关文章

iOS手机竖着拍的照片被旋转了90°的原因以及解决方案

EXIF.getData(IMG_FILE, function () { // IMG_FILE为图像数据 var orientation EXIF.getTag(this, “Orientation”); console.log(“Orientation:” orientation); // 拍照方向 }); 获取拍照方向的结果为1-8的数字&#xff1a; 注意&#xff1a;对于上面的八种方向中&a…

Docker的安装及使用摘要

本文分享一些在docker安装及使用过程中的部分要点&#xff0c;会持续更新&#xff0c;供参考。 1. docker安装 1.1 在ubuntu系统安装 安装指定版本的docker。 # 更新时间&#xff1a;2024年6月23日# docker官方的源无法安装&#xff0c;采用阿里云的源安装docker # 更新软件…

2024.7.4作业

1.梳理笔记(原创) 2. 终端输入一个日期&#xff0c;判断是这一年的第几天 scanf("%d-%d-%d",&y,&m,&d); 闰年2月29天&#xff0c;平年2月28天 #include <stdio.h> int main(int argc,const char *argv[]) { int y0,m0,d0,sum0,i0; …

[论文笔记] pai-megatron-patch Qwen2-72B-CT 后预训练 踩坑记录

经过以下修改,Qwen2-72B-CT可以正常训练,并且benchmark指标和loss正常。 Qwen2-72B-CT开长文本,256卡,16K会OOM,目前能开11K(11008)。 开context parallel需要后续测试。 [论文笔记] Pai-megatron Qwen1.5-14B-CT 后预训练 踩坑记录_pai-megatron-patch 多机-CSDN博客 …

数据库表导出到excel:前置知识1 ALL_TAB_COLS

ALL_TAB_COLS 当前用户可访问的表、视图和群集的列的相关信息 其中几个字段: OWNER&#xff1a;表&#xff0c;视图及群集的Owner   TABLE_NAME&#xff1a; 表&#xff0c;视图及聚簇的名称   COLUMN_NAME&#xff1a; 字段名   DATA_TYPE &#xff1a;字段的数据类型…

python 字典 一个key 多 value 遍历

在Python中&#xff0c;如果一个键对应多个值&#xff0c;你需要确保这些值被存储在一个容器类型&#xff08;如列表或集合&#xff09;中。你可以使用默认字典&#xff08;collections.defaultdict&#xff09;来简化这个过程。以下是一个示例代码&#xff1a; from collecti…

python vtk lod 设置

在Python中使用VTK库设置Level of Detail (LOD)可以通过vtkLODProp3D类来实现。这个类允许你为一个模型指定不同级别的细节表示&#xff0c;从而在渲染时根据模型与摄像机的距离自动切换到更适合的表示。 以下是一个简单的例子&#xff0c;展示如何使用vtkLODProp3D来设置LOD&…

万字长文MySQL Binlog 详细指南

目录 第一阶段 MySQL Binlog 基础用法1. Binlog基本概念1.1 什么是Binlog1.2 Binlog的作用1.3 Binlog格式 2. 配置和管理Binlog2.1 开启Binlog2.2 设置Binlog文件大小和保留时间2.3 查看Binlog状态 3. Binlog的实际应用3.1 数据恢复3.2 主从复制3.3 审计 4. Binlog工具使用4.1 …

收银系统源码-收银台营销功能-购物卡

1. 功能描述 购物卡&#xff1a;基于会员的电子购物卡&#xff0c;支持设置时效、适用门店、以及可用商品&#xff1b;支持售卖和充值赠送&#xff0c;在收银台可以使用&#xff1b; 2.适用场景 会员充值赠送活动&#xff0c;例如会员充值1000元&#xff0c;赠送面值100元购…

LeetCode题解:205. 同构字符串,哈希表,JavaScript,详细注释

原题链接&#xff1a; https://leetcode.cn/problems/isomorphic-strings/ 理解题意&#xff1a; s "foo"和t "bar"&#xff0c;s中的o同时映射了a和r&#xff0c;不正确s "badc"和t "baba"&#xff0c;t中的a同时映射了a和c&am…

145-四路16位125Msps AD FMC子卡模块

一、概述 该板卡可实现4路16bit 125Msps AD 功能&#xff0c;是xilinx开发板设计的标准板卡。FMC连接器是一种高速多pin的互连器件&#xff0c;广泛应用于板卡对接的设备中&#xff0c;特别是在xilinx公司的所有开发板中都使用。该AD&#xff0c;DA子卡模块就专门针对xilinx开发…

pytorch镜像如何通过dockerfile和启动脚本封装tensorboard

一&#xff1a;dockerfile文件内容&#xff0c;这里我们以pytorch/pytorch:1.13.1-cuda11.6-cudnn8-devel基础镜像为例&#xff1a; # 定义基础镜像 FROM pytorch/pytorch:1.13.1-cuda11.6-cudnn8-devel# 设置非互动模式以避免一些安装过程中的对话框 ENV DEBIAN_FRONTENDnoni…

go语言并发编程2-runtime

runtime.Gosched() 作用是让出CPU时间片&#xff0c;重新等待安排任务。执行runtime.Gosched()后&#xff0c;其他协程优先执行&#xff0c;当前所在协程最后执行。 package mainimport ("fmt""runtime" )func main() {go func(s string) {for i : 0; i …

网络爬虫之爬虫逆向的学习途径、相关网站和学习资料

网络爬虫之爬虫逆向的学习途径、相关网站和学习资料 演练和学习网站 CTFTIME 一个全球性的CTF&#xff08;Capture The Flag&#xff09;赛事信息平台&#xff0c;收录了各类CTF比赛。你可以通过参加这些比赛来提升自己的逆向工程和安全技能。 安全客 由360公司运营的安全资讯…

iview 里面的ip 组件封装_iview ipinput

</ul><div v-if"erro_ip" style"color: red;">ip格式错误!</div> </div>最终的效果图如下&#xff1a;![](https://img-blog.csdnimg.cn/20190513170751269.png)最后为了方便大家的沟通与交流请加QQ群&#xff1a; [625787746]( )…

Github 2024-07-03开源项目日报Top10

根据Github Trendings的统计,今日(2024-07-03统计)共有10个项目上榜。根据开发语言中项目的数量,汇总情况如下: 开发语言项目数量JavaScript项目3Jupyter Notebook项目2Python项目2C++项目1Rust项目1TypeScript项目1Vue项目1Go项目1Moby 项目 - 软件容器化的开源工具集 创建…

波动方程 - 波动方程是个什么方程

波动方程 - 波动方程是个什么方程 flyfish 波动方程或称波方程&#xff08;英语&#xff1a;wave equation&#xff09;是一种二阶线性偏微分方程,波动方程是双曲型偏微分方程的最典型代表. 微分方程 微分方程&#xff08;Differential Equation&#xff09;是一类包含未知…

C++语言特性层(Language Features Layer)

1.语言基础 &#xff08;1&#xff09;指针 定义&#xff1a; 指针是一个变量&#xff0c;用于存储另一个变量的内存地址。 特性&#xff1a; 可变性&#xff1a;指针可以重新指向不同的变量。空指针&#xff1a;指针可以为空&#xff08;即指向 nullptr&#xff09;。大小&am…

羊大师:羊奶养生,解锁健康之道的新密码

在探寻健康与养生的旅途中&#xff0c;我们总渴望找到那把开启健康之门的钥匙。而今&#xff0c;羊奶以其独特的营养价值和健康益处&#xff0c;正悄然成为那把解锁健康之道的新密码。 羊奶&#xff0c;自古以来便是自然赋予的珍贵礼物。它富含优质蛋白、多种维生素及矿物质&am…

nginx的重定向(rewrite)

1、location 匹配 location匹配的就是后面的URL&#xff0c;对访问的路径做访问控制或者代理转发 共有三个匹配&#xff1a;精确匹配、正则匹配、一般配 a、精确匹配 格式&#xff1a;location/ 对字符串进行完全匹配&#xff0c;必须完全合 c、正则匹配 ^~&#xff1a;前…