①单细胞学习-数据读取、降维和分群

目录

①数据读取

②计算线粒体基因比例

③分开进行质控

④两组单细胞数据合并

⑤细胞周期评分

⑥降维标准流程

降维

UMAP可视化

选择分群

⑦marker基因

分析marker基因

marker基因可视化

⑧细胞定群命名

单细胞的数据格式学习:单细胞 10X 和seurat对象学习-CSDN博客

①数据读取

数据集:GSE164241

rm(list = ls()) 
library(Seurat)
folders=list.files('./',pattern='[123]$')
folders
scList = lapply(folders,function(folder){ CreateSeuratObject(counts = Read10X(folder), project = folder,min.cells = 3, min.features = 200)
})BM <- merge(scList[[1]], y = c(scList[[2]],scList[[3]]), add.cell.ids = c("BM1","BM2","BM3"), project = "BM")
BM
GM <- merge(scList[[4]], y = c(scList[[5]],scList[[6]]), add.cell.ids = c("GM1","GM2","GM3"), project = "GM")
GM##两部分数据进行读取合并
②计算线粒体基因比例
  • 细胞基因检出数,低质量细胞基因检出数通常较低,双细胞或者同时捕获多个细胞会有很高的基因数。所以要去除低质量的,和过高的细胞。

  • 细胞检测出的分子数

  • 线粒体基因比例,一般低质量细胞或者死细胞线粒体基因检出数很高。但是特殊情况特殊对待,有些细胞功能活跃,线粒体活跃,检出数自然也会很高。所以不能一刀切。

GM[["percent.mt"]] <- PercentageFeatureSet(GM,pattern = "^MT-")#分别计算线粒体
BM[["percent.mt"]] <- PercentageFeatureSet(BM,pattern = "^MT-")
preQC_GM <- VlnPlot(GM, features = c("nFeature_RNA", "nCount_RNA", "percent.mt"), ncol = 3, group.by = "orig.ident", pt.size = 0)
preQC_BM <- VlnPlot(BM, features = c("nFeature_RNA", "nCount_RNA", "percent.mt"), ncol = 3, group.by = "orig.ident", pt.size = 0)
preQC_GM
preQC_BM
dev.off()


③分开进行质控

nCount_RNA:每个细胞的UMI数量
nFeature_RNA:基因数
percent.mt:线粒体基因百分比

GM、BM分开,涉及后期QC参数有可能不同。

单细胞质控是指对单细胞测序数据进行质量控制和过滤,保证数据的准确性和可靠性。单细胞测序数据存在许多干扰因素,如PCR扩增偏差、杂交、RNA降解等,这些因素可能导致数据出现噪声、低表达、脱落细胞等问题。因此,单细胞质控的主要任务是过滤出质量高、可靠的单细胞数据,以便进行下一步的数据分析和挖掘。

单细胞质控通常包括以下几个方面的内容:

  1. 基础质控:包括测序质量、mapping率、reads数、UMI数等基本指标的计算和分析。
  2. 细胞质控:包括细胞表达量、基因检测率、rRNA比例、mapping率等细胞特异性指标的计算和分析,用于确定是否保留该细胞的数据。
  3. 基因质控:包括基因表达量、基因检测率、基因表达异质性等指标的计算和分析,用于确定保留哪些基因。
  4. 附加质控:包括PCA分析、t-SNE分析、聚类分析等附加分析,用于确定哪些细胞可以被合并或排除。
#设定各指标阈值使用subset函数取子集
#进行质控
GM <- subset(GM, subset = nFeature_RNA > 200 & nFeature_RNA < 5000 & percent.mt < 15)
BM <- subset(BM, subset = nFeature_RNA > 200 & nFeature_RNA < 5000 & percent.mt < 15)
postQC_GM <- VlnPlot(GM, features = c("nFeature_RNA", "nCount_RNA", "percent.mt"), ncol = 3, group.by = "orig.ident", pt.size = 0)
postQC_BM <- VlnPlot(BM, features = c("nFeature_RNA", "nCount_RNA", "percent.mt"), ncol = 3, group.by = "orig.ident", pt.size = 0)
postQC_GM#这时线粒体基因比例都小于15%
postQC_BM

接下来就是将两个数据合并,去除批次效应,整合成一个seurat对象进行下游降维。

④两组单细胞数据合并
#FindIntegrationAnchors合并数据
BM <- NormalizeData(BM)
BM <- FindVariableFeatures(BM, nfeatures = 4000)#数据标准化及计算高变基因
GM <- NormalizeData(GM)
GM <- FindVariableFeatures(GM, nfeatures = 4000)#整合成一个seurat对象进行下游降维,IntegrateData去除批次效应
sampleList <- list(GM, BM)
scedata <- FindIntegrationAnchors(object.list = sampleList, dims = 1:50)
scedata <- IntegrateData(anchorset = scedata, dims = 1:50)
save(scedata, file = "scedata.RData")
⑤细胞周期评分

作者还计算了细胞周期评分,因为我们收集到的细胞可能处于不同的分裂时期,所以看周期是很有必要的,尤其是针对具体的研究目的。

rm(list = ls()) 
library(Seurat)
load("scedata.RData")
s.genes <- cc.genes$s.genes#细胞周期查看(根据需求)
g2m.genes <- cc.genes$g2m.genes
scedata <- CellCycleScoring(scedata, s.features = s.genes, g2m.features = g2m.genes,set.ident = TRUE)
VlnPlot(scedata,features = c("S.Score","G2M.Score"),group.by = "orig.ident")

⑥降维标准流程
降维

单细胞数据RunPCA降维分析前的scale操作,通常指对数据进行标准化或归一化处理,旨在消除不同基因或细胞之间的量纲或数量级的差异。这样做可以使得不同基因或细胞在降维后对PCA的贡献更加均衡,避免了那些具有高表达量的基因或高测序深度的细胞对PCA结果产生过大影响的情况。数据标准化与wilcox分析_wilcox检验的数据必须log2吗-CSDN博客

# 标准流程(vars.to.regress进行降维)
scedata <- ScaleData(scedata, vars.to.regress = c("S.Score", "G2M.Score"), verbose = FALSE)
scedata <- RunPCA(scedata, npcs = 50, verbose = FALSE)
scedata <- FindNeighbors(scedata, reduction = "pca", dims = 1:50)
scedata <- FindClusters(scedata, resolution = seq(from = 0.1, to = 1.0, by = 0.1))
scedata <- RunUMAP(scedata, reduction = "pca", dims = 1:50)library(clustree)
clustree(scedata)

PCA流程结束。我们需要挑选合适的PCs进行后续的细胞亚群鉴定
这里官方给出了两种方法,JackStraw()ElbowPlot()
实际上该选择多少并没有一个明确的规定,往往只能通过继续向下游分群注释去做,出现问题了回来调整或者往下做之后更换几个PC再做一次看结果的重复性是否良好

scedata <- JackStraw(scedata, num.replicate = 100)#JackStraw方法
scedata <- ScoreJackStraw(scedata, dims = 1:20)
JackStrawPlot(scedata, dims = 1:15) #可视化前15个PC
ElbowPlot(scedata)#肘型图
##默认都是计算50个PC可根据需求调整参数

关于单细胞降维最佳PC数量选取 - 简书 (jianshu.com)的筛选。

1.主成分累积贡献大于90%;2.PC本身对方差贡献小于5%;3.两个连续PCs之间差异小于0.1%

UMAP可视化
DimPlot(scedata)#UMAP降维进行可视化

选择分群
#选择好合适的细胞分群,将其设置为active.ident
Idents(scedata) <- "integrated_snn_res.1"#这里根据上面的分析进行选择
scedata$seurat_clusters <- scedata@active.ident
DimPlot(scedata,label = T,split.by = "orig.ident",ncol = 3)#每个样本进行性可视化

⑦marker基因
分析marker基因
#使用FindAllMarkers鉴定各个细胞群的高表达基因
DefaultAssay(scedata) <- "RNA"
all.markers  <- FindAllMarkers(scedata, only.pos = TRUE, min.pct = 0.25, logfc.threshold = 0.75)
significant.markers  <- all.markers [all.markers $p_val_adj < 0.2, ]
#write.csv(significant.markers, file = "significant.markers.csv")#保存
marker基因可视化
markers <- c("ACKR1","RAMP2","SELE","VWF","PECAM1","LUM","COL3A1","DCN","COL1A1","CFD","KRT14","KRT5","S100A2","CSTA","SPRR1B","CD69","CD52","CXCR4","PTPRC","HCST")
DotPlot(scedata,features = markers)+coord_flip()

#UMAP图
FeaturePlot(scedata,features = c("ACKR1","LUM","KRT14","CD69"))
#热图
alldata <- ScaleData(scedata, features = markers, assay = "RNA")
p <- DoHeatmap(alldata, features = markers,group.by = "seurat_clusters",assay = "RNA")
p
dev.off()

⑧细胞定群命名

#对各个细胞群命名:自动定群结果不一定完全正确
scedata <- subset(scedata, idents = c("21"), invert = TRUE)#去掉低质量细胞群
new.cluster.ids <- c("0"="Fibroblast", "1"="Endothelial", "2"="Endothelial", "3"="Endothelial", "4"="Immune", "5"="Immune", "6"="Endothelial", "7"="Fibroblast", "8"="Other", "9"="Immune", "10"="Epithelial", "11"="Endothelial", "12"="Fibroblast", "13"="Immune", "14"="Other", "15"="Immune", "16"="Fibroblast", "17"="Endothelial", "18"="Fibroblast", "19"="Epithelial", "20"="Endothelial", "22"="Immune","23"="Immune", "24"="Immune", "25"="Epithelial","26"="Immune","27"="Immune", "28"="Immune", "29"="Other")
scedata <- RenameIdents(scedata, new.cluster.ids)                        
scedata$celltype <- scedata@active.ident
DimPlot(scedata, group.by = "celltype")
save(scedata, file = "scedata1.RData")##进行降维和清洗后的数据


感谢:TS的美梦

跟着Cell学单细胞转录组分析(八):单细胞转录组差异基因分析及多组结果可视化_单细胞测序差异基因显著性怎么看-CSDN博客

感谢:Biomamba-生信基地的个人空间-Biomamba-生信基地个人主页-哔哩哔哩视频 (bilibili.com)

感谢:03.单细胞数据质控及降维标准流程 - 简书 (jianshu.com)

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

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

相关文章

SpringBoot项目文件上传校验(注解版)

需求 要实现了一个文件上传和验证的功能&#xff0c;具有以下特点&#xff1a; 1. 自定义注解&#xff1a;FileValidation注解用于标记需要进行文件验证的方法。 2. 文件验证拦截器&#xff1a;FileValidationInterceptor拦截器会在每个请求处理之前被调用。如果请求处理的方…

2024年深圳市专精特新企业申报条件-专精特新企业认定、申请时间、流程及奖励补贴

一、深圳专精特新企业申报对象 根据《优质中小企业梯度培育管理暂行办法》&#xff08;工信部企业〔2022〕63号&#xff09;和《深圳市工业和信息化局优质中小企业梯度培育管理实施细则》&#xff08;深工信规〔2022〕7号&#xff09;相关规定&#xff0c;我局组织开展2023年深…

vue2.0滚动加载组件

vue2.0滚动加载组件 一、直接上代码 需求&#xff1a;刚开始用的element-ui的滚动加载组件&#xff0c;个别电脑会在滚动加载没到底就停止了&#xff0c;怀疑是有bug,就自己写了一个 一、直接上代码 <div class"threadListAttach" ref"replyscrollDom"…

git cloen的错误

~ % git clone https://github.com/xxx/core.git Cloning into core... error: RPC failed; curl 56 Recv failure: Operation timed out error: 56 bytes of body are still expected fatal: expected flush after ref listing看起来在克隆仓库时仍然遇到了问题。错误信息显示…

【软件设计师】算法

1、算法的效率 时间复杂度:程序从开始到结束所需要的时间 空间复杂度&#xff1a;算法在运行过程中临时占用存储空间大小的度量 时间渐近复杂度&#xff1a;时间复杂度由最高次幂决定(判断大小技巧&#xff1a;将n10代入&#xff09; O(log2 n):二分查找法 O(n&#xff09;:for…

MySql8.0.25部署MGR集群

1 准备mysql单机实例 当前部署的mysql使用8.0.25&#xff0c;使用传统的方式初始化data目录&#xff0c;启动服务等。 --初始化&#xff0c;start.conf会放在当前文档目录中 ./mysqld --defaults-file/mgr/start.conf --explicit_defaults_for_timestamp --initialize-insecur…

家政预约小程序07服务分类展示

目录 1 创建服务分类页面2 侧边栏选项卡配置3 配置数据列表4 从首页跳转到分类页总结 上一篇我们开发了首页的服务展示功能&#xff0c;本篇我们讲解一下服务分类功能的开发。在小程序中通常在底部导航栏有一个菜单可以展示所有服务&#xff0c;侧边选项卡可以展示分类信息&…

Python零基础一天丝滑入门教程(非常详细)

目录 第1章 初识python 第1节 python介绍 1.为什么要学习Python&#xff1f; 2.python排名 3.python起源 4.python 的设计目标 第2节 软件安装 第2章 快速上手&#xff1a;基础知识 第1节 Python3 基础语法 Python 变量 字面量 数据类型转换 Python3 注释 数据类…

人工智能核心技术:机器学习总览

&#x1f4a1;机器学习作为人工智能的核心&#xff0c;与计算机视觉、自然语言处理、语音处理和知识图谱密切关联 &#x1f4a1;【机器学习】是实现人工智能的核心方法&#xff0c;专门研究计算机如何模拟/实现生物体的学习行为&#xff0c;获取新的知识技能&#xff0c;利用经…

垂类短视频:四川鑫悦里文化传媒有限公司

垂类短视频&#xff1a;内容细分下的新媒体力量 随着移动互联网的迅猛发展和智能手机的普及&#xff0c;短视频已成为当下最受欢迎的媒介形式之一。四川鑫悦里文化传媒有限公司而在短视频领域&#xff0c;一个新兴的概念——“垂类短视频”正逐渐崭露头角&#xff0c;以其独特…

设计模式 21 备忘录模式 Memento Pattern

设计模式 21 备忘录模式 Memento Pattern 1.定义 备忘录模式是一种行为型设计模式&#xff0c;它允许你将一个对象的状态保存到一个独立的“备忘录”对象中&#xff0c;并在之后恢复到该状态。 2.内涵 主要用于以下场景&#xff1a; 需要保存对象状态以备恢复&#xff1a; 当…

torch.matmul()的用法

这篇文章记录torch.matmul()的用法 这里仿照官方文档中的例子说明&#xff0c;此处取整数随机数&#xff0c;用于直观的查看效果&#xff1a; vector x vector 两个一维向量的matmul相当于点积&#xff0c;得到一个标量 tensor1 torch.randint(1, 6, (3,)) tensor2 torch.…

机器学习基础笔记

周志华老师的机器学习初步的笔记 绪论 知识分类 科学 是什么&#xff0c;为什么 技术 怎么做 工程 多快好省 应用 口诀&#xff0c;技巧&#xff0c;实际复杂环境&#xff0c;行行出状元 定义 经典定义 利用经验改善系统自身的性能 训练数据 模型 学习算法 分类 决策树…

Django5+React18前后端分离开发实战14 React-Router6 入门教程

使用nodejs18 首先&#xff0c;将nodejs切换到18版本&#xff1a; nvm use 18创建项目 npm create vitelatest zdpreact_basic_router_dev -- --template react cd zdpreact_basic_router_dev npm install react-router-dom localforage match-sorter sort-by npm run dev此…

nlohmann json C++ 解析

学习材料&#xff1a;nlohmann json json官方 源码解析 源码 要学习并理解这份代码&#xff0c;可以按照以下步骤进行&#xff0c;逐步梳理代码的逻辑&#xff1a; 基本步骤&#xff1a; 配置宏: 理解用于配置的宏定义&#xff0c;这些宏控制库的不同特性和行为。例如&…

Java-常见面试题收集(十五)

二十四 Elasticsearch 1 Elasticsearch 的倒排索引 传统的检索方式是通过文章&#xff0c;逐个遍历找到对应关键词的位置。 倒排索引&#xff0c;是通过分词策略&#xff0c;形成了词和文章的映射关系表&#xff0c;也称倒排表&#xff0c;这种词典 映射表即为倒排索引。 其中…

印度政策变革下,中国跨国企业如何应对?一家高科技企业的数据本地化之路

自2001年底印度加入世贸组织以来&#xff0c;印度政府一直积极采取措施促进经济的发展&#xff0c;推出相关政策吸引外资并调整产业结构&#xff0c;以推动经济实现跨越式增长。外资纷纷涌入印度&#xff0c;在各地建立大规模的企业&#xff0c;促使印度成为全球工厂之一&#…

回答网友问题:在C# 中调用非托管DLL

在一个QQ群里&#xff0c;有人在问如何“在C# 中调用非托管DLL”。 俺脑子抽抽了一下&#xff0c;就回了一句“你喜欢用那种声明方式&#xff0c;就用那种方式去调用。” 然后就有人说&#xff1a;“参数声明要和DLL的声明完全一致”。 俺脑子又抽抽了一下&#xff0c;又回了…

图论中的两种递推计数法

递推计数法 生成树计数&#xff1a; τ ( G ) τ ( G − e ) τ ( G ⋅ e ) \tau(G) \tau(G-e)\tau(G\cdot e) τ(G)τ(G−e)τ(G⋅e) G的生成树的颗数&#xff0c;可以分为两类&#xff1a;包含边e的为 τ ( G ⋅ e ) \tau(G\cdot e) τ(G⋅e)&#xff0c;不包含边e的为 …

kafka跨地区跨集群同步工具MirrorMaker2 —— 筑梦之路

MM2简介 KIP-382: MirrorMaker 2.0 - Apache Kafka - Apache Software Foundation 有四种运行MM2的方法&#xff1a; As a dedicated MirrorMaker cluster.&#xff08;作为专用的MirrorMaker群集&#xff09; As a Connector in a distributed Connect cluster.&#xff08…