基因集合打分不考虑基因高低表达 基因集合评分不考虑高低表达 绝对值 基因集合里面的正负表达影响addmodulescore的结果 构建的基因集进行评分时 基因的高低正负表达基因集合评分不考虑高低表达_YoungLeelight的博客-CSDN博客

单细胞数据分析之缺氧评分 - 腾讯云开发者社区-腾讯云

试一试这个R包做单细胞分析-scMetablism - 腾讯云开发者社区-腾讯云

偶尔逛朋友圈发现一年前跟着我们生信技能树学生信的研究生开发了自己的单细胞数据分析相关R包,4(热图,气泡图,upset图,堆叠条形图)+4(密度散点图,半小提琴,山峦图,密度热图)美图吸引了我的注意力,果断邀稿,希望可以介绍他的R包使用方法,以及开发新的体会!

下面范垂钦的投稿

一、背景

  1. 在单细胞转录组测序中,整合多样本数据进行分析逐渐成为一种趋势。
  2. 越来越多的研究者趋于使用未校正批次效应后的数据进行差异基因分析。
  3. 他们认为,校正批次效应后的数据,会掩盖部分真实的生物学差异。但校正批次效应后的数据是否能用于基因集富集分析,以及样本之间的批次效应是否会影响基因富集分析结果仍然是一个争论。
  4. 因此,我们重新审视现有的、主流的基因集富集方法,希望能找到受批次效应影响较小的基因集富集分析方法
  5. 在这里,我们重新审视了9种常见的功能集打分方法:GSEA、GSVA、PLAGE、Zscore、AddModuleScore、ssGSEA、AUCell、UCell和singscore

二、审视结果

  1. GSEA:执行前需对所有样本进行分组,然后基于分组去计算得到排序基因列表。这个过程中,我们需要考虑不同分组中样本构成的影响;
  2. GSVA:首先需要对所有样本中每个基因进行累积分布密度函数的核估计。这个过程中,需要考虑所有样本,容易受到样本背景信息的影响
  3. PLAGE 和 z-score:首先需要对基因表达矩阵执行标准化处理。同样的,这个过程容易受样本构成的影响;
  4. AddModuleScore:Seurat包中的AddModuleScore函数,需要先计算基因集中所有基因的平均值,再根据平均值把表达矩阵切割成若干份,然后从切割后的每一份中随机抽取对照基因(基因集外的基因)作为背景值。因此,在整合不同样本的情况下,即使使用相同基因集为相同细胞打分,也会产生不同的富集评分;
  5. AUCell:基于单个样本中的基因表达排名(gene expression rank),使用曲线下面积来评估输入基因集是否在单个样本的前5%表达基因内富集;
  6. UCell:基于单个样本的基因表达排名,使用Mann-Whitney U统计量计算单个样本的基因集富集评分;7.singscore:基于单个样本的基因表达排名,评估基因集远离中心的程度从而计算基因集富集评分;8.ssgsea:基于单个样本的基因表达排名,通过计算单个样本中基因集内和基因集外的经验累积分布函数之间的差值进而生成富集分数。然后,根据基因表达谱最大表达值和最小表达值的差值对富集分数进行标准化。这一步容易受样本构成的影响。

三、原理解析

1.纳入合适的方法:

我们淘汰了所有需要考虑样本组成的基因集富集分析方法,选取了基于单个样本的基因表达排名的基因集分析方法:AUCell、UCell和singscore。同时,我们也部分改进了ssGSEA,取消最后的标准化步骤,使之更贴近于单个样本的基因集富集分析。

2.构建基因集:

为了方便用户获取MSigDB数据库中预先定义好的基因集,我们内置了MSigDB包进行基因集的获取。同时,我们也支持多个物种的基因集获取,以及多种基因格式的表达矩阵的输入。除此之外,我们还允许用户自定义自己的基因集进行打分计算如果用户的基因集中既包括正向的基因,也包括负向的基因。该基因集将默认只使用UCell包和singscore包进行打分处理

3.输入对象和数据清洗:

我们允许直接输入单细胞表达矩阵或者Seurat对象。我们内置了Seurat包,可以将多种基因集的富集分数矩阵直接保存到Seurat对象中。同时,我们也支持过滤单细胞表达矩阵中所有细胞表达量为0的基因。当然,用户也可以自定义自己的过滤标准。合适的过滤指标可以改善富集分析的结果

4.差异分析和综合评估:

为了评估基因集在某个细胞亚群中是否富集我们通过多种基因集富集方法分别对单个细胞进行打分,并生成多个基因集富集分数矩阵。接着,我们通过wilcox检验计算各个基因集富集分数矩阵中每个细胞亚群差异表达的基因集单一的基因集富集分析方法不仅只能反映有限的信息,而且也容易带来误差。我们期待从多个角度解释复杂的生物学问题,并找到生物学问题中的共性部分。简单地为多种基因集富集分析方法的结果取共同交集,不仅容易得到少而保守的结果,而且忽略了富集分析方法中很多的其他信息,例如不同基因集的相对富集程度信息。单纯地取共同交集无法体现基因集的相对富集程度,而我们希望目标基因集在大部分富集分析方法中都是富集且富集程度没有明显差异。因此,我们通过RobustRankAggreg包中的秩聚合算法(robust rank aggregation, RRA)对差异分析的结果进行综合评估,筛选出在大部分基因集富集分析方法中都显著富集的基因集

5.可视化展示:

我们内置了多种可视化的函数,不仅允许用户通过热图、气泡图、柱状图和upset图展示它们综合结果,而且允许用户通过密度散点图、半小提琴图、山峦图和密度热图展示目标基因集在具体富集分析方法中的表达水平和数据分布。其中,热图、upset图、密度热图主要由ComplexHeatmap包生成;气泡图和柱状图主要由ggplot2包、ggtree包、aplot包生成;密度散点图主要由Seurat包和Nebulosa包生成,半小提琴图主要由Seurat包和gghalves包生成;山峦图主要由Seurat包和ggridges包生成。最后,为了方便用户将可视化结果与其他的ggplot2对象进行拼图操作,我们也通过ggplotify包把输出结果转换为ggplot2对象。最后,为了方便用户的实现,我们开发了R包irGSEA (https://github.com/chuiqin/irGSEA/),集成了上述工作流程。

四、使用教程

1.安装前的准备

在安装irGSEA 需要查看自己是否安装了以下这些包

# install packages from CRAN
cran.packages <- c("msigdbr", "dplyr", "purrr", "stringr","magrittr","RobustRankAggreg", "tibble", "reshape2", "ggsci","tidyr", "aplot", "ggfun", "ggplotify", "ggridges","gghalves", "Seurat", "SeuratObject", "methods","devtools", "BiocManager","data.table","doParallel","doRNG")
if (!requireNamespace(cran.packages, quietly = TRUE)) {install.packages(cran.packages, ask = F, update = F)
}# install packages from Bioconductor
bioconductor.packages <- c("GSEABase", "AUCell", "SummarizedExperiment","singscore", "GSVA", "ComplexHeatmap", "ggtree","Nebulosa")
if (!requireNamespace(bioconductor.packages, quietly = TRUE)) {BiocManager::install(bioconductor.packages, ask = F, update = F)
}
# install packages from Github
if (!requireNamespace("UCell", quietly = TRUE)) {devtools::install_github("carmonalab/UCell")
}
if (!requireNamespace("irGSEA", quietly = TRUE)) {devtools::install_github("chuiqin/irGSEA")
}

复制

2.加载示例数据

我们通过SeuratData包加载示例数据集(注释好的PBMC数据集)

# devtools::install_github('satijalab/seurat-data')
library(SeuratData)
# view all available datasets
View(AvailableData())
# download 3k PBMCs from 10X Genomics
# 示例数据有点大,建议在网速好的时候下载,如果失败了,请重启R后多尝试几遍
InstallData("pbmc3k")
# the details of pbmc3k.final
?pbmc3k.final

复制

library(Seurat)
library(SeuratData)
# loading dataset
data("pbmc3k.final")
pbmc3k.final <- UpdateSeuratObject(pbmc3k.final)
# plot
DimPlot(pbmc3k.final, reduction = "umap",group.by = "seurat_annotations",label = T) + NoLegend()

复制

PBMC的UMAP图

顺手设置一下数据集的ident

# set cluster to idents
Idents(pbmc3k.final) <- pbmc3k.final$seurat_annotations

复制

3.加载R包

这一步出错的话,要看一下前面的包有没有装好

library(UCell)
library(irGSEA)

复制

4.计算富集分数

当你的ncore设置大于1的时候,发生下面的错误:Error (Valid ‘mctype’: ‘snow’ or ‘doMC’),你应该检查一下你的AUCell 版本,确保版本大于等于1.14 。如果你比较懒,那你直接把ncore设置为1也是可以的,只是运行速度会稍微慢一点。

pbmc3k.final <- irGSEA.score(object = pbmc3k.final, assay = "RNA",slot = "data", seeds = 123, ncores = 1,min.cells = 3, min.feature = 0,custom = F, geneset = NULL, msigdb = T,species = "Homo sapiens", category = "H",  subcategory = NULL, geneid = "symbol",method = c("AUCell", "UCell", "singscore","ssgsea"),aucell.MaxRank = NULL, ucell.MaxRank = NULL,kcdf = 'Gaussian')
#> Validating object structure
#> Updating object slots
#> Ensuring keys are in the proper strucutre
#> Ensuring feature names don't have underscores or pipes
#> Object representation is consistent with the most current Seurat version
#> Calculate AUCell scores
#> Warning: Feature names cannot have underscores ('_'), replacing with dashes
#> ('-')#> Warning in .filterFeatures(expr, method): Feature names cannot have underscores
#> ('_'), replacing with dashes ('-')
#> Finish calculate ssgsea scores# 返回一个Seurat对象,富集分数矩阵存放在RNA外的assay中
Seurat::Assays(pbmc3k.final)
#> [1] "RNA"       "AUCell"    "UCell"     "singscore" "ssgsea"

复制

5.整合差异基因集

# Wlicox test is perform to all enrichment score matrixes and gene sets
# with adjusted p value &lt; 0.05 are used to integrated through RRA.
# Among them, Gene sets with p value &lt; 0.05 are statistically
# significant and common differential in all gene sets enrichment analysis
# methods. All results are saved in a list.
result.dge <- irGSEA.integrate(object = pbmc3k.final,group.by = "seurat_annotations",metadata = NULL, col.name = NULL,method = c("AUCell","UCell","singscore","ssgsea"))
#> Calculate differential gene set : AUCell
#> Calculate differential gene set : UCell
#> Calculate differential gene set : singscore
#> Calculate differential gene set : ssgsea
class(result.dge)
#> [1] "list"

复制

6.可视化展示

1). 全局展示

①.热图

热图展示了综合评价中具体基因集在每个细胞亚群是否具有统计学意义差异;其中,浅蓝色的格子无统计学差异,红色的格子具有统计学差异。格子中的星号越多,格子的P值越小;左边的聚类树代表不同基因集在不同细胞亚群中表达模式的相似性;上方的条形图分别代表不同的细胞亚群,以及差异基因集在细胞亚群中是呈现上调还是下调趋势;你还可以把method从'RRA"换成“ssgsea”,展示特定基因集富集分析方法中差异上调或差异下调的基因集;

irGSEA.heatmap.plot <- irGSEA.heatmap(object = result.dge,method = "RRA",top = 50,show.geneset = NULL)
irGSEA.heatmap.plot

复制

热图

②.气泡图

气泡图展示了综合评价中具体基因集在每个细胞亚群是否具有统计学意义差异;其中,浅蓝色的点无统计学差异,红色的点具有统计学差异。点越大,P值越小;左边的聚类树代表不同基因集在不同细胞亚群中表达模式的相似性;上方的条形图分别代表不同的细胞亚群,以及差异基因集在目标细胞亚群中是呈现上调还是下调趋势;你还可以把method从'RRA"换成“ssgsea”,展示特定基因集富集分析方法中差异上调或差异下调的基因集;如果运行中提示“ error (argument “caller_env” is missing, with no default)”,请卸载了ggtree包,并重新安装remotes::install_github(”YuLab-SMU/ggtree“)

irGSEA.bubble.plot <- irGSEA.bubble(object = result.dge,method = "RRA",top = 50)
irGSEA.bubble.plot

复制

气泡图

③.upset plot

upset图展示了综合评估中每个细胞亚群具有统计学意义差异的基因集的数目,以及不同细胞亚群之间具有交集的差异基因集数目;左边不同颜色的条形图代表不同的细胞亚群;上方的条形图代表具有交集的差异基因集的数目;中间的气泡图单个点代表单个细胞亚群,多个点连线代表多个细胞亚群取交集.;如果运行代码发生以下warning“the condition has length > 1 and only the first element will be used”。这个警告可忽略。

irGSEA.upset.plot <- irGSEA.upset(object = result.dge,method = "RRA")
#> Warning in if (as.character(ta_call[[1]]) == "upset_top_annotation") {: the
#> condition has length > 1 and only the first element will be used
irGSEA.upset.plot

复制

upset图

④.堆叠条形图

堆叠柱状图具体展示每种基因集富集分析方法中每种细胞亚群中上调、下调和没有统计学差异的基因集数目;上方的条形代表每个亚群中不同方法中差异的基因数目,红色代表上调的差异基因集,蓝色代表下调的差异基因集;中间的柱形图代表每个亚群中不同方法中上调、下调和没有统计学意义的基因集的比例;

irGSEA.barplot.plot <- irGSEA.barplot(object = result.dge,method = c("AUCell", "UCell", "singscore","ssgsea"))
irGSEA.barplot.plot

复制

image.png

2). 局部展示

①.密度散点图

密度散点图将基因集的富集分数和细胞亚群在低维空间的投影结合起来,展示了特定基因集在空间上的表达水平。其中,颜色越黄,代表富集分数越高;

scatterplot <- irGSEA.density.scatterplot(object = pbmc3k.final,method = "UCell",show.geneset = "HALLMARK-INFLAMMATORY-RESPONSE",reduction = "umap")
scatterplot

复制

image.png

②.半小提琴图

半小提琴图同时以小提琴图(左边)和箱线图(右边)进行展示。不同颜色代表不同的细胞亚群;

halfvlnplot <- irGSEA.halfvlnplot(object = pbmc3k.final,method = "UCell",show.geneset = "HALLMARK-INFLAMMATORY-RESPONSE")
halfvlnplot

复制

image.png

③.山峦图

山峦图中上方的核密度曲线展示了数据的主要分布,下方的条形编码图展示了细胞亚群具体的数量。不同颜色代表不同的细胞亚群,而横坐标代表不同的表达水平;

ridgeplot <- irGSEA.ridgeplot(object = pbmc3k.final,method = "UCell",show.geneset = "HALLMARK-INFLAMMATORY-RESPONSE")
ridgeplot
#> Picking joint bandwidth of 0.00533

复制

image.png

④.密度热图

密度热图展示了具体差异基因在不同细胞亚群中的表达和分布水平。颜色越红,代表富集分数越高;

densityheatmap <- irGSEA.densityheatmap(object = pbmc3k.final,method = "UCell",show.geneset = "HALLMARK-INFLAMMATORY-RESPONSE")
densityheatmap

复制

image.png

五、拓展:自己如何做一个R包

如果这篇推文只到教程结束,那一切都将索然无味。

我也没有必要在前面水那么多字数的背景,并详细介绍每一步我是怎么做的。

因为我的单细胞R包框架是可以被模仿的,并且只需要以下3步

  • 1.针对一个生物学问题收集多种解决方法;
  • 2.对多种解决方法的结果进行一个综合评估;
  • 3.对综合评估的结果进行可视化展示;

举个例子,你可以做一个计算细胞亚群差异基因的单细胞R包

  • 1.收集多种单细胞差异基因集计算方法(例如wilicox,AUC,MAST等);
  • 2.分别计算,每种方法中每种细胞亚群中差异表达的基因,用一种综合评估的方法(例如RRA算法)评估这些差异分析方法的结果,找出在大部分方法中都具有差异的基因;
  • 3.寻找合适的可视化方法(例如热图、散点图、气泡图等)可视化展示综合评估的结果;

再举个例子,你可以做一个计算TCGA中不同分子亚型免疫浸润差异的R包

  • 1.收集多种免疫浸润的计算方法(例如cibersortX,estimate,xcell,ssgsea,等);
  • 2.分别计算,每种方法中每种分子亚型中差异表达的免疫细胞群,用一种综合评估的方法(例如RRA算法)评估这些差异分析方法的结果,找出在大部分方法中都具有差异的基因;
  • 3.寻找合适的可视化方法(例如热图、散点图、气泡图等)可视化展示综合评估的结果;

再再再举个例子,收集多种单细胞转录因子调控网络的计算方法(例如SCENIC,DoRothEA等),收集多种转录组差异分析的方法(例如limma-voom,edgeR, DESeq2等)...然后按照上文介绍的流程,如法炮制即可。

简单总结一下这个套路:多种方法 + 综合评估 + 可视化 = R包

如果你已经厌倦了GEO/TCGA的数据挖掘,可以尝试一下开发一两个有趣的R包。

当然,只有框架还是不够,你还需要有一定的R语言基础,以及R包开发的经验。

关于R包开发,推荐你们去看生信技能树在B站的视频,以及Hadley Wickham写的《R包开发》这本书,相信会让你减少很多摸索的时间。在此,我也要感谢一下生信技能树组织的《R包开发》课程,以及永和大神的指点。

评分addmodule 绝对值评分 8种方法可视化你的单细胞基因集打分gsea 缺氧评分相关推荐

  1. JS实现星星评分功能实例代码(两种方法)

    转载自   JS实现星星评分功能实例代码(两种方法) 一.方法1 1.用到图片 2.结构和样式 <!DOCTYPE html> <html lang="en"&g ...

  2. java实现星级评分功能_JS实现星星评分功能实例代码(两种方法)

    一.方法1 1.用到图片 2.结构和样式 Document ul { padding-left: 0; overflow: hidden; } ul li { float: left; list-st ...

  3. Python学习之求绝对值的几种方法

    import mathdef abs_value1():a = float(input('1.请输入一个数字:'))if a >= 0:a = aelse:a = -aprint('绝对值为:% ...

  4. java中获取绝对值的方法_Java完美判断绝对值的两种方法 | 彬菌

    版权声明:转载原创文章请以超链接形式请注明原文章出处,尊重作者,尊重原创! 恰饭广告 if-else语句判断: import java.util.Scanner; public class Absol ...

  5. 统计图的连通块的个数的两种方法

    @算法学习 两种方法 DFS遍历法 并查集法 1. DFS遍历计算连通块 先上代码: #include <stdio.h> #include <vector>using nam ...

  6. Matlab 计算均方误差MSE的三种方法

    Matlab 计算均方误差MSE的三种方法 数据说明: ytest 测试集y,真实的y值,是一维数组: ytest_fit 基于测试集 x 预测的y值,是一维数组: test_error 是预测误差. ...

  7. 【Elasticsearch】Elasticsearch自定义评分的N种方法

    1.概述 首先参考文章:[Elasticsearch]Elasticsearch 相关度评分 TF&IDF 然后转载文章:实战 | Elasticsearch自定义评分的N种方法 2.三个问题 ...

  8. python与excel做数据可视化-用Python进行数据可视化的10种方法

    原标题:用Python进行数据可视化的10种方法 2015-11-19 关于转载授权 大数据文摘作品,欢迎个人转发朋友圈,自媒体.媒体.机构转载务必申请授权,后台留言"机构名称+转载&quo ...

  9. 带圆圈大小的散点图_Python数据可视化,Matplotlib绘制“散点图”的两种方法!...

    前言 散点图是Matplotlib常用图形之一,与线形图类似.但是这种图形不再由线段连接,而是由独立的点.圆圈或其他形状构成.那么怎么画散点图呢?Matplotlib给出了两种不同的方法,去画散点图. ...

最新文章

  1. SQL Server中查看SQL句子执行所用的时间
  2. Windows驱动开发 - 设备对象初步学习
  3. vs代码补全的快捷键_一款Python编程的自动补全插件神器——kite
  4. ubuntu让/etc/hosts修改后立刻生效
  5. webpack 异步加载配置文件_详解webpack异步加载业务模块
  6. 【飞秋】位运算与组合搜索(二)
  7. Docker快速安装ZooKeeper开源分布式协调服务器
  8. 系统提示服务器响应错误,Win10系统无法打开软件提示“服务器没有及时响应或控制请求”错误的解决方法...
  9. JIRA中设置[描述]字段的默认值
  10. 这个国庆,我去佛山看舞狮,太惊艳!
  11. H5前端实现微信分享(处理二次分享问题)
  12. 集线器、交换机、路由器、中继器及网关、网桥之间的区别
  13. vue动态调节背景图片
  14. 三峡学院计算机调剂,2018年重庆三峡学院考研预调剂公告
  15. 程序员薪酬到底有多高?来看硅谷的工程师统计
  16. 【新闻推荐系统】(task3)自动化构建用户及物料画像
  17. Tensorflow2.10 Object Detetcion安装教程
  18. Ubuntu测试使用速腾RS-Lidar-16
  19. mysql备份恢复项目_mysql备份恢复之xtrabackup (XBK、Xbackup)
  20. 学校学生计算机教室解决方案,学校学生计算机教室解决方案设计.docx

热门文章

  1. 故障频出的摩拜单车,背后是野蛮生长的原罪
  2. diff算法中的概念
  3. 适合 Kubernetes 初学者的一些实战练习(二)
  4. Linux命令jar包操作
  5. Linux 运行jar包命令如下:
  6. “华为手机”和“荣耀手机”的区别?今天全都告诉你
  7. 实拍当贝F5与坚果J10S,莫要听信网上水军
  8. (二)SQL大小写规范和 sql_mode
  9. ZBLOG 翻译插件
  10. CentOS7.9安装twemproxy,实现redis集群