单细胞转录组:差异基因分析和富集分析 - 教程

news/2025/10/20 15:04:38/文章来源:https://www.cnblogs.com/wzzkaifa/p/19152689

一、单细胞与普通转录组差异基因分析的区别

单细胞转录组(scRNA-seq)与普通转录组(bulk RNA-seq)在差异基因分析方面存在本质区别:

1. 分析单位不同:

  • 普通转录组:以组织或细胞混合物为单位,结果反映的是平均表达水平

  • 单细胞转录组:以单个细胞为单位,能够揭示细胞间的异质性

2. 数据特性不同:

  • 普通转录组:信号强度高,噪音相对较低

  • 单细胞转录组:存在大量零值(dropout),数据稀疏性高,技术噪音大

3. 分析方法不同:

  • 普通转录组:DESeq2、edgeR等方法直接比较组间差异

  • 单细胞转录组:需要先分群,再进行细胞类型特异性分析,可使用MAST、Wilcoxon等特殊方法

4. 解析能力不同:

  • 普通转录组:无法区分细胞类型特异性变化

  • 单细胞转录组:可以揭示特定细胞类型的变化,避免少数细胞类型的信号被稀释

5. 批次效应处理不同:

  • 单细胞数据批次效应更为复杂,需要特殊的整合方法(如Harmony、Seurat整合)

二、实践操作:基于Seurat进行单细胞转录组差异基因分析

以下是针对已完成Seurat标准流程并注释好细胞类型的数据,进行不同组间同一细胞类型差异基因分析的完整流程:

我使用的数据是Seurat处理胰腺单细胞转录组数据(panc8数据集),把”celseq2”和“smartseq2”两种技术之间的差异作为分组。所以seurat_obj$group在本次处理过程就是seurat_obj$tech。

1. 准备工作:加载包和数据

# 加载必要的R包
library(Seurat)
library(dplyr)
library(ggplot2)
library(EnhancedVolcano)
# 加载处理好的Seurat对象
seurat_obj <- readRDS("path/to/your_seurat_object.rds")
# 查看数据结构
print(table(seurat_obj$cell_type, seurat_obj$group))

这一步加载了必要的R包和Seurat对象。table函数帮助我们检查每个细胞类型在不同组中的细胞数量,确保有足够的细胞进行比较。

2. 设置差异分析参数

# 设置感兴趣的细胞类型
cell_types_of_interest <- c("acinar")  # 这里仅分析acinar,可扩展到其他类型
# 设置分组变量和比较组
group_var <- "group"  # 元数据中的分组变量名
ident.1 <- "treatment"  # 处理组名称
ident.2 <- "control"    # 对照组名称
# 设置差异分析参数
min.pct <- 0.1  # 至少在一个组中有10%的细胞表达
logfc.threshold <- 0.25  # 对数倍变化阈值
test.use <- "wilcox"  # 使用Wilcoxon秩和检验

这一步设置了分析参数,包括感兴趣的细胞类型、分组变量和比较组,以及差异分析的参数。

3.处理选中的细胞类型

# 创建结果存储列表
deg_results <- list()
#DEGs
Fibroblast <- subset(seurat_obj, celltype=="acinar")
diff_Fibroblast <- FindMarkers(Fibroblast,
                               group.by = group_var,
                               ident.1 = ident.1,
                               ident.2 = ident.2,
                               min.pct = min.pct,
                               logfc.threshold = logfc.threshold,
                               test.use = test.use)
write.csv(diff_Fibroblast, "DEG_acinar_celseq2_vs_smartseq2.csv", row.names = FALSE)

这一步根据前面设置感兴趣的细胞类型,提取子集,进行差异基因分析,并将结果存储在列表中。

4.绘制火山图

sig_genes <- diff_Fibroblast[diff_Fibroblast$p_val_adj < 0.05 & abs(diff_Fibroblast$avg_log2FC) >= 0.5,]
p <- EnhancedVolcano(
  diff_Fibroblast,
  lab = rownames(diff_Fibroblast),
  x = "avg_log2FC",
  y = "p_val_adj",
  pCutoff = 0.05,
  FCcutoff = 0.5,
  pointSize = 2.0,
  labSize = 3.0,
  title = "DEGs in Acinar Cells: celseq2 vs smartseq2",
  subtitle = "Wilcoxon Test",
  caption = paste("Total DEGs:", length(sig_genes)),
  legendPosition = "right",
#selectLab = sig_genes  # 只标记显著基因
)
print(p)
ggsave("volcano_acinar_celseq2_vs_smartseq2.pdf", p, width = 10, height = 8)

这一步为指定的细胞类型绘制火山图,用于可视化差异基因的分布。火山图x轴表示对数倍变化,y轴表示调整后的p值。

图片

5.筛选显著差异基因并进行功能富集分析

# 加载包
library(clusterProfiler)
library(org.Hs.eg.db)
#筛选显著差异基因
#sig_genes <- diff_Fibroblast[diff_Fibroblast$p_val_adj < 0.05 & abs(diff_Fibroblast$avg_log2FC) >= 0.5,]
# 从 sig_genes 中提取基因符号(行名)
sig_gene_symbols <- rownames(sig_genes)
# 将基因符号转换为 Entrez ID
sig_entrez <- bitr(sig_gene_symbols,
                   fromType = "SYMBOL",
                   toType = "ENTREZID",
                   OrgDb = org.Hs.eg.db)
# 查看转换结果
head(sig_entrez)

这一步筛选显著差异基因,因为在绘制火山图之前我已经执行过这一步,这里我注释掉了。可以查看一下sig_genes的情况,提取之后还有5863个基因。

图片

富集分析:

# 4. 进行 GO 富集分析(Biological Process)
go_enrich <- enrichGO(gene = sig_entrez$ENTREZID,
                      OrgDb = org.Hs.eg.db,
                      ont = "BP",  # Biological Process,可改为 "MF" 或 "CC"
                      pAdjustMethod = "BH",  # Benjamini-Hochberg 校正
                      pvalueCutoff = 0.05,
                      qvalueCutoff = 0.2)
# 查看 GO 富集结果前几行
head(go_enrich)
#进行 KEGG 通路富集分析
kegg_enrich <- enrichKEGG(gene = sig_entrez$ENTREZID,
                          organism = "hsa",  # 人类
                          pAdjustMethod = "BH",
                          pvalueCutoff = 0.05,
                          qvalueCutoff = 0.2)
# 查看 KEGG 富集结果前几行
head(kegg_enrich)
# 保存富集分析结果为 CSV 文件
write.csv(as.data.frame(go_enrich), "go_enrichment.csv", row.names = FALSE)
write.csv(as.data.frame(kegg_enrich), "kegg_enrichment.csv", row.names = FALSE)
#  可视化(可选)
# GO 富集结果条形图
barplot(go_enrich, showCategory = 10)  # 显示前 10 个显著类别
ggsave("go_enrichment_barplot.pdf", width = 10, height = 8)
# KEGG 富集结果点图
dotplot(kegg_enrich, showCategory = 10)  # 显示前 10 个显著通路
ggsave("kegg_enrichment_dotplot.pdf", width = 10, height = 8)

图片

图片

注意事项

  1. 确保Seurat对象中的cell_typegroup变量名与自己的数据一致,必要时调整代码中的变量名

  2. 根据数据特点调整参数,如min.pctlogfc.thresholdtest.use

  3. 单细胞数据中,每个细胞类型的细胞数量对结果可靠性有重要影响,建议在分析前检查细胞数量

  4. 考虑使用特定于单细胞数据的差异分析方法,如MAST

  5. 对于稀有细胞类型,可能需要调整参数或使用伪时间分析(pseudobulk)方法

通过以上步骤,可以系统地分析不同组间特定细胞类型的差异基因表达,并通过火山图直观展示结果。这种细胞类型特异性的分析是单细胞转录组技术的独特优势,可以揭示在普通转录组分析中可能被掩盖的细胞类型特异性变化。当然,你也可以选择分析不同类群(亚群)间的差异基因表达。后面富集分析参照普通转录组的做法。

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

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

相关文章

DBA必备脚本:Oracle获取绑定变量的字面SQL文本版版本替代

我们的文章会在微信公众号IT民工的龙马人生和博客网站( www.htz.pw )同步更新 ,欢迎关注收藏,也欢迎大家转载,但是请在文章开始地方标注文章出处,谢谢! 由于博客中有大量代码,通过页面浏览效果更佳。脚本的获取请…

联通光猫烽火吉比特HG6145F获取超级密码

第一步 获取光猫MAC地址 这是两种方法看看光猫实体后面的标注 arp -a 192.168.1.1返回的物理地址第二步 开启Telnet 浏览器输入http://192.168.1.1/telnet?enable=1&key=************ 「************」处为光猫的…

083_尚硅谷_多分支基本使用

083_尚硅谷_多分支基本使用1.分支控制if_else基本语法及说明 2.多分支的流程图 3.多分支案例1 4.多分支案例2

为什么制造业的仓库经验,放到电商就行不通?

我的一些朋友,原先是在制造业做过仓库管理的人,后面转行了电商企业,就直接把原来的一套经验搬过来继续管理仓库,但真正实行后,很快就发现了问题:在制造业里运转得很顺的一套方法,在电商场景下却水土不服。 我个…

Oracle案例:grid环境关于asm diskpath是否需要一致

我们的文章会在微信公众号IT民工的龙马人生和博客网站( www.htz.pw )同步更新 ,欢迎关注收藏,也欢迎大家转载,但是请在文章开始地方标注文章出处,谢谢! 由于博客中有大量代码,通过页面浏览效果更佳。Oracle案例:…

宠物去哪啦小程序系统:智能宠物管理与定位解决方案

一、概述总结 “宠物去哪啦” 是一款聚焦宠物管理与安全的智能小程序系统,支持微信端部署使用。核心依托北斗 + GPS 双模定位的宠物项圈硬件,搭配微擎系统交付的小程序平台,为宠物主人提供实时定位、安全围栏、轨迹…

Windows 如何关闭 dep数据执行保护 - 软件双击没反应的解决办法

1、按 win + s 输入 高级系统设置2、单击“高级”选项卡,然后单击“性能”下面的“设置”;3、单击“数据执行保护”选项卡 选择仅为基本Windows程序和服务启用DEP

2025年整平机厂家推荐排行榜,整平机/校平机/矫平机/开平机/平板机/矫直机/校直机,高精度/精密/液压式/数控/金属/高效稳定/多种规格/全自动整平机公司推荐

2025年整平机厂家推荐排行榜:权威解析行业领军企业 在金属加工行业持续升级的背景下,整平设备作为提升产品质量的关键装备,其技术水平和性能表现直接影响着生产效率和产品精度。随着制造业向高质量方向发展,市场对…

一佳旅游票务系统:旅游行业数字化一体化解决方案

一、概述总结 一佳旅游票务系统是一款聚焦旅游行业的专业化互联网平台,以微信小程序为核心应用载体,提供从景点门票、旅游线路到旅游商品的全链条服务解决方案。系统采用微擎系统交付模式,支持无代码拖拽制作与极速…

2025年10月洗碗机品牌推荐:海信领衔五大机型对比评测榜。

一、引言 对于计划升级厨房、降低餐后劳动强度的家庭用户,以及需要控制后厨人力成本的餐饮创业者而言,洗碗机已从“可选电器”变为“效率刚需”。2025年第四季度,国内能效新标与海外出口补贴同步落地,品牌方集中释…

广告敏感词图文检测微信小程序:高效合规检测解决方案

一、概述总结 广告敏感词图文检测微信小程序是西安立云体网络科技有限公司自主研发的合规检测工具,软著登记号为 2019SR1312226,具备官方正品保障与合法知识产权。该小程序基于微擎系统交付,支持多版本 PHP 环境,提…

2025年10月油烟机品牌推荐:海信领衔静音技术榜对比评测

一、引言 厨房空气管理已成为家庭健康与居住舒适度的关键指标,对于正在装修或焕新厨电的消费者而言,油烟机不仅是排烟工具,更是长期陪伴的能耗设备。用户普遍关注“吸力是否足够、噪声能否更低、清洁是否省心、售后…

Newtonsoft.Json笔记 -JToken、JObject、JArray详解

Newtonsoft.Json(Json.NET)核心概念与高级用法笔记 一、JSON 解析的核心过程 当我们调用: JObject obj = JObject.Parse(jsonString);Json.NET 会执行以下过程:解析 JSON 文本字符串(如 "{ "name"…

软件测试流程-入门

一.测试需求文档 产品需求文档、产品原型图、用户使用手册 重点理解业务需求: 了解熟悉业务,分析需求测试点 二.测试用例 设计测试用例是整个测试中最重要的部分,复杂性也最高。 需要充分理解测试需求和业务流程,才…

什么是人工智能?——AI的定义、发展历程与主要分类

什么是人工智能?——AI的定义、发展历程与主要分类pre { white-space: pre !important; word-wrap: normal !important; overflow-x: auto !important; display: block !important; font-family: "Consolas"…

CF2110F Faculty

这让我想到了信友队的一个题。 首先一个经典结论是,\(x, y\) 其中必有之一为最大值,证明就不证了。 然后你发现如果其他数和 \(maxi\) 的倍数关系不超过 \(2\),这个是很好计算的,如果超过了 \(2\),我们重构一下,…

国产0.38mm超小22pF/50V/C0G电容HLCC2250G,77GHz实测S参数公开,可pin-to-pin替换

国产0.38mm超小22pF/50V/C0G电容HLCC2250G,77GHz实测S参数公开,可pin-to-pin替换各位射频/高速/车载雷达朋友, 恒利泰HLCC2250G,尺寸0.380.380.15 mm,NP0特性,-55 ℃~+125 ℃,22 pF 20 %。 网络仪20 GHz对比Mur…

微信消息管理桌面提醒版:桌面提醒与AI回复的完美结合

在这个信息爆炸的时代,微信已然成为我们日常沟通交流、工作协作以及生活分享不可或缺的工具。每天,大量的微信消息如潮水般涌来,让人应接不暇。错过重要消息、无法及时回复、在海量消息中难以快速定位关键内容等问题…

pip会读取 pyproject.toml 的 project.dependencies 字段进行依赖安装吗?

这是一个非常重要且常见的问题,答案是:通常不会 —— pip 本身不会直接从 pyproject.toml 的 project.dependencies 字段安装依赖,除非你正在安装当前项目本身(例如 pip install . 或 pip install -e .)。详细解释…

昇腾npu架构运行deepseek

Atlas 800-3010 部署 deepseek模型导读 这东西写了有半年了,一直在仓库里吃灰,主要是过程有些不尽人意。 本意是想着能像Windows那种方式,用内存虚拟显存把671B的ds给啃下来,后来发现这条路走不通,设备不支持,生…