单细胞分析(一)——seurat包单个样本处理

2023-10-08 07:40

本文主要是介绍单细胞分析(一)——seurat包单个样本处理,希望对大家解决编程问题提供一定的参考价值,需要的开发者们随着小编来一起学习吧!

seurat包单个样本处理

  • 10X genomics的基本原理
  • 创建对象
    • 加载数据
    • 创建 Seurat 对象
    • count matrix是什么样子?
  • 预处理流程
    • 计算线粒体含量
    • 质控展示
    • 正式筛选
  • 数据标准化
  • 特征筛选
    • 寻找高变基因
    • 高变基因绘图
  • 缩放数据
  • PCA线性降维
    • 降维
    • 降维展示
  • UMAP非线性降维
  • 聚类分群
  • 参考文章

10X genomics的基本原理

大致如下
10X genomics基本步骤在这个教程中,主要将分析 10X Genomics 免费提供的外周血单核细胞 (PBMC) 数据集。在 Illumina NextSeq 500 上对 2,700 个单细胞进行了测序。可以在https://cf.10xgenomics.com/samples/cell/pbmc3k/pbmc3k_filtered_gene_bc_matrices.tar.gz此处找到原始数据。我们从读取数据开始。 Read10X() 函数从 10X 读取 cellranger 管道的输出,返回一个唯一的分子识别 (UMI) 计数矩阵。此矩阵中的值表示在每个单元格(列)中检测到的每个特征(即基因;行)的分子数。

Read10X() 函数是针对于整理好的10X Genomics 数据,如果手头的不是类似文件,可以将其进行转换,成为格式一致的文件。

在这里插入图片描述

接下来使用计数矩阵创建一个 Seurat 对象。该对象用作包含单细胞数据集的数据(如计数矩阵)和分析(如 PCA 或聚类结果)的容器。例如,count matrix存储在 pbmc[[“RNA”]]@counts 中。

library(dplyr)
library(Seurat)
library(patchwork)

创建对象

加载数据

# Load the PBMC dataset
scobj <- Read10X(data.dir = "../data/pbmc3k/filtered_gene_bc_matrices/hg19/")

创建 Seurat 对象

### 2.创建Seurat对象
### counts 输入的是数据,行是基因,列是细胞
### project 参数输入的是项目名称,出现在metadata的orig.ident这一列
### min.cells 限定的是基因:每个基因在至少多少个细胞中出现
### min.features 限定的是细胞: 每个细胞中最少有多少个基因
scobj <- CreateSeuratObject(counts = scdata, project = "pbmc3k", min.cells = 3, min.features = 200)

count matrix是什么样子?

count矩阵是稀松矩阵,可以减少占用空间

pbmc.data[c("IGF2BP2", "TCL1A", "MS4A1"), 1:30]dense.size <- object.size(as.matrix(pbmc.data))
dense.sizesparse.size <- object.size(pbmc.data)
sparse.sizedense.size/sparse.size

预处理流程

计算线粒体含量

这是质控的重要步骤,使用PercentageFeatureSet函数

### 主要PercentageFeatureSet函数计算线粒体含量
### 人类使用pattern = "^MT-",小鼠使用pattern = "^mt-"
scobj[["percent.mt"]] <- PercentageFeatureSet(scobj, pattern = "^MT-")### 该操作会在metadata数据里面增加一列叫做percent.mt
metadata <- scobj@meta.data

一般情况下,可以认为线粒体含量多,意味着细胞可能趋于死亡,这样的细胞就应该剔除。但是如果本身研究的就是和线粒体相关的内容,例如药物干预会引起线粒体的变化,那么就要根据具体情况来分析是否需要剔除。

质控展示

### 质控数据可视化,使用VlnPlot函数
### nFeature_RNA, number of Feature, 每个细胞中有多少个基因
### nCount_RNA, number of counts, 每个细胞中有多少个counts
### percent.mt, 我们自己增加的列,  每个细胞中线粒体基因的比例
VlnPlot(scobj, features = c("nFeature_RNA", "nCount_RNA", "percent.mt"), ncol = 3)

质控展示
大部分情况下,线粒体的含量应该去除较大的,这个可以根据具体的情况进行分析,一般10%以内可以接受。但是依然是要根据相应的干预和特点进行取舍。

也可以展示特征之间的关系

# FeatureScatter is typically used to visualize feature-feature relationships, but can be used
# for anything calculated by the object, i.e. columns in object metadata, PC scores etc.plot1 <- FeatureScatter(scobj, feature1 = "nCount_RNA", feature2 = "percent.mt")
plot2 <- FeatureScatter(scobj, feature1 = "nCount_RNA", feature2 = "nFeature_RNA")
plot1 + plot2

在这里插入图片描述

正式筛选

这一步依然是根据metadata的数据进行筛选,这一部分依然是根据相关的具体情况进行筛选。

### 正式筛选,筛选的是细胞,最终细胞减少
### nFeature_RNA > 200
### nFeature_RNA < 2500
### percent.mt < 5
scobj <- subset(scobj, subset = nFeature_RNA > 200 & nFeature_RNA < 2500 & percent.mt < 5)

筛选之后,再次用小提琴图展示,是上面小提琴图的缩减版本。

至于到底应该确定到多少的筛选指标,可以先按照常规指标进行设定,在后续的细胞聚类中,看看是否有与线粒体相关的群聚类出来,然后再更改筛选指标,聚类指标确实消失,那么就可以对筛选指标确定。

数据标准化

每个细胞分别进行检测,所以如果要比较各个细胞,还是要进行标准化

### 先除以总数,再乘以系数,然后取log
scobj <- NormalizeData(scobj, normalization.method = "LogNormalize", scale.factor = 10000)
### 默认参数可以省略
scobj <- NormalizeData(scobj)

特征筛选

如果某些基因的表达再各个细胞之间表达是恒定的,所以要区分各个细胞,就要使用变化差异大的基因。

寻找高变基因

scobj <- FindVariableFeatures(scobj, selection.method = "vst", nfeatures = 2000)
### 默认参数可以省略
scobj <- FindVariableFeatures(scobj)

高变基因绘图

### 使用VariableFeatures 函数提取高变基因
### 等同于 scobj@assays$RNA@var.features
top10 <- head(VariableFeatures(scobj), 10)
### 使用VariableFeaturePlot 画图
plot1 <- VariableFeaturePlot(scobj)
plot2 <- LabelPoints(plot = plot1, points = top10, repel = TRUE)
plot1 + plot2

缩放数据

缩放的数据达到的效果是:基因的平均值为0,方差为1。便于不同类型之间的比较。

### 降维之前的必备操作
scobj <- ScaleData(scobj, features = rownames(scobj))
### 如果不限定参数,只会缩放高变基因
### scobj <- ScaleData(scobj)
### 缩放后的数据存放在scobj@assays$RNA@scale.data,会很大
scale.data <- scobj@assays$RNA@scale.data

这一步缩放数据会花费比较长的时间,原因是这一步的操作是在所有基因上展开。但是,其实我们最需要的是变异大的基因,也就是可能是标志基因的数据。所以,这个函数默认是可以使用前2000的基因进行缩放的,可以省略features参数的内容。

scobj <- ScaleData(scobj)

需要说明的是,虽然缩放会花费一些时间,数据也会变大,但是,因为后续如果要绘制热图,使用的依然是缩放数据,所以依然在这一步将所有基因的数据进行缩放。
如果,保存数据的时候,还是可以把缩放数据去除后,再进行保存。下次需要的时候再次缩放即可

可以这样操作

scobj@assays$RNA@scale.data <-matrix()

PCA线性降维

PCA操作的对象为缩放过的数据,因此,需要保证前面的缩放数据。结果会保存在对象的reduction中,可以从中提取数据。

为什么要进行主成分分析,主要还是用少量的数据和维度,尽可能的保留原来数据特点。

降维

scobj <- RunPCA(scobj, features = VariableFeatures(object = scobj))
DimPlot(scobj, reduction = "pca")### PCA降维数据存放在scobj@reductions$pca中
pcadata = as.data.frame(scobj@reductions$pca@cell.embeddings)
ggplot(pcadata,aes(PC_1,PC_2,color="red"))+geom_point()### 选择合适的PCA维度
ElbowPlot(scobj)

降维展示

# Examine and visualize PCA results a few different ways
print(scobj[["pca"]], dims = 1:5, nfeatures = 5)
VizDimLoadings(scobj, dims = 1:2, reduction = "pca")


使用热图进行展示

DimHeatmap(scobj, dims = 1, cells = 500, balanced = TRUE)

dims 是确定需要展示维度数,可以尝试

UMAP非线性降维

如何将高维数据压缩到二维平面中,使用UMAP进行处理。其原理可以如下展示:
Understanding UMAP

### 依赖PCA的结果
### dims = 1:10 由上一步确定
scobj <- RunUMAP(scobj, dims = 1:10)
DimPlot(scobj, reduction = "umap")
### UMAP降维数据存放在scobj@reductions$umap中
umapdata = as.data.frame(scobj@reductions$umap@cell.embeddings)
ggplot(umapdata,aes(UMAP_1,UMAP_2,color="red"))+geom_point()

聚类分群

二维空间中看到了细胞分群,那么在这个基础上,怎么进行亚群分类,这就是这一步的意义。

### 找紧邻,dims = 1:10 跟UMAP相同
scobj <- FindNeighbors(scobj, dims = 1:10)
### 分群
scobj <- FindClusters(scobj, resolution = 0.5)
### 会在metadata中增加两列数据"RNA_snn_res.0.5" "seurat_clusters"
metadata <- scobj@meta.data

resolution是这个主要参数,也就是亚群切割的分辨率,这个具体怎么挑选,可以按照以下的方法可视化查看。

### 设置多个resolution选择合适的resolution
scobj <- FindClusters(scobj, resolution = seq(0.2,1.2,0.1))
metadata <- scobj@meta.data
library(clustree)
clustree(scobj)

根据上面的可视化结果,可以大致确定一个分辨率

### 选择特定分辨率得到的分群此处为RNA_snn_res.0.5
scobj@meta.data$seurat_clusters <- scobj@meta.data$RNA_snn_res.0.5
Idents(scobj) <- "seurat_clusters"
DimPlot(scobj, reduction = "umap", label = T)

大致得到了相应的分群,后续就需要对这些分群进行注释。注释是一个与个人知识背景很相关的分析过程,后续补充。

参考文章

Seurat - Guided Clustering Tutorial

Visualization of gene expression with Nebulosa (in Seurat)

Use regularized negative binomial regression to normalize UMI count data

Tutorial: Integrating stimulated vs. control PBMC datasets to learn cell-type specific responses

Understanding UMAP

这篇关于单细胞分析(一)——seurat包单个样本处理的文章就介绍到这儿,希望我们推荐的文章对编程师们有所帮助!



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

相关文章

Python使用getopt处理命令行参数示例解析(最佳实践)

《Python使用getopt处理命令行参数示例解析(最佳实践)》getopt模块是Python标准库中一个简单但强大的命令行参数处理工具,它特别适合那些需要快速实现基本命令行参数解析的场景,或者需要... 目录为什么需要处理命令行参数?getopt模块基础实际应用示例与其他参数处理方式的比较常见问http

Java Response返回值的最佳处理方案

《JavaResponse返回值的最佳处理方案》在开发Web应用程序时,我们经常需要通过HTTP请求从服务器获取响应数据,这些数据可以是JSON、XML、甚至是文件,本篇文章将详细解析Java中处理... 目录摘要概述核心问题:关键技术点:源码解析示例 1:使用HttpURLConnection获取Resp

Java中Switch Case多个条件处理方法举例

《Java中SwitchCase多个条件处理方法举例》Java中switch语句用于根据变量值执行不同代码块,适用于多个条件的处理,:本文主要介绍Java中SwitchCase多个条件处理的相... 目录前言基本语法处理多个条件示例1:合并相同代码的多个case示例2:通过字符串合并多个case进阶用法使用

Java实现优雅日期处理的方案详解

《Java实现优雅日期处理的方案详解》在我们的日常工作中,需要经常处理各种格式,各种类似的的日期或者时间,下面我们就来看看如何使用java处理这样的日期问题吧,感兴趣的小伙伴可以跟随小编一起学习一下... 目录前言一、日期的坑1.1 日期格式化陷阱1.2 时区转换二、优雅方案的进阶之路2.1 线程安全重构2

Java程序进程起来了但是不打印日志的原因分析

《Java程序进程起来了但是不打印日志的原因分析》:本文主要介绍Java程序进程起来了但是不打印日志的原因分析,具有很好的参考价值,希望对大家有所帮助,如有错误或未考虑完全的地方,望不吝赐教... 目录Java程序进程起来了但是不打印日志的原因1、日志配置问题2、日志文件权限问题3、日志文件路径问题4、程序

Python处理函数调用超时的四种方法

《Python处理函数调用超时的四种方法》在实际开发过程中,我们可能会遇到一些场景,需要对函数的执行时间进行限制,例如,当一个函数执行时间过长时,可能会导致程序卡顿、资源占用过高,因此,在某些情况下,... 目录前言func-timeout1. 安装 func-timeout2. 基本用法自定义进程subp

Java字符串操作技巧之语法、示例与应用场景分析

《Java字符串操作技巧之语法、示例与应用场景分析》在Java算法题和日常开发中,字符串处理是必备的核心技能,本文全面梳理Java中字符串的常用操作语法,结合代码示例、应用场景和避坑指南,可快速掌握字... 目录引言1. 基础操作1.1 创建字符串1.2 获取长度1.3 访问字符2. 字符串处理2.1 子字

Java字符串处理全解析(String、StringBuilder与StringBuffer)

《Java字符串处理全解析(String、StringBuilder与StringBuffer)》:本文主要介绍Java字符串处理全解析(String、StringBuilder与StringBu... 目录Java字符串处理全解析:String、StringBuilder与StringBuffer一、St

浅析Java中如何优雅地处理null值

《浅析Java中如何优雅地处理null值》这篇文章主要为大家详细介绍了如何结合Lambda表达式和Optional,让Java更优雅地处理null值,感兴趣的小伙伴可以跟随小编一起学习一下... 目录场景 1:不为 null 则执行场景 2:不为 null 则返回,为 null 则返回特定值或抛出异常场景

深入理解Apache Kafka(分布式流处理平台)

《深入理解ApacheKafka(分布式流处理平台)》ApacheKafka作为现代分布式系统中的核心中间件,为构建高吞吐量、低延迟的数据管道提供了强大支持,本文将深入探讨Kafka的核心概念、架构... 目录引言一、Apache Kafka概述1.1 什么是Kafka?1.2 Kafka的核心概念二、Ka