AUCell和AddModuleScore函数进行基因集评分

2024-06-22 17:44

本文主要是介绍AUCell和AddModuleScore函数进行基因集评分,希望对大家解决编程问题提供一定的参考价值,需要的开发者们随着小编来一起学习吧!

AUCellAddModuleScore 分析是两种主流的用于单细胞RNA测序数据的基因集活性分析的方法。这些基因集可以来自文献、数据库或者根据具体研究问题进行自行定义。

AUCell分析原理:

1、AUCell分析可以将细胞中的所有基因按表达量进行排序,生成一个基因排名列表,表达量越高的基因排名越靠前。

2、接下来对每个基因集中的基因找到它们在每个细胞的基因排名列表中的位置,这些位置则反映了基因集内基因在特定细胞群中的表达情况。

3、计算基因集在每个细胞中的活性评分。基于基因集内基因的排名,通过计算这些基因排名的累计面积(AUC,Area Under the Curve)来评估基因集的活性。AUC值越大,表明该基因集在该细胞中的表达越活跃。

AddModuleScore分析原理:

1、计算每个基因集中的基因在每个细胞中的平均表达值

2、选择与目标基因集大小相似的背景基因集,通过目标基因集的平均表达值减去背景基因集的平均表达值,得到基因模块评分(这个评分代表)。这个背景基因集是有函数内部把表达矩阵自行切割若干份之后随机抽取每一份中的基因作为背景基因集。

应用场景
  • 细胞类型鉴定:识别不同细胞类型或细胞亚群的特征基因集活性。

  • 功能状态分析:分析细胞的功能状态,例如细胞周期、免疫反应等。

  • 生物学过程探索:探索特定生物学过程中基因集的表达活性。

AUCell分析步骤:

1.读取基因集数据(采用了Autophagy基因数据)

rm(list = ls())
autophagy_genes <- read.xlsx("Autophagy.xlsx",colNames = T) 
g <- autophagy_genes$Symbol
head(autophagy_genes)
#GeneId                                                                    Name Symbol
#1  55626                                          autophagy/beclin-1 regulator 1 AMBRA1
#2   8542                                                     apolipoprotein L, 1  APOL1
#3    405                          aryl hydrocarbon receptor nuclear translocator   ARNT
#4    410                                                         arylsulfatase A   ARSA
#5    411                                                         arylsulfatase B   ARSB
#6    468 activating transcription factor 4 (tax-responsive enhancer element B67)   ATF4

2.读取R包和需要分析的数据(采用了单细胞pbmc数据)

library(Seurat)
library(tidyverse)
library(openxlsx)
load("step1.final.Rdata") #pbmc数据
sce <- step1.final
#check一下
DimPlot(sce,group.by = "celltype",label = T)+ NoLegend() +ggsci::scale_color_d3()

3、AUCell分析

library(GSEABase)
geneSets <- GeneSet(g, setName="autophagy")
geneSets
#BiocManager::install("AUCell")
library(AUCell)
exprMatrix = sce@assays$RNA@data
rownames(exprMatrix) = Features(sce)
colnames(exprMatrix) = Cells(sce)
#对矩阵中的每个细胞里,给基因进行排序(可见下图1)
cells_rankings <- AUCell_buildRankings(exprMatrix,, plotStats=TRUE) 
#把每个细胞中的基因表达量从高到底排列并计算数量
#Quantiles for the number of genes detected by cell: 
#(Non-detected genes are shuffled at the end of the ranking. Keep it in mind when choosing #the threshold for calculating the AUC).
#    min      1%      5%     10%     50%    100% 
# 491.00  633.15  806.00  897.30 1323.00 5418.00 
cells_AUC <- AUCell_calcAUC(geneSets, cells_rankings,aucMaxRank = nrow(cells_rankings)*0.05)
#为了计算AUC,默认情况下只使用排名中前5%的基因(即检查基因集中的基因是否在前5%之内)。
#这样可以在更大的数据集上更快地执行,并减少排名底部噪音的影响(例如,许多基因可能表达量为0)。要考虑的百分比可以通过参数 aucMaxRank 参数进行修改。可以通过AUCell_buildRankings提供的直方图进行辅助判断。
#Genes in the gene sets NOT available in the dataset: 
# autophagy:  26 (12% of 222)
set.seed(123)
#见下图2
cells_assignment <- AUCell_exploreThresholds(cells_AUC, plotHist=TRUE, assign=TRUE) 
auc_thr = cells_assignment$HMMR$aucThr$selected 
auc_thr

不同细胞中基因表达情况呈偏态分布

这里重点提一下如何看这个AUC histogram,我这边采用autophagy基因集在pbmc数据集中进行分析发现,AUC大于0.11的细胞为活性细胞,但是pbmc中没有细胞符合要求~ 这也提示了我们在做AUCell分析前,需要仔细考虑纳入分析的基因集和单细胞数据是否”合适“。

接下来我以AUCell github上的资料为例子,这些AUC柱状图是对“自定义/活性基因集的细胞”的直观展示。开发者采用了神经细胞的数据集进行分析和展示。

1) 图片的标题是指不同的细胞亚群和基因数量,比如Astrocyte_Cahoy (526g),星形胶纸细胞,代表这群细胞在研究者纳入分析的数据集中存在的基因数量为526个。其中Random(50g), 代表研究者随机提取细胞和基因。

2) X轴代表了AUC值的大小,Y轴代表了不同AUC值下的细胞数量,图形里边的AUC值代表了阈值,AUC阈值下边的具体数值代表了达到要求的“活性”细胞数量。

3) 理想的情况图形呈双峰分布,数据集中大多数细胞是呈现较低的AUC值,而少部分细胞则呈现较高的AUC值。比如Oligodendrocyte_Cahoy分析结果。比较类似的结果是Neuron和Microglia_lavin,虽然它们的基因集细胞数量很少或者符合要求的细胞数很少,但结果仍旧呈现出了双峰形态,或许侧面说明了他们的基因集或者筛除的细胞能够明显的表征某些特性。

4) 如果基因集是数据集中的高比例细胞的标记(比如 Neuron结果),那么分布图形可能类似于管家基因的分析结果(HK-like)。

5) 基因集的大小会影响结果。更大的基因集(100-2k)可会使结果更稳定,更易评估,随着基因数目的减少,AUC = 0的细胞数目可能增加。当然如果选定的基因集非常给力,那么也可能呈现较好的结果 (即Neuron_Lein结果)。

4、图形可视化

#由于我这里使用的autophagy基因并不能很好的区分高低活性细胞,因此我更换了数据集
#数据集采用了GSE162025鼻咽癌数据集中的上皮细胞
#基因集采用了"SCNM1","CDC42SE1","ZNF687" (随意指定)
Idents(sce) <- sce$celltype
sce$auc_score = as.numeric(getAUC(cells_AUC))
sce$auc_group = ifelse(sce$auc_score>auc_thr,"high_A","low_A") #自行修改dat<- data.frame(sce@meta.data, sce@reductions$umap@cell.embeddings,seurat_annotation = sce@active.ident)
class_avg <- dat %>%group_by(seurat_annotation) %>% #按照seurat_annotation列(即细胞的分类)对数据进行分组。summarise(umap_1 = median(umap_1),umap_2 = median(umap_2) #对每个分组计算UMAP坐标的中位数 画label)library(ggpubr)
ggplot(dat, aes(umap_1, umap_2))  +geom_point(aes(colour  = auc_score)) +viridis::scale_color_viridis(option="A") +ggrepel::geom_label_repel(aes(label = seurat_annotation),data = class_avg,label.size = 0,segment.color = NA)+theme_bw()
ggsave(filename="Aucell.pdf",width = 12,height = 7)# 高低分组
ggplot(dat,aes(x = umap_1,y = umap_2))+geom_point(aes(color = auc_group),size = 0.5)+theme_classic()
ggsave(filename="Aucell2.pdf",width = 8,height = 7)

然后接下来可以根据基因集的高低分组进行更多的个性化分析啦~

AddModuleScore分析步骤:

1、读取数据和基因集

rm(list = ls())
library(Seurat)
g <- c("SCNM1","CDC42SE1","ZNF687")
load("sce_epi.Rdata")  #数据集采用了GSE162025鼻咽癌数据集中的上皮细胞
#DimPlot(sce,group.by = "celltype",label = T)+ NoLegend() +ggsci::scale_color_d3()

2、AddModuleScore分析

#AddmoduleScore函数是seraut包自带的很方便
sce =  AddModuleScore(object = sce,features = g)
colnames(sce@meta.data)
p =FeaturePlot(sce,'Cluster1') #默认名称cluster1
p dat<- data.frame(sce@meta.data, sce@reductions$umap@cell.embeddings,seurat_annotation = sce@active.ident)
class_avg <- dat %>%group_by(seurat_annotation) %>% #按照seurat_annotation列(即细胞的分类)对数据进行分组。summarise(umap_1 = median(umap_1),umap_2 = median(umap_2) #对每个分组计算UMAP坐标的中位数 画label)library(ggpubr)
ggplot(dat, aes(umap_1, umap_2))  +geom_point(aes(colour  = Cluster1)) + #修改这里的colourviridis::scale_color_viridis(option="A") +ggrepel::geom_label_repel(aes(label = seurat_annotation),data = class_avg,label.size = 0,segment.color = NA)+theme_bw()
ggsave(filename="Aucell-addmodulescore.pdf",width = 12,height = 7)

获得cluster评分之后就可以按照中位值或者其他的值进行分组,从而进行后续的个性化分析啦~

:若对内容有疑惑或者有发现明确错误的朋友,请联系后台(希望多多交流)。更多内容可关注公众号:生信方舟

- END -

这篇关于AUCell和AddModuleScore函数进行基因集评分的文章就介绍到这儿,希望我们推荐的文章对编程师们有所帮助!



http://www.chinasem.cn/article/1084980

相关文章

hdu1171(母函数或多重背包)

题意:把物品分成两份,使得价值最接近 可以用背包,或者是母函数来解,母函数(1 + x^v+x^2v+.....+x^num*v)(1 + x^v+x^2v+.....+x^num*v)(1 + x^v+x^2v+.....+x^num*v) 其中指数为价值,每一项的数目为(该物品数+1)个 代码如下: #include<iostream>#include<algorithm>

【Prometheus】PromQL向量匹配实现不同标签的向量数据进行运算

✨✨ 欢迎大家来到景天科技苑✨✨ 🎈🎈 养成好习惯,先赞后看哦~🎈🎈 🏆 作者简介:景天科技苑 🏆《头衔》:大厂架构师,华为云开发者社区专家博主,阿里云开发者社区专家博主,CSDN全栈领域优质创作者,掘金优秀博主,51CTO博客专家等。 🏆《博客》:Python全栈,前后端开发,小程序开发,人工智能,js逆向,App逆向,网络系统安全,数据分析,Django,fastapi

业务中14个需要进行A/B测试的时刻[信息图]

在本指南中,我们将全面了解有关 A/B测试 的所有内容。 我们将介绍不同类型的A/B测试,如何有效地规划和启动测试,如何评估测试是否成功,您应该关注哪些指标,多年来我们发现的常见错误等等。 什么是A/B测试? A/B测试(有时称为“分割测试”)是一种实验类型,其中您创建两种或多种内容变体——如登录页面、电子邮件或广告——并将它们显示给不同的受众群体,以查看哪一种效果最好。 本质上,A/B测

C++操作符重载实例(独立函数)

C++操作符重载实例,我们把坐标值CVector的加法进行重载,计算c3=c1+c2时,也就是计算x3=x1+x2,y3=y1+y2,今天我们以独立函数的方式重载操作符+(加号),以下是C++代码: c1802.cpp源代码: D:\YcjWork\CppTour>vim c1802.cpp #include <iostream>using namespace std;/*** 以独立函数

函数式编程思想

我们经常会用到各种各样的编程思想,例如面向过程、面向对象。不过笔者在该博客简单介绍一下函数式编程思想. 如果对函数式编程思想进行概括,就是f(x) = na(x) , y=uf(x)…至于其他的编程思想,可能是y=a(x)+b(x)+c(x)…,也有可能是y=f(x)=f(x)/a + f(x)/b+f(x)/c… 面向过程的指令式编程 面向过程,简单理解就是y=a(x)+b(x)+c(x)

遮罩,在指定元素上进行遮罩

废话不多说,直接上代码: ps:依赖 jquer.js 1.首先,定义一个 Overlay.js  代码如下: /*遮罩 Overlay js 对象*/function Overlay(options){//{targetId:'',viewHtml:'',viewWidth:'',viewHeight:''}try{this.state=false;//遮罩状态 true 激活,f

利用matlab bar函数绘制较为复杂的柱状图,并在图中进行适当标注

示例代码和结果如下:小疑问:如何自动选择合适的坐标位置对柱状图的数值大小进行标注?😂 clear; close all;x = 1:3;aa=[28.6321521955954 26.2453660695847 21.69102348512086.93747104431360 6.25442246899816 3.342835958564245.51365061796319 4.87

OpenCV结构分析与形状描述符(11)椭圆拟合函数fitEllipse()的使用

操作系统:ubuntu22.04 OpenCV版本:OpenCV4.9 IDE:Visual Studio Code 编程语言:C++11 算法描述 围绕一组2D点拟合一个椭圆。 该函数计算出一个椭圆,该椭圆在最小二乘意义上最好地拟合一组2D点。它返回一个内切椭圆的旋转矩形。使用了由[90]描述的第一个算法。开发者应该注意,由于数据点靠近包含的 Mat 元素的边界,返回的椭圆/旋转矩形数据

Unity3D 运动之Move函数和translate

CharacterController.Move 移动 function Move (motion : Vector3) : CollisionFlags Description描述 A more complex move function taking absolute movement deltas. 一个更加复杂的运动函数,每次都绝对运动。 Attempts to

✨机器学习笔记(二)—— 线性回归、代价函数、梯度下降

1️⃣线性回归(linear regression) f w , b ( x ) = w x + b f_{w,b}(x) = wx + b fw,b​(x)=wx+b 🎈A linear regression model predicting house prices: 如图是机器学习通过监督学习运用线性回归模型来预测房价的例子,当房屋大小为1250 f e e t 2 feet^